O(n log n) technique for PAR2?

[animetosho]

This branches off discussion from this post, where Parfast has seemingly implemented an FFT style technique on PAR2's Vandermonde matrix.

I have no idea of the math and don't know how it works, but if anyone wants to comment or provide info/speculation/discoveries etc, here's a place which won't pollute the existing PAR3 discussion.

[Yutaka-Sawada]

Thanks Ike Devolder and Anime Tosho for discussion place.

PAR2 Encoding is done as multiplication of matrix and vector. (Decoding is multiplication of inverted matrix and vector.)
I search about "vandermonde matrix and fast multiplication" over the internet.

Fast Matrix-Vector Multiplication for Structured Matrices
https://www.netlib.org/utk/people/JackDongarra/etemplates/node384.html
From this article, some kinds of matrices (Vandermonde, Toeplitz, Hankel, and Cauchy) have faster solution.

I found some papars and C source code for Vandermonde matrix. But, it's difficult for me.

The fast transposed Vandermonde solver and its implementation in C
https://summit.sfu.ca/_flysystem/fedora/2024-05/etd23012.pdf

Implementing Kaltofen and Yagati’s fast transposed Vandermonde solver
https://mapletransactions.org/index.php/maple/article/download/21229/17556

Some papars are available for Reed-Solomon Codes, too. Because PAR2's generator matrix may not be full-rank, I'm not sure that fast multiplication works always or fails rarely. Anyway, it seems that PAR2 has a chance to improve speed.

Fast En/Decoding of Reed-Solomon Codes for Failure Recovery
https://ieeexplore.ieee.org/document/9359493

Efficient Erasure Decoding for Generalized Reed Solomon Codes
https://lihaoxu.eng.wayne.edu/NISL/Papers/Tech/GRS-Dec.pdf

New Decoding of Reed-Solomon Codes Based on FFT and Modular Approach
https://arxiv.org/abs/2207.11079

A Fast Decoding Algorithm for Generalized Reed-Solomon Codes and Alternant Codes
https://arxiv.org/pdf/2502.02356

[maksverver]

Thanks for sharing these @Yutaka-Sawada. When reading these papers, it's imporant to check which (finite) field they operate on. Some of them work only in prime fields (GF(p) ) while the ones most useful to us are those that work on fields of characteristic 2 (GF(2^k), like GF(2^16) that is used by PAR 2.0). It may sometimes be possible to adapt an algorithm from one field to another but it's far from trivial.

[maksverver]

Introducting PAR2AFFT

Yesterday I cheated a bit and asked ChatGPT to help me develop an O(n log n) algorithm for multiplication of a transpose Vandermonde matrix by a vector.

The idea was to develop a fast algorithm for multiplying any vector by the entire 65536x65536 transpose Vandermonde matrix defined as VT_ij = j^i for (0 ≤ i, j < 2^16), similar to what Frédéric Didier did in Efficient erasure decoding of Reed-Solomon codes, but using the Additive Fast Fourier Transform (AFFT) to make it run in O(n log n) time.

The AFFT requires converting the input to a different basis that has been called the Cantor basis in Leopard 1.0, and is also called “novel polynomial basis” or LCH basis after the 2014 paper by Lin, Chung, and Han.

I have no idea if what parfast does is similar; the author talked about “Forney transforms” on reddit, but I cannot find anything about how that works, and while Wikipedia does have an article on Forney's algorithm it's not obvious how it applies here, and the AI-generated code is too difficult for me to navigate.

Show me the code!

Back to practical matters: I was able to develop this into a working tool that generates Par2 files with high numbers of blocks much faster than Parpar. The code is here: https://github.com/maksverver/par2afft

The core algorithm was generated by ChatGPT but it's not too complex: gf16_vandermonde_transpose.c

It looks even simpler if I strip away the code used for initialization and testing: core-algorithm.c

This is very little code but it's where 99% of CPU cycles are spent!

It's not hard to see that gf16_vandermonde_transpose_multiply() runs in O(q log q) time, where q is the field size, in our case q=2^16.

The C code is also quite simple, and has lots of room for potentional optimizations, for example using SIMD, or unrolling some of the recursive calls.

Performance

To check how well this runs in practice, I included a benchmark which processes about 40 MiB/s of input data per CPU core. Note that for PAR 2.0 we can use at most half of that since only 32768 constants are used for inputs. By comparison: Leopard 1.0 does about 600 MiB/s. That's 15x faster, but Leopard was heavily optimized with SIMD instructions.

I also did a simple test generating an actual par2 file from 10×100 MiB input files with 32768 input blocks and 32768 recovery blocks using parfast-1.6.0, my new par2afft code, and parpar 0.4.5 for reference. Measuring just CPU time I get the following results (details here:

• parfast 1.6.0: 19.48s (~51 MiB/s)
• par2afft (my tool): 70.04s (~12 MiB/s)
• parpar 0.4.5: 4488.84s (~0.2 MiBs/)

Note that parfast and parpar actually finish faster than this because they use multiple threads, but this is a fairly trivial optimization, so I only measured CPU seconds to keep things fair.

My tool really struggles with I/O because I simply mmap() all input and output files into memory, which puts a lot of pressure on the page cache. It works well on a system with a lot of free memory and a fast SSD. This confirms my suspicion that with fast algorithms like Leopard, a lot of the difficulty of making fast tools shifts from optimizing the CPU logic to optimizing the I/O to make sure the CPU has access to the data it needs without reading the entire input into RAM.

The conclusion is that parfast is really fast, but only within a factor 4x of the AFFT code I wrote in a day. That makes me pretty optimistic about the possibility of optimizing it further to match parfast performance.

At the very least, this should be an explicit demonstration that O(n log n) performance is possible while remaining compatible with the the existing PAR 2.0 file format.

Future work

There are a number of things I may want to explore in the future:

• Optimization: see if the algorithm can be sped up using optimizations like PCLMULQDQ instead of lookup tables for the field arithmetic.

• Decoding: the current code only implements an encoder. Of course a complete tool needs a decoder too. I think a similar approach should be feasible but some details need to be worked out.

• Efficiency at smaller block counts: The current code always performs the full 65536x65536 matrix multiplication which is obviously inefficient when encoding with smaller block counts. An interesting question is whether it's possible to maintain O(n log n) performance when n < q.

[animetosho]

Yes, and that's not even the only O(q log^2 q) approach proposed, but I dismissed these as impractical because for log q = 16, so q^2 = 256 is already quite a big number and the constant factors involved probably would make it infeasible

If it can be based off n instead of q, I imagine a O(n log^2 n) should beat O(q log q) if n is much smaller than q.
In which case, it sounds like an optimised client would implement three different approaches:

• O(n^2) for when n is small
• O(n log^2 n) for when n is larger
• O(q log q) for when n is very large

That's assuming a O(n log^2 n) approach is possible.
Does that sound right, or have I completely missed the mark?

The distinction matters because addition in GF(2^16) is just a XOR operation which is super cheap, while multiplications are expensive

Ah I see.

This looks ugly but can be optimized quite well by the compiler, increasing throughput from ~20 MiB/s to ~100 MiB/s on my computer

Nice!

You might object that the call to subspace_poly_eval() is O(log q) too, but you'll find it amortizes to 4q operations in total, and I added a lookup table to make it run in O(1) time anyway

Actually, if you make additive_fft_transpose iterative, you could solve subspace_poly_eval by taking the value from the previous iteration, and just a gf16_sqr + xor. The squaring implementation already is a lookup.
The other thing is, you'd likely do more than just a 16-bit int at a time, for an optimised implementation, and subspace_poly_eval doesn't need to be repeated when processing more input.

[maksverver]

If it can be based off n instead of q, I imagine a O(n log^2 n) should beat O(q log q) if n is much smaller than q.
In which case, it sounds like an optimised client would implement three different approaches [..]

Yeah, that makes sense to me.

Actually, if you make additive_fft_transpose iterative, you could solve subspace_poly_eval by taking the value from the previous iteration, and just a gf16_sqr + xor. The squaring implementation already is a lookup.

You're right! I hadn't thought of that. I have since eliminated this part entirely by just precalculating the results of those functions for each iteration (there are only 65535 of them, so it doesn't take any more work than the sqr lookup table). Actually the table stores the logarithm of the consant, which saves one more table lookup when doing the multiplication.

But you're correct (again) that the bigger gains at this point are probably in processing multiple elements from each slice during a single transform. I did some experiments multiplying a vector of random elements with a single constant, comparing my current lookup table based GF(2^16) multiplication with an implementation that performs 16 multiplications using in AVX2 parallel using 4 nibble lookup tables, and found that there is a lot of potential for performance improvements there:

Elements:           16      32      64     128     256     512    1024    2048    4096    8192   16384   32768   65536  131072  262144  524288 1048576
Lookup table:    26.48   52.74  101.02  195.96  354.58  550.85  790.72  977.11  728.91  411.67  334.68  600.19  762.15  820.00  826.58  825.97  829.39 MiB/s
Nibble AVX2:     25.17   49.86   98.61  196.93  375.18  721.85 1284.47 2117.21 3200.94 4204.32 5078.58 5649.80 5956.97 5938.60 6171.72 6305.23 6131.87 MiB/s

The lookup table approach tops out at almost 1 GiB/s around 2048 elements (4096 bytes), but the AVX2-accelerated nibble multiplication reaches a peak of 6 GiB/s around 256K elements (512 KiB).

(Also note an interesting effect with the lookup table version where 8192 is much slower than 4096, I'm guessing due to cache effects, which the nibble multiplication avoids because it keeps all its lookup tables in registers.)

Again this shifts the bottleneck to memory: processing e.g. blocks of 32768 elements at a time means keeping 32768 * (32768 + 65536) * 2 bytes = 6 GiB in memory, and double that if you want to overlap I/O with computation. But even at e.g. 2048 elements there is a nice speedup. I'll have to test how well this works in practice.

[animetosho]

Also note an interesting effect with the lookup table version where 8192 is much slower than 4096, I'm guessing due to cache effects, which the nibble multiplication avoids because it keeps all its lookup tables in registers

Your data does show that performance recovers at around 64KB though.
Are you sure it's not just your CPU throttling for part of the benchmark?

The lookup table approach tops out at almost 1 GiB/s around 2048 elements (4096 bytes), but the AVX2-accelerated nibble multiplication reaches a peak of 6 GiB/s around 256K elements (512 KiB).

With the SIMD nibble lookup technique, you generally want to separate low/high bytes into two registers, so that with AVX2 you can do 32 multiplies per iteration (instead of 16).
You can also pre/post-arrange the data so that you only need to de/interleave the bytes once, instead of every time you multiply.

e.g. blocks of 32768 elements at a time means keeping 32768 * (32768 + 65536) * 2 bytes = 6 GiB in memory

If you're reading 64KB chunks, it sounds like it'd be 4GB of memory?
You don't have to read all blocks, then put them into the FFT tree, you can just read them directly into the FFT tree at the appropriate locations.
You'd need a lookup that's a combination of your InputSlice.constant and inverse of node, to direct a given input block to the location in the tree.

and double that if you want to overlap I/O with computation

Doubling sounds rough. I'd probably just buffer a few extra blocks in memory whilst computation is ongoing, instead of trying to construct a second tree at the same time*. But getting good concurrency without blowing up memory usage is a challenge with this scheme.

* I suppose you could just do it to some extent on modern OSes regardless, but relying on overcommit doesn't sit right with me
You could also read data in tree order, instead of sequentially, which might be fine on an SSD.

---

I have been thinking about how I'd implement this technique. The main consideration I had was around memory layout and the order of processing for cache blocking and threading.
Since the fast region multiplication techniques need some constants computed per multiplicand, you ideally want to be multiplying a block of several KB with the same number (to amortize the cost of computing the constants). But you don't want the block to be too large, as that would hurt caching.
65536 is a lot of blocks, so I'd probably operate on 256 lots of 256 blocks instead, for more efficient cache usage.

I also may have multiple smaller trees by splitting up each input block, which might help with threading.

One thing I don't quite get is why processing would get faster when n is larger. I'd expect it to be roughly the same or slightly slower (since you may not have to process the whole tree in monomial_to_novel_transpose).

[maksverver]

Thanks for this comment. Lots of good suggestions for improvement.

Are you sure it's not just your CPU throttling for part of the benchmark?

If you mean thermal throttling, I'm fairly sure that isn't a factor, because the same thing happens when I fix CPU frequency to 2.4 GHz which it can easily maintain especially when loading ony a single core as I do here.

I think there is some cache-related thing going on (maybe the strides getting to wide to fit in L2 cache or something like that) but haven't really checked in detail because it's clear that the avx2 algorithm is much faster at those sizes anyway.

With the SIMD nibble lookup technique, you generally want to separate low/high bytes into two registers, so that with AVX2 you can do 32 multiplies per iteration (instead of 16).

Thanks for the pointer. I'll have to look at it in detail to understand how this works.

My current implementation already separates the low/high result bytes, but duplicates the nibble tables because the shuffle operation treats a single 256 bit register as two separate 128 bit lanes. Probably that can be avoided; I'll have to look into it.

You don't have to read all blocks, then put them into the FFT tree, you can just read them directly into the FFT tree at the appropriate locations.

I think you're right if the algorithm is a sequential:

1. read a chunk of input
2. perform O(q log q) matrix multiplication
3. write a chunk of output
4. goto 1

but I think to get better performance I need to switch to a model where I/O happens in parallel with computation, which implies separate buffers for I/O and computation.

There is no strict requirement that the chunk size used for multiplying is the same as for I/O. This saves a little memory. But in my current implementation the input block count should be close to 32768 for maximum throughput means you either need a lot of memory or do small reads at a time, neither of which is ideal.

One thing I don't quite get is why processing would get faster when n is larger.

Because we're multiplying n elements at a time by the same constant. For the lookup table approach that does essentialy exp[log[a] + log[b]] that means log[b] is the same so you only have to look it up once. For the nibble table it means you only have to load the right nibble table once.

To be clear the benchmarks were for the vector multiplication alone. The matrix multiplication (which uses it) also increases in throughput up as the batch size increases, but the effect is less pronounced and tops at around 8192 elements (16 KiB).

[animetosho]

If you mean thermal throttling, I'm fairly sure that isn't a factor

Oh interesting.

My current implementation already separates the low/high result bytes, but duplicates the nibble tables because the shuffle operation treats a single 256 bit register as two separate 128 bit lanes

Duplicating the tables across the two lanes in AVX2 registers is the typical approach. Though since you're using a lookup table, I wouldn't duplicate them in the LUT - AVX2 does have an efficient way to do a load-broadcast 128->256b (VBROADCASTI128 instruction) (or you could load both 128-bit halves and permute to 2x 256b registers) so there's no need to double memory usage.

I personally wouldn't use a LUT at all as that uses a fair bit of memory on something that isn't that expensive to compute (especially when you're amortizing it across 64KB).
The approach I take is to generate a vector for x*0, x*1 ... x*15 (you can do this via multiplying by 2 and a few xor operations), then do three rounds of multiply by 16. The code might be a bit hard to read, but multiplying by 16 is just ((v << 4) & 0xfff0) ^ reduction[v >> 12] and the nibble table lookup can be used to perform the reduction lookup there.

But what I was referring to earlier though, isn't the lookup tables, but how you perform the lookup operations themselves. You're extracting 4 bits per 16-bit word so you're only processing 25% of the available data at a time.
This can be increased to 50% if you split the low/high bytes, where you'd then extract 4 bits per 8-bit half-word.

but I think to get better performance I need to switch to a model where I/O happens in parallel with computation, which implies separate buffers for I/O and computation.

So earlier when you said "32768 * (32768 + 65536) * 2 bytes", I take it that the "32768 ... * 2 bytes" is the block size (64KB), the 65536 is the tree where all the data goes, and the second 32768 is the input_blocks structure you're using?
Then to overlap I/O with computation, you'd double the memory usage?

If so, what I'm saying is that the input_blocks structure could be removed, and thus the memory usage would be "64KB * 65536" = 4GB, not 6GB.
Then if you want to overlap I/O with computation, you could pre-read, say, 1024 blocks, instead of al 32768, to save memory at the expense of being only a partial overlap. (the overlap factor is easily tunable though)

To be clear the benchmarks were for the vector multiplication alone

Oh that makes sense.

[Yutaka-Sawada]

At the very least, this should be an explicit demonstration that O(n log n) performance is possible while remaining compatible with the the existing PAR 2.0 file format.

Thanks Maks Verver for nice improvement. Using FFT will be a new trend for Reed-Solomon Codes.

Decoding: the current code only implements an encoder. Of course a complete tool needs a decoder too. I think a similar approach should be feasible but some details need to be worked out.

Forney algorithm seems to be used in decoder. Because PAR2 requires Erasure Correction only (no need Error Locator), the decoder would be simpler than classic Reed-Solomon Codes. Anyway, I cannot understand the math.