Reed–Solomon error correction算法

Reed-solomon 为前向纠错码（forward error correcting）的子类-线性分组码（line block code）的一种，通过添加冗余数据来抵抗信号传输过程中发生的翻转（bit flip）和丢失（erasure）。从信息论角度看，FEC属于信道编码（channel coding），与信源编码（source coding）的目的正好相反。

通信

通信(Communication)自古以来是人类社会进步不可或缺的过程，一般来说，通信不外乎以下几个阶段。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-aiaLxNP3-1633441699581)(https://km.woa.com/gkm/api/img/cos-file-url?url=https%3A%2F%2Fkm-pro-1258638997.cos.ap-guangzhou.myqcloud.com%2Ffiles%2Fphotos%2Fpictures%2F202110%2F1633344043-7143-615ada2bae69c-609080.png&is_redirect=1)]

比如在原始社会，人类将要传达的消息（input message）通过肢体动作、手势等来进行编码（coding），然后经过可见光反射至眼球表面并投影至视网膜中（channel），视网膜捕获到光信号后，接着通过大脑解码成图像信息，最后大脑对图像信息进行解析。

而对于现代数字通信，为了使数据传输更高效。对上述编码解码过程进行了更细致的划分如下（参考Coding theroy）。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-dFo18suK-1633441699584)(https://km.woa.com/gkm/api/img/cos-file-url?url=https%3A%2F%2Fkm-pro-1258638997.cos.ap-guangzhou.myqcloud.com%2Ffiles%2Fphotos%2Fpictures%2F202110%2F1633344072-0074-615ada4801d21-816124.png&is_redirect=1)]

其中信源编码（source coding）就是我们常用的压缩解压缩技术，用于移除数据中的冗余信息从而减少数据体积。密码编码（cyptographic coding）用于保证数据的安全性。信道编码（channel coding）也称为FEC，用于抵抗传输信道中的发生的数据错误。线路编码（line coding）用于将传输数据通过物理世界中的某些信号来表示。

本文要介绍的 Reed-Solomon Code 便是处于信道编码阶段。

信道编码之线性分组码

一般来说，待编码的原始数据称为消息（message），编码后的数据称为码字（codeword），信道编码主要分为线性分组码和卷积码（convolutional code），其中线性分组码是指对原始数据划分为固定大小的一个分组，然后对一个分组的数据进行编码或解码的操作。

最简单的线性分组码为重复码（repetition code），比如待发送的消息为“Hello”，经过三次重复码可以编码为码字“HHHeeelllooo”，对于$n$次重复码，其有效数据占比为$\frac{1}{n}$。显然重复码编码效率低下。

再考虑异或码（xor code），异或相当于每个bit位的${\mathbb{Z}}_2$加法运算。假设待发消息为“101”和“010”，则异或生成冗余数据“111”。这三个数据丢失任意一个，都可以通过剩下两个异或还原出来。假如一个分组有$M$个消息，通过异或生成一份冗余数据，其有效数据占比为$\frac{M}{M+1}$。其冗余数据开销非常小，但是其缺点也非常明显，仅能恢复出分组里的一个消息。

于是希望有这样一种线性分组码，一个分组里有$M$个消息，可以生成$N$个冗余数据，丢失任意小于等于$N$个消息都可以恢复出来。为了形式化描述该问题，假设一个消息由长度为$k$的字符串组成，消息的每个字符的取值范围属于某个集合$F$，给定一个长度为$k$的消息字符串$m$，希望通过某种算法或函数$f(x)$生成字符串长度为$n$的码字$c$
$$m⟼f(m)$$
，使得码字$c$中的任意$k$个字符都能还原出剩下的$n-k$个字符。

构建Reed-Solomon码

学过数值分析，我们知道在多项式插值中，如果我们观测到$k$个二维平面上的点集 { x i , y i ∣ 0 ≤ i < k } \{x_i,y_i|0\le i< k\} {xi​,yi​∣0≤i<k}，我们可以通过拉格朗日插值（Lagrange interpolation）找到最高次幂不超过$k-1$的多项式$f(x)=a_0+a_{1}x+\cdots a_{k-1}{x}^{k-1}$。写成矩阵的形式如下
$$
\begin{pmatrix} 1 & x_0 & \cdots & x_0^{k-1} \
1 & x_1 & \cdots & x_1^{k-1} \
\vdots & \vdots & \ddots & \vdots \
1 & x_{k-1} & \cdots & x_{k-1}^{k-1}
\end{pmatrix}
\cdot
\begin{pmatrix}
a_0 \
a_1 \
\vdots \
a_{k-1}
\end{pmatrix}

\begin{pmatrix}
y_0 \
y_1 \
\vdots \
y_{k-1}
\end{pmatrix}
KaTeX parse error: Can't use function '$' in math mode at position 9: 记上述方阵为 $̲V_{k\times k}$，…
\det(V_{k\times k})=\prod_{0\le i<j\le k}(x_j-x_i)
$$
因此当$x_i,x_j$两两不同时，此方阵可逆。

于是多项式系数$f(x)$可以得出
$$
\begin{pmatrix}
a_0 \
a_1 \
\vdots \
a_{k-1}
\end{pmatrix}

(V_{k\times k})^{-1}
\begin{pmatrix}
y_0 \
y_1 \
\vdots \
y_{k-1}
\end{pmatrix}
KaTeX parse error: Can't use function '$' in math mode at position 22: …们可以构造一种编码方案，多项式$̲f(x)$的系数看作需要编码的…
\begin{pmatrix}
a_0 \
a_1 \
\vdots \
a_{k-1}
\end{pmatrix}

\begin{pmatrix} 1 & x_{i_1} & \cdots & x_{i_1}^{k-1} \
1 & x_{i_2} & \cdots & x_{i_2}^{k-1} \
\vdots & \vdots & \ddots & \vdots \
1 & x_{i_k} & \cdots & x_{i_k}^{k-1}
\end{pmatrix}^{-1}
\begin{pmatrix}
y_{i_1} \
y_{i_2} \
\vdots \
y_{i_k}
\end{pmatrix}
$$
实际上，这个思路正是Irving Stoy Reed和Gustave Solomon于1960发表的论文（特定有限域上的多项式编码）POLYNOMIAL CODES OVER CERTAIN FINITE FIELDS。

另外证明此编码方式满足线性，满足线性的本质是因为有基（basis）的存在。不妨记向量${\alpha }_i^T=(x_0^{i-1},x_1^{i-1},\cdots ,x_k^{i-1})$为$V_{k\times k}$的第$i$列，于是码字$y_0{y}_1\cdots y_{n-1}$为向量${\alpha }_i$的线性组合
$$a_0{\alpha }_0+a_1{\alpha }_1+\cdots a_{k-1}{\alpha }_{k-1}=\vec{y}$$
由于$\det (V_{k\times k})\ne 0$，所以此$k$个向量线性无关。记消息集合为$M$，码字集合为$C$，那么对于任意两个码字${\vec{y}}_m=m_0{\alpha }_0+m_1{\alpha }_1+\cdots m_{k-1}{\alpha }_{k-1}\in C$和${\vec{y}}_n=n_0{\alpha }_0+n_1{\alpha }_1+\cdots n_{k-1}{\alpha }_{k-1}\in C$，其线性组合$k_1{\vec{y}}_m+k_2{\vec{y}}_n=(m_0+n_0){\alpha }_0+(m_1+n_1){\alpha }_1+\cdots (m_{k-1}+n_{k-1}){\alpha }_{k-1}$仍为向量${\alpha }_i$的线性组合。即$k_1{\vec{y}}_m+k_2{\vec{y}}_n\in C$。

实际上，可逆矩阵的本质就是线性变换。而线性分组码的可以看做是将低维空间的消息通过线性变换映射至高维空间，使得码字之间的最小汉明距离大于1。

几何视角

如果把多项式$f(x)=a_0+a_{1}x+\cdots a_{k-1}{x}^{k-1}$重写成
$$-f(x)+a_0+a_{1}x+\cdots a_{k-1}{x}^{k-1}=0$$
并记向量
$${\vec{a}}^T=(a_0,a_1,\cdots ,a_{k-1})$$
和
$${\vec{x}}^T=(-f(x),1,x,\cdots ,x^{k-1})$$
于是有
$${\vec{a}}^T\cdot \vec{x}=0$$
可以看作是由法向量$\vec{x}$确定的$k$维超平面，且点$\vec{a}$就在此超平面上。类比二维空间中任意两条不平行的直线相交于一点和三维空间中任意三个两两不平行的平面相交于一点（两个平面交于一条线，该线与第三个平面交于一点），可以知道$k$维空间中任意$k$个线性无关的超平面相交于一点。即有$k$个线性无关的由法向量${\overrightarrow{x_i}}^T=(-f(x_i),1,x_i,\cdots ,x_i^{k-1})$确定的$k$维超平面，这些超平面必定相交于一点，此点即为系数$\vec{a}$。

系统编码（Systematic encoding）

实际应用中，通常希望生成的码字的前$k$个字符正好就是原始消息。这样如果收到前$k$个字符就无需解码，满足此特性的编码也叫系统编码。同样可以根据拉格朗日插值来构造编码矩阵。这里假设编码多项式为
$$g(x)=b_0+b_{1}x+\cdots b_{k-1}{x}^{k-1}$$
消息
$$m=a_0{a}_1\cdots a_{k-1}$$
约定好的$n$个求值点为
{ x 0 , x 1 , x 2 , ⋯   , x n − 1 } \{x_0,x_1,x_2,\cdots,x_{n-1}\} {x0​,x1​,x2​,⋯,xn−1​}
于是构造$g(x)$，使之满足
$$\{\begin{array}{l}g(x_0)=a_0\\ g(x_1)=a_1\\ ⋮\\ g(x_{k-1})=a_{k-1}\end{array}$$
等价的矩阵表示为
$$
\begin{pmatrix} 1 & x_0 & \cdots & x_0^{k-1} \
1 & x_1 & \cdots & x_1^{k-1} \
\vdots & \vdots & \ddots & \vdots \
1 & x_{k-1} & \cdots & x_{k-1}^{k-1}
\end{pmatrix}
\begin{pmatrix}
b_0 \
b_1 \
\vdots \
b_{k-1}
\end{pmatrix}

\begin{pmatrix}
a_{0} \
a_{1} \
\vdots \
a_{k-1}
\end{pmatrix}
$$
求出系数${\vec{b}}^T=(b_0,b_1,\cdots ,b_{k-1})$后，再对$g(x)$分别代入$x_{k+1},\cdots ,x_{n-1}$得出冗余字符$g(x_{k+1})\cdots g(x_{n-1})$，最终得到码字$a_0{a}_1\cdots a_{k-1}g(x_{k+1})\cdots g(x_{n-1})$。

有限域与域扩张

提到Reed-Solomon码，不得不提到域（field），域的定义很简单，记元素集合为 S = { s 0 , s 1 , ⋯   } S=\{s_0,s_1,\cdots\} S={s0​,s1​,⋯}，且此集合中存在两种二元运算符 { △ , ◯ } \{\bigtriangleup,\bigcirc\} {△,◯}，满足交换律，结合律，分配律，且两种算符下均存在单位元$e_{△}$和$e_{◯}$，任意元素均存在逆元且唯一，则称$S$为定义在 { △ , ◯ } \{\bigtriangleup,\bigcirc\} {△,◯}上的域，假设算符优先级$△<◯$，如果在算符$◯$上不满足任意元素均在逆元的，则称为环（ring）。如果集合$S$的元素个数有限，则称为有限域（finite field），也称作Galois field。一般记$q$个元素的有限域为$F_q$，常见的有限域为模素数集合 ${\mathbb{Z}}_p$。

为什么是素数？反证法，不妨假设模合数集合${\mathbb{Z}}_n$，对于元素$m\in {\mathbb{Z}}_n$，且记最大公约数$GCD(m,n)=t>1$，显然对于任何元素$s\in {\mathbb{Z}}_n,(s\cdot m)\mathrm{m}\mathrm{o}\mathrm{d}n|t$。即${\mathbb{Z}}_n$中不存在$m$的逆元，${\mathbb{Z}}_n$只能是环。另外当$n=p$为素数时，可通过扩展欧几里得算法（EEA）求出其唯一的逆元。

通常来说，前面推导Reed-Solomon的时候默认的数域为定义在算符 { + , × } \{+,\times\} {+,×}上的实数集合$\mathbb{R}$，由于计算机表示数字的精度有限（理论上可以用字符串存储数字实现接近无限精度的计算，但是计算开销会很大）。一般实际会在有限域上进行操作。

考虑一个$n$维空间，该空间的每个维度的取值范围均为域${\mathbb{Z}}_p$，则此空间的每个向量$\vec{v}\in {\mathbb{Z}}_p^n$，共有$p^n$个元素，如果定义向量在算符 { × } \{\times\} {×}上的操作和${\mathbb{Z}}_p$上一样，那么显然该向量的每个元素在 { × } \{\times\} {×}上均有逆元，于是通过有限域$F_p$构造了一个新的有限域$F_{p^n}$。

如果对上述$n$维空间的每个元素${\vec{v}}^T=(v_0,v_1,\cdots ,v_{n-1}),v_i\in {\mathbb{Z}}_p$映射为多项式$v(x)=v_0+v_{1}x+\cdots +v_{n-1}{x}^{n-1}$，并记此多项式集合为$F_p[x]$（多项式的系数在数域$F_p$上），并假设$p(x)$为定义在$F_p$上的最高次幂为$n-1$的不可约多项式（irreducible polynomial），即不存在解$a\in F_p$使得$p(x)=(x-a)h(x)$。则$F_p[x]/p(x)$构成了域。证明如下

对于任意$v(x)\in F_p[x]$，采用长除法都可以找到余数$r(x)\in F_p[x]$，使得$v(x)=d(x)p(x)+r(x)$，所以映射
$$v(x)⟼r(x)$$

$$F_p[x]⟼F_p[x]/p(x)$$

一一对应，这里余数$r(x)$的集合即为$F_p[x]/p(x)$，显然$F_p[x]/p(x)$是环，其在算符 { + } \{+\} {+}上可逆有单位元，现在需要证明其在 { × } \{\times\} {×}上也存在逆元。由于$p(x)$为不可约多项式，所以$GCD(p(x),r(x))=1$，即可以使用扩展欧几里得算法找到$s(x),t(x)\in F_p[x]/p(x)$，使得
$$GCD(p(x),r(x))=s(x)p(x)+t(x)r(x)=1$$
即$t(x)r(x)\mathrm{m}\mathrm{o}\mathrm{d}p(x)=1$，也即$t(x)$为$r(x)$在$F_p[x]/p(x)$上的逆元，所以$F_p[x]/p(x)$也是域，其元素个数为$p^n$。

实际上引入不可约多项式$p(x)$的本质就是引入了一个新元素$\beta \notin F_p$，且满足$p(\beta )=0$，即$p(x)$在加入$\beta$这个元素之后可以分解出因子$(x-\beta )$。即通过引入一个根$\beta$替换$x$，域$F_p$扩张成了新域$F_p[x]/p(x)$，也叫域扩张（field extension）。新域元素个数为$p^n$。

举一个简单的例子，考虑一个${\mathbb{Z}}_2$上的多项式$p(x)=x^3+x+1$，因为$p(0)=1,p(1)=1$，所以$p(x)$在${\mathbb{Z}}_2$上不可分解，即为不可约多项式， F 2 [ x ] = { b 0 + b 1 x + b 2 x 2 + ⋯ + ∣ b i ∈ Z 2 } F_2[x]=\{b_0+b_1x+b_2x^2+\cdots+|b_i\in\mathbb{Z_2}\} F2​[x]={b0​+b1​x+b2​x2+⋯+∣bi​∈Z2​}为多项式构成的环（$F_2[x]$不是域，并非每个元素在$\times$上存在逆元），考虑模$p(x)$
$$F_2[x]⟼F_2[x]/p(x)$$
并代入根$\beta$

$\beta$	$F_2[x]/p(x)$
0	0
$\beta$	x
${\beta }^2$	$x^2$
${\beta }^3$	$x^3\mathrm{m}\mathrm{o}\mathrm{d}{x}^3+x+1=x+1$
${\beta }^4$	$(x+1)x=x^2+x$
${\beta }^5$	$(x^2+x)x\mathrm{m}\mathrm{o}\mathrm{d}{x}^3+x+1=x^2+x+1$
${\beta }^6$	$(x^2+x+1)x\mathrm{m}\mathrm{o}\mathrm{d}{x}^3+x+1=x^2+1$
${\beta }^7$	$(x^2+1)x\mathrm{m}\mathrm{o}\mathrm{d}{x}^3+x+1=1$

可以看到$F_2[x]/p(x)$上除了算符$\times$上的零元，所有元素均可以通过$\beta$的幂次来表示，这里$\beta$称为生成元（generator）。根据此表，可以快速算出$F_{p^n}$域上任意两个元素的乘积。现在我们知道，给定一个模素整数域$Z_p$，只要找到该域上一个最高次幂为$n-1$的不可约多项式，便可以快速的构造出域$F_{p^n}$。由于计算机内存表示数据是以字节为最小单位，因此在真正工程实现Reed-Solomon码时，一般采用有限域$F_{2^8}$。

Cauchy-RS码

回顾前面的介绍，一个有效的编码矩阵，只需要满足其行列式不为0，且任意子矩阵也保持该性质，除了Vandarmonde矩阵，实际上，Cauchy矩阵也具有该性质，Cauchy矩阵的结构如下
$$\left(\begin{array}{cccc}(x_1-y_1{)}^{-1}& (x_1-y_2{)}^{-1}& \cdots & (x_1-y_n{)}^{-1}\\ (x_2-y_1{)}^{-1}& (x_2-y_2{)}^{-1}& \cdots & (x_2-y_n{)}^{-1}\\ ⋮& ⋮& \ddots & ⋮\\ (x_m-y_1{)}^{-1}& (x_m-y_2{)}^{-1}& \cdots & (x_m-y_n{)}^{-1}\end{array}\right)$$
记为$C_{m\times n}$其中元素$x_i,y_i$为某个域$F_q$中的元素，且$x_i-y_j\ne 0,1\le i\le m,1\le j\le n$。当Cauchy矩阵为方阵时，其行列式为
$$\det (C)=\frac{{\prod }_{i=2}^n{\prod }_{j=1}^{i-1}(x_i-x_j)(y_k-y_i)}{{\prod }_{i=1}^n{\prod }_{j=1}^n(x_i-y_j)}$$
从行列式表达式可以看出只要$x_i,y_i$两两不相同，行列式非0，因此Cauchy方阵可逆。在工程实现中，经常会选择Cauchy矩阵作为编码矩阵，因为其逆矩阵计算复杂度低，只需$\mathcal{O}(n^2)$。其证明过程可以参考论文An XOR-Based Erasure-Resilient Coding Scheme

参考文献

[1] Reed, Irving S., and Gustave Solomon. “Polynomial codes over certain finite fields.” Journal of the society for industrial and applied mathematics 8.2 (1960): 300-304.

[2] Bloemer, Johannes, et al. “An XOR-based erasure-resilient coding scheme.” (1995).