New code matrix packet types

[maksverver]

Currently, the Par3 spec supports three types of code matrices:

1. Cauchy Matrix Packet
2. Sparse Random Matrix Packet (generated from a specific PRNG)
3. Explicit Matrix Packet

In practice, these require an O(N^3) matrix inversion for recovery, and an O(N^2) matrix multiplication for both generation and recovery. (The sparse random matrix is more efficient during generation; not sure if this also holds for recovery.)

However it looks like the state of the art has improved, with some fast implementations of Reed-Solomon codes being proposed:

• FastECC (O(N log n) RS)
• Leopard-RS (O(N log n) MDS RS)
• Wirehair “O(N) Fountain Code for Large Data”
• Raptor codes

The first two provide the same redundancy as the Cauchy matrix, but can be support generation and recovery in O(N log N) time, making them much faster for large number of blocks. Leopard-RS claims to encode at 450 MB/s and decode at 190 MB/s, which is really fast.

The latter two should be even faster, though like sparse random matrices, they might need slightly more recovery blocks than the number of errors in a cohort.

It might be possible to support some of these coding schemes with an explicit matrix packet (I haven't delved deeply enough into the implementation to figure it out; the math is a bit intimidating), but even if true, that leaves two problems:

1. The explicit matrix packet takes up a lot of space, especially when interleaving is involved (see Missing feature: subslice hashes #11).
2. It may not be trivial to recover an efficient algorithm from an explicit packet.

Which brings me to the main question: should Par3 support any of these coding schemes explicitly? If so, which ones? And if so, should we have have explicit support for interleaving?

It's worth noting that @Yutaka-Sawada already implemented a Leopard matrix packet:

The FFT matrix packet has a type value of "PAR FFT\0" (ASCII). The packet's body contains the following:

*Table: FFT Packet Body Contents*

| Length (bytes) | Type         | Description                          |
|---------------:|:-------------|:-------------------------------------|
|        8       | unsigned int | Index of first input block           |
|        8       | unsigned int | Index of last input block plus 1     |
|        1       | signed int   | Max number of blocks (as power of 2) |
|      0 ~ 8     | byte[]       | Number of interleaving blocks        |

One solution is to just add this to the spec. But that leads to a side question: if Leopard-RS is an option, how useful are the other coding schemes? Are e.g. the sparse random codes sufficiently faster to justify using them instead?

[malaire]

Leopard-RS claims to encode at 450 MB/s and decode at 190 MB/s, which is really fast.

It can do much better with SIMD, see some numbers at reed-solomon-simd which is a Rust implementation based on Leopard-RS.

EDIT: Those numbers use a bit old hardware. Some single-core results on my computer (AMD Ryzen 9 7950X, DDR5-6400):

32 768 source blocks + 32 768 parity blocks

• block size 1 kB: encoding 3.0 GB/s, decoding 1.0 GB/s
• block size 100 kB: encoding 2.4 GB/s, decoding 870-890 MB/s

8 192 source blocks + 57 344 parity blocks

• block size 1 kB: encoding 3.4 GB/s, decoding 1.1 GB/s
• block size 100 kB: encoding 2.6 GB/s, decoding 910 MB/s

57 344 source blocks + 8 192 parity blocks

• block size 1 kB: encoding 2.8 GB/s, decoding 1.0 GB/s
• block size 100 kB: encoding 2.3 GB/s, decoding 840 MB/s

[Yutaka-Sawada]

In practice, these require an O(N^3) matrix inversion for recovery

Cauchy Matrix may not require slow Gaussian Elimination for recovery. You may refer wiki pages Cauchy matrix and Inverse of Cauchy Matrix. Though I cannot understand the mathematical theory, I implemented the method in my par3cmdline sample. This is why I suggested to use Cauchy Matrix instead of other invertible matrices.

Then, I think that current Sparse Random Matrix construction is bad. Because it will require O(N^3) matrix inversion for recovery. The sparse random matrix is made from a specific PRNG. It requires large memory space to keep the matrix, even when only a small part is used. Sparse Random Matrix consumes large memory and is slow to recovery. This is why I didn't implement Sparse Random Matrix in my par3cmdline at this time.

Thus, I suggested to use interleaving method for constructing sparse matrix. The idea seems to be known as "Interleaved Reed–Solomon Codes". The interleaving is used to support many blocks for Leopard-RS in my par3cmdline sample. The same method can be adaptable for Cauchy Matrix, too. I explain the usage below.

Example of creating 4 parity blocks from 10 source blocks as normal RS Codes.
[s0][s1][s2][s3][s4][s5][s6][s7][s8][s9] -> [p0]
[s0][s1][s2][s3][s4][s5][s6][s7][s8][s9] -> [p1]
[s0][s1][s2][s3][s4][s5][s6][s7][s8][s9] -> [p2]
[s0][s1][s2][s3][s4][s5][s6][s7][s8][s9] -> [p4]
Each parity block can recover one lost source block.

Now, this is an example of making a sparce matrix of 50% density.
Create 4 parity blocks from 10 source blocks as Interleaved RS Codes.
[s0][__][s2][__][s4][__][s6][__][s8][__] -> [p0]
[__][s1][__][s3][__][s5][__][s7][__][s9] -> [p1]
[s0][__][s2][__][s4][__][s6][__][s8][__] -> [p2]
[__][s1][__][s3][__][s5][__][s7][__][s9] -> [p3]

The sparce matrix is made form independent 2 matrices.
Create 2 parity blocks from 5 source blocks in each matrix.
[s0][s2][s4][s6][s8] -> [p0]
[s0][s2][s4][s6][s8] -> [p2]
and
[s1][s3][s5][s7][s9] -> [p1]
[s1][s3][s5][s7][s9] -> [p3]
Each parity block can recover one lost source block in their matrix.

I call each sub-matrix as cohort. Each cohort is made by Cauchy Matrix. When sparce matrix is 50% density, there are 2 cohorts. Creation and Recovery can be done in each cohort independently. It's same as using 2 Cauchy Matrices to construct 1 sparce matrix. The speed is 2 times faster than normal dense matrix. By setting more cohorts (interleave many times), it becomes more sparce. As it becomes more sparce, it will be faster. And it can support more blocks with 16-bit Reed–Solomon Codes, too. I think that "Interleaved Reed–Solomon Codes" is better than current Sparse Random Matrix.

[maksverver]

Cauchy Matrix may not require slow Gaussian Elimination for recovery.

Thanks for clarifying. I did not know this was possible! But even so, it seems like both generating and recovering using the Gauchy matrix is still a lot slower than Leopard, both in theory and practice. So the general question stands: why would I want to use a Gauchy matrix when Leopard is an option?

Thus, I suggested to use interleaving method for constructing sparse matrix.

Right. To be clear, this works for any coding matrix, including the sparse random matrix, because it is essentially equivalent to subslicing. I think that's more obvious if you write it from top to bottom, where each cohort is a column:

[s0][s1]
[s2][s3]
[s4][s5]
[s6][s7]
[s8][s9]
--------
[p0][p1]
[p2][p3]

Ideally all supported matrix types should be able to use interleaving.

[Yutaka-Sawada]

why would I want to use a Gauchy matrix when Leopard is an option?

Reed-Solomon Codes over Cauchy Matrix is simple and flexible. FFT-based Reed-Solomon Codes is fast for many blocks. Their properties are different. Even though Leopard-RS is good for most usage, it may not be useful in rare case.

1. Leopard-RS needs to load all source blocks to calculate parity blocks.

When Leopard-RS calculates parity blocks, it needs to read all source blocks at first. This isn't a problem for small data or SSD. When file data is larger than RAM size, it causes random access on data storage. It may be a problem for disc drives like HDD. Leopard-RS isn't suitable for streaming input or incremental reading.

2. Leopard-RS calculates all parity blocks at once.

You may create additional parity blocks after the first creation. For example, you made 1000 parity blocks. If you want to make 10 more parity blocks later, Leopard-RS calculates whole 1010 parity blocks and saves the last 10 blocks. It may not be so fast to make a few blocks, as it just discards preceding 1000 parity blocks. Also, Leopard-RS sets max number of blocks at first, and it cannot exceed the limit in next time. When you plan to add more parity blocks in future, Leopard-RS may not be useful.

[maksverver]

Reed-Solomon Codes over Cauchy Matrix is simple and flexible.

Using the Cauchy matrix is simpler by itself, but if Par3 implementations must support both Cauchy and Leopard encoding, that's more complex than supporting Leopard alone. And if they only support Cauchy, then they miss out on large performance improvements.

1. Leopard-RS needs to load all source blocks to calculate parity blocks. [..] Leopard-RS isn't suitable for streaming input or incremental reading.

This is a fair point, but the counterpoint is that if you generate Cauchy (or sparse random) recovery bocks in a streaming fashion, generating only one recovery slice at a time, then you must stream the same input files over and over, once for each recovery block generated, which may well end up being a lot slower.

I don't think the memory requirements for Leopard-RS are prohibitive. If I understood correctly, the library can work on buffers as small as 64 bytes per piece. Using the maximum 65536 pieces, that's only 64 × 64 KiB = 4 MiB total memory needed.

If instead we process 4 KiB at a time (the default Linux page size), that would require 256 MiB in RAM at a time with the maximum number of 65536 pieces. That's not a tiny amount of memory, but also not that much for a modern system.

2. Leopard-RS calculates all parity blocks at once. [..] For example, you made 1000 parity blocks. If you want to make 10 more parity blocks later, Leopard-RS calculates whole 1010 parity blocks and saves the last 10 blocks.

It's not clear to me if that's an intrinsic limitation of the encoding scheme, or just a limitation of the library API, but even assuming the former: Leopard's time complexity scales with O(N log M) where N and M are the number of input and output blocks respectively. Generating 1010 output blocks would take log(1010)/log(10) = ~3 times as long as generating the first 10 blocks alone. It's slower, but not terribly slow. It may even still be faster than using the Cauchy matrix: log(1010)/log(2) = 9.98 = ~10, so which one is faster in this case depends on constant factors, and generating the first 1000 blocks with a Gauchy matrix would definitely be slower!

So I don't think these counter-examples are a strong argument against using Leopard-RS for everything. While conversely, the Gauchy matrix is a lot slower in very common cases, such as generating a large number of recovery blocks, or recovering from large numbers of errors.

[malaire]

If I understood correctly, the library can work on buffers as small as 64 bytes per piece

That is for SIMD optimizations, algorithm itself can work down to 2 bytes per piece.

Leopard's time complexity scales with O(N log M) where N and M are the number of input and output blocks respectively

I believe it's O(n log k) for (n,k) erasure code, i.e. n is total number of blocks (input + output) and k is number of input blocks.

[maksverver]

That is for SIMD optimizations, algorithm itself can work down to 2 bytes per piece.

Fair enough, but presumably performance suffers when not using SIMD instructions (and I'm not sure the linked library supports this, though perhaps your rust crate does). My point was that we can easily do 64 or even 4096 bytes at a time, so we can achieve great performance without using an enormous amount of memory.

I believe it's O(n log k) for (n,k) erasure code, i.e. n is total number of blocks (input + output) and k is number of input blocks.

I think you might be right. I got confused by the Leopard README which says: “Operations are done O(K Log M), where K is the original data size, and M is up to twice the size of the recovery set.”, which seems wrong.

The conclusion that generating additional blocks is pretty fast still stands; if anything, the math is more favorable: Let's say we have 2000 input blocks. The cost of generating 10, 1000 and 1010 recovery blocks is 2010 log 2000, 3000 log 2000, 3010 log 2000 respectively; the log 2000 factors are the same, so generating 10 additional blocks after the fact is only about 50% more expensive than generating those 10 blocks up front would have been (3010/2010 = ~1.4975).

[animetosho]

When you're doing 4KB at a time, does this mean the I/O pattern is to read 4KB, seek to the next block and read the next 4KB etc? If so, 4KB is awfully tiny for an I/O buffer, especially if it's non-sequential.

I recall PAR3 was enabling more than GF16, which would mean 65536 is not the block limit, so that's just a limitation of Leopard?

but the counterpoint is that if you generate Cauchy (or sparse random) recovery bocks in a streaming fashion, generating only one recovery slice at a time

I think 'streaming' here refers to streaming the input, not recovery blocks, in which case, you'd be generating all recovery at once, not one slice at a time. 'Streaming' the recovery can be done (in the sense of writing to a pipe/stream), just that you have to buffer all the output slices in memory.

[malaire]

I recall PAR3 was enabling more than GF16, which would mean 65536 is not the block limit, so that's just a limitation of Leopard?

Yes that's current limitation of Leopard. Author mentions some challenges of increasing that to 2^32 in Request: support >65536 pieces (2017), e.g. that tables used in algorithm could then take several GBs of RAM.

[malaire]

When you're doing 4KB at a time, does this mean the I/O pattern is to read 4KB, seek to the next block and read the next 4KB etc? If so, 4KB is awfully tiny for an I/O buffer, especially if it's non-sequential.

Algorithm requires corresponding pieces from all source blocks in RAM, so yes it needs to do non-sequential reads across all source blocks. With 4 KB buffer algorithm starts by reading first 4 KB of each source block, uses those to calculates first 4 KB of each parity block and writes those to disk, then reads next 4 KB of each source block and so on.

So the problem is that selected buffer size will be multiplied by block count (upto 65536), so using e.g. 32 KB buffer would require upto 32 KB * 65536 = 2 GB total buffer size.

So the algorithm is OK for SSD or with a lot of RAM, but with non-SSD and little RAM it can get slow because of all the seeking.

I believe this is same for all implementations, not just Leopard.

p.s. This is also one challenge in increasing block count beyond 2^16.

[maksverver]

When you're doing 4KB at a time, does this mean the I/O pattern is to read 4KB, seek to the next block and read the next 4KB etc? If so, 4KB is awfully tiny for an I/O buffer, especially if it's non-sequential.

Yes, that was the idea, and I agree the access pattern is pretty bad. But the question is: how does it compare to the alternative?

A quick test with par3 suggests it processes around 2 MB/s of input when encoding 5000 input blocks + 1000 recovery blocks with the Gauchy matrix, but probably this can be optimized/parallelized a lot.

Ignoring caches/readahead entirely, a hard disk that can do 100 IOPS would do only 0.4 MB/s with a 4096 random read speed. So you'd need to increase the buffer size to 20 KiB just to break even.

I see on my system the readahead is 128 KiB. Using that as the buffer size would bring the speed up to 12.5 MB/s taking 750 MB, which is a lot faster than Gauchy but still pretty slow (orders of magnitudes away from the rate that Leopard could sustain, so we're purely I/O bound now). And at the maximum 65536 pieces in a cohort this would take 8 GiB of RAM.

But notice that this is assuming an old fashioned HDD; I bet SSDs are actually a lot faster even at relatively small random reads.

I think 'streaming' here refers to streaming the input, not recovery blocks, in which case, you'd be generating all recovery at once, not one slice at a time.

Fair, but then you still need enough RAM to buffer all output slices at least. Admittedly, that's usually a fraction of the input size (say, 10%) which saves a considerable amount of RAM.

OK, I'm starting to become convinced of the potential advantages of the Cauchy encoding scheme vs Leopard-RS. One is CPU-bound but relatively efficient when it comes to I/O, the other is very fast CPU-wise but requires awkward I/O which may be slow.

[malaire]

But notice that this is assuming an old fashioned HDD; I bet SSDs are actually a lot faster even at relatively small random reads.

SSDs are benchmarked by manufacturers using 4KB random reads. Even older SATA SSDs can sustain 10k - 100k IOPS with 4 KB reads, and newer M.2 SSDs can reach over 1M IOPS. So yeah, this is non-issue for SSDs.

[animetosho]

Author mentions some challenges of increasing that to 2^32 in Request: support >65536 pieces (2017), e.g. that tables used in algorithm could then take several GBs of RAM

That's interesting. I don't understand the math behind all of it, but I'll point out that GF32 multiplication absolutely can be done using much smaller lookup tables, just that you may need to do several lookups per multiply.
All the fast GF region multiplication techniques use SIMD instead of lookups though, so it's rather moot. For single multiplies, all modern CPUs have 64-bit carryless multiply instructions (PCLMULQDQ on x86, PMULL on ARM).

I believe this is same for all implementations, not just Leopard.

For PAR2, implementations typically buffer recovery slices instead of input blocks, since you often have less recovery data. But yes, if there's not enough memory to buffer all recovery, you'll need to work with it in chunks, but you can choose not to do that and just run multiple passes with fewer recovery slices per pass.
Ignoring I/O, too small chunks actually impacts CPU as well, due to the overhead for setting up SIMD constants for a region multiply, so doing multiple passes over the data with larger chunks may be beneficial.

A quick test with par3 suggests it processes around 2 MB/s of input when encoding 5000 input blocks + 1000 recovery blocks with the Gauchy matrix, but probably this can be optimized/parallelized a lot.

You can test optimised PAR2 implementations (such as MultiPar and ParPar) as a rough guide on what GF16, O(N*M) can do.

[Yutaka-Sawada]

I tested Leopard-RS, Wirehair, and nanorq (RaptorQ RFC 6330 library for gcc) on my PC. I made similar sample applications and compare speed.

1. Leopard-RS (Reed-Solomon Codes with FFT)

Leopard-RS library can use OpenMP for multi-threading. Because I don't enable OpenMP in my par3cmdline sample, it uses single thread only now. When I enables OpenMP for tests, it seems to be slower. I don't know why it can become slow by using more threads. It may happen to be busy by using too many threads. Because it's difficult to determine good number of using threads, I disable OpenMP at this time.

2. Wirehair

Wirehair uses Fountain Codes (a kind of Tornade Codes). Wirehair seems to be faster than Leopard-RS at decoding. From the nature of Reed-Solomon Codes, its decoding complexity is higher than encoding. At encoding, Leopard-RS and Wirehair are similar level. One is faster in a case, and another is faster in another case. The most problem of using Wirehair is patent issue. Though basic patent of Fountain Codes might expire, I'm not sure others.

3. RaptorQ (Rapid Tornade Codes for Qualcomm ?)

I uses nanorq as RaptorQ implementation. Because it was made for gcc, I modified it for MS Visual Studio. Though RFC 6330 sets interleaving by default, I disabled the feature in test. RaptorQ supports max 56403 source blocks (source symbols in RFC 6330). While nanorq is a little slower than Wirehair, the property is same. More optimization may improve speed of RaptorQ in nanorq.

Caution ! RFC 6330 words are different from PAR3 definition.

RFC 6330         PAR3
source symbol  = source block
encoded symbol = recovery block
symbol size    = block size
source block   = interleaved cohort

Some test results on Intel Core i5 10400

[ 32000 source blocks, 8000 parity blocks, block size = 8000 bytes ]

Leopard Encoder(256 MB in 32000 pieces, 8000 parities): Input= 1080.17 MB/s, Output= 270.042 MB/s
Leopard Decoder(256 MB in 32000 pieces, 8000 losses): Input= 317.618 MB/s, Output= 79.4045 MB/s

Wirehair Encoder(256 MB in 32000 pieces, 8000 parities): Input= 679.045 MB/s, Output= 169.761 MB/s
Wirehair Decoder(256 MB in 32000 pieces, 8000 losses): Input= 519.27 MB/s, Output= 129.817 MB/s

nanorq Encoder(256 MB in 32000 pieces, 8000 parities): Input= 368.876 MB/s, Output= 92.219 MB/s
nanorq Decoder(256 MB in 32000 pieces, 8000 losses): Input= 382.09 MB/s, Output= 95.5224 MB/s

[ 32000 source blocks, 16000 parity blocks, block size = 8000 bytes ]

Leopard Encoder(256 MB in 32000 pieces, 16000 parities): Input= 879.725 MB/s, Output= 439.863 MB/s
Leopard Decoder(256 MB in 32000 pieces, 16000 losses): Input= 303.318 MB/s, Output= 151.659 MB/s

Wirehair Encoder(256 MB in 32000 pieces, 16000 parities): Input= 622.871 MB/s, Output= 311.436 MB/s
Wirehair Decoder(256 MB in 32000 pieces, 16000 losses): Input= 486.692 MB/s, Output= 243.346 MB/s

nanorq Encoder(256 MB in 32000 pieces, 16000 parities): Input= 345.479 MB/s, Output= 172.74 MB/s
nanorq Decoder(256 MB in 32000 pieces, 16000 losses): Input= 347.354 MB/s, Output= 173.677 MB/s

[ 32000 source blocks, 32000 parity blocks, block size = 8000 bytes ]

Leopard Encoder(256 MB in 32000 pieces, 32000 parities): Input= 648.101 MB/s, Output= 648.101 MB/s
Leopard Decoder(256 MB in 32000 pieces, 32000 losses): Input= 276.757 MB/s, Output= 276.757 MB/s

Wirehair Encoder(256 MB in 32000 pieces, 32000 parities): Input= 530.021 MB/s, Output= 530.021 MB/s
Wirehair Decoder(256 MB in 32000 pieces, 32000 losses): Input= 455.516 MB/s, Output= 455.516 MB/s

nanorq Encoder(256 MB in 32000 pieces, 32000 parities): Input= 318.012 MB/s, Output= 318.012 MB/s
nanorq Decoder(256 MB in 32000 pieces, 32000 losses): Input= 303.318 MB/s, Output= 303.318 MB/s

[ 20000 source blocks, 20000 parity blocks, block size = 16000 bytes ]

Leopard Encoder(320 MB in 20000 pieces, 20000 parities): Input= 547.009 MB/s, Output= 547.009 MB/s
Leopard Decoder(320 MB in 20000 pieces, 20000 losses): Input= 191.962 MB/s, Output= 191.962 MB/s

Wirehair Encoder(320 MB in 20000 pieces, 20000 parities): Input= 560.42 MB/s, Output= 560.42 MB/s
Wirehair Decoder(320 MB in 20000 pieces, 20000 losses): Input= 481.203 MB/s, Output= 481.203 MB/s

nanorq Encoder(320 MB in 20000 pieces, 20000 parities): Input= 341.515 MB/s, Output= 341.515 MB/s
nanorq Decoder(320 MB in 20000 pieces, 20000 losses): Input= 331.606 MB/s, Output= 331.606 MB/s

[ 10000 source blocks, 10000 parity blocks, block size = 32000 bytes ]

Leopard Encoder(320 MB in 10000 pieces, 10000 parities): Input= 598.131 MB/s, Output= 598.131 MB/s
Leopard Decoder(320 MB in 10000 pieces, 10000 losses): Input= 197.044 MB/s, Output= 197.044 MB/s

Wirehair Encoder(320 MB in 10000 pieces, 10000 parities): Input= 588.235 MB/s, Output= 588.235 MB/s
Wirehair Decoder(320 MB in 10000 pieces, 10000 losses): Input= 500 MB/s, Output= 500 MB/s

nanorq Encoder(320 MB in 10000 pieces, 10000 parities): Input= 354.374 MB/s, Output= 354.374 MB/s
nanorq Decoder(320 MB in 10000 pieces, 10000 losses): Input= 338.624 MB/s, Output= 338.624 MB/s

[ 5000 source blocks, 5000 parity blocks, block size = 64000 bytes ]

Leopard Encoder(320 MB in 5000 pieces, 5000 parities): Input= 613.027 MB/s, Output= 613.027 MB/s
Leopard Decoder(320 MB in 5000 pieces, 5000 losses): Input= 205.656 MB/s, Output= 205.656 MB/s

Wirehair Encoder(320 MB in 5000 pieces, 5000 parities): Input= 570.41 MB/s, Output= 570.41 MB/s
Wirehair Decoder(320 MB in 5000 pieces, 5000 losses): Input= 494.59 MB/s, Output= 494.59 MB/s

nanorq Encoder(320 MB in 5000 pieces, 5000 parities): Input= 356.745 MB/s, Output= 356.745 MB/s
nanorq Decoder(320 MB in 5000 pieces, 5000 losses): Input= 336.488 MB/s, Output= 336.488 MB/s

Does someone know current patent issues of Wirehair or RaptorQ ? Qualcomm seems to have patents about RaptorQ. I found a post "What is the patent situation with Raptor codes nowadays?". There is "Legal status" item on wikipedia's Raptor code page. That may be "QUALCOMM Incorporated's Statement about IPR related to RFC 6330". Legal words are too difficult for my English skill. I don't know that adding Wirehair or RaptorQ to PAR3 is good or not.

If someone knows other fast recovery codes or useful library, post here. I will test it.