Reed-Solomon Codes

MDS codes

An MDS (Maximum Distance Separable) code is a linear code of dimension $k$ and length $n$ over a finite field $\mathbb{F}$ (a vector subspace of dimension $k$ of ${\mathbb{F}}^n$), which reaches the Singleton bound (the minimal Hamming distance between any two codewords is $d=n-k+1$), and so provides the maximum erasure-correcting capability possible (as a linear code cannot verify $d>n-k+1$, again by the Singleton bound). So we can correct up to $d-1=n-k$ erasures by finding the closest codeword under the Hamming distance. A generating matrix $G\in {\mathbb{F}}^{n\times k}$ of a linear code $\mathcal{C}$ is defined by $\mathcal{C}=\text{Column-span}(G)=\{G{\mathbf{x}}^T:\mathbf{x}\in {\mathbb{F}}^k\}$. The property that any $k$ coefficients $\tilde{\mathbf{c}}$ of a codeword $\mathbf{c}$ determine it, comes from the fact that any set of $k$ rows of its generating matrix forms a full rank $k\times k$ matrix $A$ (the equation $A{\mathbf{c}}^T={\tilde{\mathbf{c}}}^T$ where $\mathbf{c}$ is the unknown has a unique solution, as $A$ is invertible).

Reed-Solomon codes

Let $\mathbb{F}$ be a prime field of order $r$. Let $n=\alpha k$ such that $n|r-1$. Let $\omega \in \mathbb{F}$ be a primitive $n$-th root of unity. Aside: how do we find such primitive roots of unity in practice?

We consider the Reed-Solomon code of parameters $[n,k,d=n-k+1]$ and evaluation points, $\mathbf{\omega }=(1,\omega ,{\omega }^2,\dots ,{\omega }^{n-1})$ the subgroup of the $n$-th roots of unity. We denote it $\text{RS}(n,k,\mathbf{\omega }):=\{(f({\omega }^i){)}_{i\in [[0,n-1]]}|f\in \mathbb{F}[x]\wedge \mathrm{deg}f<k\}$.

RS has rate $R=k/n=1/\alpha$. To a vector $\mathbf{f}$ we associate the polynomial $f(x)={\sum }_i{f}_i{x}^i$, and reciprocally. RS is a linear code, whose generating matrix is the following Vandermonde matrix, for which every square submatrix is invertible, so RS is indeed MDS:

$$\left[\begin{array}{ccccc}1& {\omega }_1& {\omega }_1^2& \dots & {\omega }_1^{k-1}\\ 1& {\omega }_2& {\omega }_2^2& \dots & {\omega }_2^{k-1}\\ ⋮& & & \ddots & ⋮\\ 1& {\omega }_n& {\omega }_n^2& \dots & {\omega }_n^{k-1}\end{array}\right].$$

Encoding

Encoding a message $\mathbf{m}=(m_0,\dots ,m_{k-1})\in {\mathbb{F}}^k$ amounts to evaluate its associated polynomial $m(x)={\sum }_{i=0}^{k-1}{m}_i{x}^i$ at the evaluation points $\mathbf{\omega }$. This can be done with an $n$-points discrete Fourier transform supported by $\mathbb{F}$ in time $\mathcal{O}(n\log n)$:

Input: $\mathbf{m}=(m_0,\dots ,m_{k-1})\in {\mathbb{F}}^k$

Output: $\mathbf{c}=(c_0,\dots ,c_{n-1})\in \text{RS}(n,k,\mathbf{\omega })$

Return ${\text{FFT}}_n({\text{IFFT}}_k(\mathbf{m})\Vert 0_{{\mathbb{F}}^{n-k}})$

Decoding from erasures

As we saw earlier, we can decode a codeword $\mathbf{c}\in \text{RS}(n,k,\mathbf{\omega })$ with at most $d-1=n-k$ erasures, i.e. from at least $k$ components of $\mathbf{c}$.

Without loss of generality, let $\tilde{\mathbf{c}}=(c_0,\dots ,c_{k-1})$ be the received codeword with erasures, the first $k$ components of a codeword. We can retrieve the original message $\mathbf{m}$ with any $k$ components of $\mathbf{c}$ thanks to the Lagrange interpolation polynomial, where the $x_i$ are the evaluation points of $\tilde{\mathbf{c}}$

$m(x)={\sum }_{i=0}^{k-1}{c}_i{\prod }_{\begin{array}{c}j=0\\ j\ne i\end{array}}^{k-1}\frac{x-x_j}{x_i-x_j}.$

As detailed in https://arxiv.org/pdf/0907.1788v1.pdf, the idea is to rewrite $m(x)$ as a product of two polynomials $A(x)$ and $B(x)$ so that the convolution theorem allows us to recover $m(x)$ using FFTs in $O(n\log n)$. (The authors consider sums to $n-1$ while we can consider sums to $k-1$ since $m(x)$ has degree $k$.)

To do so, let $A(x):={\prod }_{i=0}^{k-1}(x-x_i),A_i(x):={\prod }_{\begin{array}{c}j=0\\ j\ne i\end{array}}^{k-1}(x-x_j).$

Let $n_i:=\frac{c_i}{A_i(x_i)}$.

The interpolation polynomial becomes: $m(x)=A(x){\sum }_{i=0}^{k-1}\frac{c_i}{(x-x_i)A_i(x_i)}=A(x){\sum }_{i=0}^{k-1}\frac{n_i}{x-x_i}.$

Note that $A_i(x_i)\ne 0$ by definition, so it is invertible in $\mathbb{F}$.

In order to replace the costly product $A_i(x)$ in this expression, we use the fact that the formal derivative $A^{\prime }(x)$ of $A(x)$ satisfies for all $i\in [[0,k-1]]$: $A^{\prime }(x_i)=A_i(x_i)$. So we can compute $(A_i(x_i){)}_i$ by evaluating $A^{\prime }(x)$ at the points $\mathbf{\omega }$ with an FFT.

Indeed: $A^{\prime }(x)=({\prod }_{i=0}^{k-1}(x-x_i){)}^{\prime }={\sum }_{i=0}^{k-1}(x-x_i{)}^{\prime }{\prod }_{\begin{array}{c}j=0\\ j\ne i\end{array}}^{k-1}(x-x_j)={\sum }_{j=0}^{k-1}{A}_j(x).$ So $A^{\prime }(x_i)={\sum }_{j=0}^{k-1}{A}_j(x_i)=A_i(x_i)$ as the other polynomials $A_j(x)$ have $x_i$ as root.

The resulting expression corresponds to the barycentric formula.

Writing the fraction $\frac{1}{x_i-x}={\sum }_{j=0}^{\infty }\frac{x^j}{x_i^{j+1}}$ as a formal power series, we obtain

$$\begin{array}{rl}m(x)/A(x)=\sum _{i=0}^{k-1}\frac{n_i}{x-x_i}\mathrm{mod}{x}^k& =-\sum _{i=0}^{k-1}\left(\sum _{j=0}^{k-1}\frac{n_i}{x_i^{j+1}}{x}^j\right).\end{array}$$

But if we let $N(x):={\sum }_{i=0}^{k-1}\frac{n_i}{x_i}{x}^i$ then ${\sum }_{i=0}^{k-1}\frac{n_i}{x-x_i}\mathrm{mod}{x}^k=-{\sum }_{j=0}^{k-1}N(x_i^{-j})x^j=:-B(x).$

$B(x)$ is thus given by the first $k$ components of $n\times {\text{IFFT}}_n(\mathbf{N})$.

So the product is given by the convolution theorem, and by linearity of the DFT and pairwise product:

$$\begin{array}{rl}\mathbf{m}& =\mathbf{A}\ast (-\mathbf{B})=-{\text{IFFT}}_{2k}({\text{FFT}}_{2k}(\mathbf{A})\odot {\text{FFT}}_{2k}(\mathbf{B})).\end{array}$$

The total cost is $O(k{\log }^{2}k+n\log n)$: the first term accounts for the computation of the product for $A(x)$ with a divide and conquer approach for the multiplication of its factors with FFT multiplication, and the second one for the other steps.

Sharding

In some applications, such as database sharding, it can be advantageous to split the codewords into chunks (aka shards). For this purpose, let $s$ be the number of shards, $l=n/s$ the length of a shard, and $\omega$ a primitive $n$-th root of unity.

The domain of evaluation is then split into cosets: $⟨\omega ⟩={⨆}_{i\in [[0,s-1]]}{\Omega }_i$, for ${\Omega }_0=\{{\omega }^{sj}{\}}_{j\in [[0,l-1]]}$ and ${\Omega }_i={\omega }^i{\Omega }_0$.

For a set of $k/s$ shard indices $Z\subseteq \{0,s-1\}$, we reorganize the product $A(x)={\prod }_{i=0}^{k-1}(x-x_i)$ into $A(x)={\prod }_{\begin{array}{c}i\in Z\\ |Z|=\frac{k}{s}\end{array}}\underset{Z_i}{\underbrace{\prod _{{\omega }^{\prime }\in {\Omega }_i}(x-{\omega }^{\prime })}}.$ We notice that $Z_0(x)=x^{|{\Omega }_0|}-1$ (as its roots are the elements of a group of order dividing $|{\Omega }_0|$) entails $Z_i(x)=x^{|{\Omega }_0|}-{\omega }^{i|{\Omega }_0|}$ (multiplying all terms by a constant ${\omega }^i$ in an integral domain), which is as sparse as it can be. More formally: every element of ${\Omega }_i$ is of the form ${\omega }^i{\omega }^{sj}$ for $j\in [[0,l-1]]$. Thus $Z_i({\omega }^i{\omega }^{sj})=({\omega }^i{\omega }^{sj}{)}^{|{\Omega }_0|}-{\omega }^{i|{\Omega }_0|}=({\omega }^i{)}^{|{\Omega }_0|}({\omega }^{sj}{)}^{|{\Omega }_0|}-{\omega }^{i|{\Omega }_0|}={\omega }^{i|{\Omega }_0|}1-{\omega }^{i|{\Omega }_0|}=0.$ So every element of ${\Omega }_i$ is a root of $Z_i(x)$. Moreover, $Z_i(x)$ is a degree $|{\Omega }_0|=l$ polynomial so has at most $l$ roots: $Z_i(x)$’s only roots are ${\Omega }_i$.

With this little observation, we’ve reduced the number of leaves of the recursion tree for the divide-and-conquer multiplication of the factors of $A(x)$ from $k$ to $k/s$.