精确线性代数

基于向量内积算法的Winograd加速算法
Strassen算法

基于中国剩余定理的模算法
$p$ -adic算法
数值算法求精确解
Wiedemann算法与线性代数中的黑箱方法

精确线性代数主要包含两个互相关联的方面:一是一般PID和域 (特别地如有限域 $F_p$ )上的线性代数问题,二是作为特例,整数环 $\mathbb{Z}$ 和有理数域 $\mathbb{Q}$ 上的线性代数问题.精确线性代数的很多问题同数值线性代数是相对应的,如线性方程组求解,矩阵的特征值与特征向量以及矩阵的各种标准型等.但由于精确计算与数值计算的本质不同,二者在算法设计思想等方面则有根本的不同.简单举例来说,如果对 $\mathbb{Z}$ 上的矩阵进行直接运算,为了保持计算的精确性,将不得不采用有理数的代数运算或引入Bezout等式进行消元.但两种方法随着计算规模的增长，都会出现矩阵元素的大小快速增加的现象，这将使得多数直接算法 (如通常的Gauss消元法)复杂度很高.尽管可以证明,在考虑了矩阵元素规模增长造成的计算复杂度增长之后,多数矩阵精确计算(包括求解线性方程组,计算行列式等)复杂度仍是多项式型的[1],但由于多项式指数相当大,对于实用计算造成很大的困难,难以直接应用.因此,在涉及整数与有理系数矩阵的计算中,一种典型的思路是借助整数的模运算,将主要的计算化归到有限域( $\mathbb{Z}/p\mathbb{Z}$ ,往往也不限于 $p$ 为素数)上进行,再将结果恢复为整数或有理数的形式;另一种典型的思路是通过某种近似,将主要的计算视为数值型计算,通过数值运算得到主要结果后再恢复为精确结果.对于后者,由于数值型算法先于精确算法得到了比较充分的发展,故在精确线性代数中常常要借用数值线性代数的思想与方法,我们在本章中将多次遇到.

本章集中介绍精确线性代数中有普遍应用的几个问题的算法:

快速矩阵乘法,这是很多矩阵算法的基础,在数值线性代数中也有广泛的应用.
线性方程组求解,包括系数矩阵非奇异的情形与一般系数矩阵下 (多相应于不定方程组)的Diophantine解.
行列式的计算,这在计算矩阵的各种标准形时有广泛的应用.
矩阵特征多项式与极小多项式的计算.
矩阵的标准形式,包括Hermite标准形,Smith标准形与Jordan标准形等.

由于精确线性代数的算法很多,且仍在迅速发展中,我们只能介绍一些基本的算法.更多的算法可以参考我们给出的文献以及它们的文献目录.

快速矩阵乘法

在有关矩阵的计算中,矩阵乘法具有基础性的意义.对于 $n$ 阶矩阵的乘法,常规算法具有 $O(n^3)$ 的复杂度.自1968年Strassen[2]发现一种基于分治策略的快速矩阵乘法算法以来,矩阵乘法复杂度的阶数已由3降到2.376[1].下面我们回顾两个经典算法,它们在实际中有着重要的应用.而更多算法虽然渐进复杂度更低,但由于算法过于复杂,且对于有限规模的问题所需运算更多,因而并不实用,可参考[3].[4]

基于向量内积算法的Winograd加速算法

以下讨论主要来自文献[5].

算法1(Winograd内积算法) 设 $x=(x_1,\cdots,x_n)^T$ , $y=(y_1,\cdots,y_n)^T$ ,记 $\xi=\sum\limits_{j=1}^{\lfloor n/2\rfloor}x_{2j-1}x_{2j}$ , $\eta=\sum\limits^{\lfloor n/2\rfloor}_{j=1}y_{2j-1}y_{2j}$ ,则内积 $(x,y)$ 可由下式给出:

$\begin{equation*} (x,y)= \begin{cases} \sum\limits^{\lfloor n/2\rfloor}_{j=1}(x_{2j-1}+y_{2j})(x_{2j}+y_{2j-1})-\xi-\eta,& \text{$n$是偶,} \\ \sum\limits^{\lfloor n/2\rfloor}_{j=1}(x_{2j-1}+y_{2j})(x_{2j}+y_{2j-1})-\xi-\eta+x_ny_n,& \text{$n$是奇.} \end{cases} \end{equation*}$

将这种算法用于 $C=AB$ 的矩阵元素运算时,由于减少重复计算 $\xi$ , $\eta$ ,可使计算所需的乘法次数减半,但同时使所需的加法运算增加.Winograd算法也是 $O(N^3)$ 的算法,仅适用于小规模的矩阵求积运算,且由于该算法破坏了向量内积的整体间运算,同时增加了内存开销,因而其算法改进价值并不很大.

Strassen算法

Strassen算法(1968)是一种分治策略的算法.它以分块矩阵运算为基础.

下面介绍改进型Strassen算法,它较原始算法[2]需要更少的矩阵加法运算[1].

算法2(Strassen算法) 设 $A$ , $B$ 为 $n$ 阶矩阵,必要时通过补充零行零列加以扩充,将其分块:

$\begin{gather*} A= \begin{bmatrix} A_{11}& A_{12}\\ A_{21}& A_{22} \end{bmatrix} ,B= \begin{bmatrix} B_{11}& B_{12}\\ B_{21}& B_{22} \end{bmatrix} ,C=AB= \begin{bmatrix} C_{11}& C_{12}\\ C_{21}& C_{22} \end{bmatrix}. \end{gather*}$ 进行如下递归运算:

若 $n\leq l$ ( $l$ 为递归下界),做直接乘法.
计算 $\begin{gather*} S_1=A_{21}+A_{22},S_2=S_1-A_{11},S_3=A_{11}-A_{21},S_4=A_{12}-S_2,\\ T_1=B_{12}-B_{11},T_2=B_{22}-T_1,T_3=B_{22}-B_{12},T_4=T_2-B_{21}. \end{gather*}$
计算 $\begin{gather*} P_1=A_{11}B_{11},P_2=A_{12}B_{21},P_3=S_4B_{22},P_4=A_{22}T_4,\\ P_5=S_1T_1,P_6=S_2T_2,P_7=S_3T_3. \end{gather*}$
计算 $\begin{gather*} U_1=P_1+P_2,U_2=P_1+P_6,U_3=U_2+P_7,U_4=U_2+P_5,\\ U_5=U_4+P_3,U_6=U_3-P_4,U_7=U_3+P_5. \end{gather*}$
返回 $\begin{equation*} \begin{bmatrix} U_1& U_5\\ U_6& U_7 \end{bmatrix}. \end{equation*}$

以上算法的正确性直接代入即可验证.可以看出,每次递归需要7次乘法与15次加法,从而其算法复杂度是 $O(N^{\log_27})\simeq O(N^{2.808})$ .

下面考虑一个技术细节,即对于阶数不是 $2^k$ 的矩阵添加零行零列的问题.很容易想到两种方案,一是在必要时才考虑添加,即在递归过程中,遇到矩阵阶数为奇数的情形则给它添加一个零行(或零列);二是统一添加,即在计算的开始首先考察矩阵的阶数,若它不满足 $2^k$ 的要求,即给它添加若干个零行零列,使之满足,而在正式计算过程中则不需再考虑矩阵阶数的问题.直观的理论分析可以知道,由于第一种方式一方面将添加零行零列的工作分成许多次完成,增加了很多冗余的判断,另一方面,它逐次添加零行零列的结果是使得矩阵在内存中的存储位置非常零碎,存储结构十分混乱,这两方面因素造成其计算效率大大下降.经实际测试,对于 $1000\times 100$ 阶与 $100\times 1000$ 阶的系数在 $0$ 到 $10$ 之间的整系数矩阵乘法,第二种方案的效率平均较第一种提高了30倍左右.

我们还对随机生成的浮点数矩阵进行了测试,并与经典算法给出的结果做对比,并未见数值稳定性有明显下降.

Strassen算法在之后有许多推广,最优渐进复杂度可以降到 $O(N^{2.376})$ .但在实际中,仅当 $N$ 极大时才有价值,故通常并不采用.可参考[4][3].

线性方程组

线性方程组的重要性是显而易见的,它不但是许多实际问题的数学模型,而且构成许多算法的基础.在"数值线性代数"部分中我们已经介绍了工程计算中使用的数值算法,在这里我们将介绍线性方程组的精确求解,主要讨论整数系数与多项式系数线性方程组的有理解,一些方法可以推广到一般整区上.

基于中国剩余定理的模算法

求解一般线性方程组的一般方法在任何一本线性代数教材中都可以找到(例如,参见[6]).其典型做法是对线性方程组的系数矩阵与增广矩阵进行行初等变换,将其化为行相抵的行既约阶梯形阵(row-reduced echelon,RRE),即如下形式(最后的0可能是子方阵,也可能没有): $\begin{equation*} \begin{bmatrix} 0&\cdots&0&*&\cdots&0&\cdots&0&\cdots\\ & & & & &*&\cdots&0&\cdots\\ & & & & & &\cdots&\cdots&\cdots\\ & & & & & & &*&\cdots\\ & & & & & & & &0 \end{bmatrix} \end{equation*}$ 即具有如下特点:

非0行的最左非0元素非0,且此1所在行的其余元素均为0;
各非0行最左非0元素的位置,随行号增加而右移,若有零行均排在最后.

详言之,若非零矩阵 $B$ 满足:存在一列整数 $1\le k_1<k_2<\cdots<k_r\le n$ ,其中 $1\le r\le m$ ( $r$ 即矩阵 $B$ 的秩),使得 $B(i,k_i)\neq 0$ , $B(i,j)=0$ , $1\le j<k_i$ ,且若 $r<m$ ,则第 $r+1,\ldots,m$ 行均为0.称 $K=\{k_1,k_2,\cdots,k_r\}$ 为 $B$ 的既约阶梯(reduced echelon,RE)序列, $B(i,k_i),1\le i\le r$ 称为RRE矩阵的对角元素.若 $A$ 与 $B$ 行相抵,则也称 $B$ 为 $A$ 的RRE, $K$ 为 $A$ 的行RE序列.对于已经化为这种形式的系数矩阵与增广矩阵,我们很容易判断线性方程组是否有解并求出其一般解.在下面我们考察的情形中，非0行往往以一个公共元素 $d$ 开始,对于整系数线性方程组的求解,由此可以避免中间计算过程出现分数.

在本节中,我们将介绍一种适用于整系数与多元多项式系数线性方程组求解的算法[7].在如下算法中,我们需要判断采用模同态的可用性,其判断标准由如下定义的矩阵的正则RE序列表征.为了叙述方便,我们首先引入[6]中的子阵的记号:设 $1\le i_1<i_2\cdots<i_p\le m$ , $1\le j_1<j_2<\cdots<j_q\le n.$ $m\times n$ 阶矩阵 $A=(a_{ij})$ 中位于第 $i_1,\cdots,i_p$ 行和第 $j_1,\cdots,j_p$ 列交叉处的元素按原序排成的方阵称为 $A$ 的一个 $p\times q$ 阶子阵,记为 $A\begin{pmatrix}i_1&\cdots&i_p\\j_1&\cdots&j_q\end{pmatrix}.$ 记 $M_j$ 为矩阵 $A$ 的前 $j$ 列构成的子矩阵.定义序列 $J=\{j_1,\cdots,j_r\}$ ,其中 $j_h$ 为最小的整数 $j$ 使得 $\mathrm{rank}(M_j)=h$ .由行变换不改变列向量之间的相关性质可知, $J$ 即上文定义的行RE序列,而且是唯一的.对于非零矩阵 $A$ ,可以找到一列互异整数 $h_1,\cdots,h_r$ , $1\le h_t\le m,1\le t\le r$ 满足 $A\begin{pmatrix}h_1&\cdots&h_s\\j_1&\cdots&j_s\end{pmatrix}\neq 0,1\le s\le r.$ 若记 $h_{r+1},\cdots,h_m$ 为其余的整数,则 $H=(h_1,\cdots,h_m)$ 构成 $1,\cdots,m$ 的一个排列.记所有这样的序列 $H$ 构成集合 $\mathcal{P}_A$ ,则 $\mathcal{P}_A$ 中,存在一个(按照字典序)最小的序列 $I_A=\{i_1,\cdots,i_m\}$ .这相当于说,对于 $1\le t\le r$ , $i_t$ 为依次选出来的最小的整数使得这些行向量线性无关,而对于 $r+1\le t\le m$ , $h_t$ 被按照原序排好.若 $A$ 为零矩阵,我们很自然地定义 $I_A=\{1,\cdots,m\}$ .在下面的算法中,我们将用 $J_A\times I_A$ 的字典序来判断实施的模同态是否改变了某些"本质的"内容,即我们认为 $(J_A,I_A)>(J_B,I_B)$ ,当且仅当 $J_A>J_B$ ,或 $J_A=J_B,I_A>I_B$ .

在定义了矩阵的正则RE序列之后,我们引入如下的正则RRE矩阵.对于 $m\times n$ 阶矩阵 $A$ ,定义 $\begin{equation*} \bar{A}_H(k,j)= \begin{cases} A \begin{pmatrix} h_1&\hdotsfor{5}&h_r\\ j_1&\cdots&j_{k-1}&j&j_{k+1}&\cdots&j_r \end{pmatrix},&1\le k\le r,\\ 0,&\text{其他}. \end{cases} \end{equation*}$ 对于利用 $J_A,H_A$ 的定义可知, $\bar{A}_H$ 为RRE矩阵,且可以证明, $\bar{A}_H$ 与 $A$ 行相抵.特别地,对于 $H_A=I_A$ ,将其记为 $\bar{A}$ ,并称为 $A$ 的正则RRE.很显然,这一标准形式是唯一确定的.对于 $H=(h_1,\cdots,h_m)$ ,定义 $\delta_H(A)=A\begin{pmatrix}h_1&\cdots&h_r\\j_1&\cdots&j_r\end{pmatrix},$ 若 $A$ 为零矩阵,则定义 $\delta(A)=0$ .我们看到 $\bar{A}_H$ 的对角线元素都等于 $\delta(A)_H$ .特别地,对于 $H=I_A$ ,将 $\delta_H(A)$ 简记为 $\delta(A)$ 或 $d$ .

该算法的整体思路,就是计算线性方程组增广矩阵的CRRE,并利用其得到线性方程组的一般解.其中,前者是算法最核心的部分.我们首先讨论后者的算法,即已知增广矩阵的CRRE,求得线性方程组的一般解的问题.为了讨论最一般的情形,我们设要求解的方程组为 $AX=B$ ,其中 $A$ 为 $m\times n$ 阶矩阵, $B$ 为 $m\times q$ 阶矩阵,未知元排列成 $n\times q$ 的矩阵 $X$ ,并设 $A$ 与 $B$ 为 $s$ 元多项式系数的矩阵对于整系数矩阵只要将 $s$ 取为0.

根据线性代数的知识,线性方程组的的一般解集合由一个特解与系数矩阵的零空间表征，即该方程组的任意解都可表达为该特解与零空间中向量的和.因此,我们要求线性方程组的一般解,就是要同时求出其特解和零空间的一组基.

首先考察一个RRE矩阵 $E$ 的零空间.设 $E$ 为 $m\times n$ 阶非零RRE矩阵,其RE序列为 $J_E=\{j_1,\cdots,j_r\},r<n$ ,公共的对角元素为 $d$ .记1到 $n$ 中除RE序列之外的元素为 $1\le k_1<\cdots<k_{n-r}$ .如下构造 $n\times n-r$ 矩阵 $Z$ , $\begin{equation*}\tag{1} Z(j_i,j)=E(i,k_j),Z(k_j,u)=\begin{cases}-d,&u=j,\\0,&u\neq j.\end{cases} \end{equation*}$ 容易证明 $EZ=0$ 且 $Z$ 为列独立阵.当 $E$ 作为矩阵 $A$ 的RRE形时,也即存在 $m\times m$ 阶可逆阵 $U$ 使得 $A=UE$ ,则 $AZ=(UE)Z=0$ ,从而 $Z$ 给出了 $A$ 的零空间的一组基.利用这一结论,可以对线性方程组的增广矩阵 $C=(A,B)$ 做出如下结论:

定理1 设 $C=(A,B)$ 为线性方程组 $AX=B$ 的增广矩阵,其中 $A$ 为 $m\times n$ 阶矩阵, $B$ 为 $m\times q$ 阶矩阵.并设 $C$ 的秩为 $r$ , $E$ 为 $C$ 的RRE形,公共的对角元为 $d$ .当线性方程组有解时,设 $Z'$ 为根据 $E$ 与 $J_C$ 按如上方式构造出的矩阵,设 $Z$ 为 $Z'$ 的前 $n$ 行与前 $n-r$ 列的子阵, $Y$ 为 $Z'$ 的前 $n$ 行与后 $q$ 列的子阵,则 $(d,Y,Z)$ 是 $AX=B$ 的一般解,即 $AY=dB$ 且 $Z$ 张成 $A$ 的零空间.

证明依定理所述,将 $Z'$ 划分为 $\begin{bmatrix}Z&Y\\U&V\end{bmatrix},$ 则应用前述结论可得 $CZ'= (AZ+BU,AY+BV)=(AZ,AY-dB)=0.$ $Z$ 的列无关性容易验证. □

然而,对于整系数线性方程组,如果直接对整系数矩阵进行行变换,并保持中间表达式为整系数,则很容易出现类似"高精度运算"部分中提到的中间表达式膨胀(intermediate expression swell),即经过多次乘法后,中间表达式的位数远远超过计算机的硬件处理上限而不得不采取高精度运算,而这将严重影响程序执行的效率.多项式系数的矩阵也有类似情形.为了避免整数与多项式系数矩阵在计算中出现的中间表达式膨胀问题,自20世纪60年代起已经发展了系统的模算法,其基本思想在引言中已经有所介绍,即通过如下两种同态映射(1) $\mathbb{Z}$ 到 $\mathbb{Z}_p$ 的模同态与(2) $\mathbb{Z}_p[x_1,\ldots,x_s]$ 到 $\mathbb{Z}_p[x_1,\ldots,x_{s-1}]$ 的计值运算,这也等价于 $\mathbb{Z}_p[x_1,\ldots,x_s]$ 模 $x_s-a_s$ 的同态映射将整数环与多项式环上的问题化到有限环 $\mathbb{Z}_p$ 上,这样就限制了主要计算的规模.由于矩阵的初等变换归根到底只涉及矩阵元素的乘法与加减法¹,因此这样的模映射保持如下的图交换: $\begin{equation*} \tag{2} \begin{CD} M_n(\mathbb{Z}[x_1,\ldots,x_s]) @>\mod p>> M_n(\mathbb{Z}_p[x_1,\ldots,x_s])@>\text{计值}>>M_n(\mathbb{Z}_p) \\ @VV{\text{行变换}}V @VV{\text{行变换}}V @VV{\text{行变换}}V\\ M_n(\mathbb{Z}[x_1,\ldots,x_s]) @>\mod p>> M_n(\mathbb{Z}_p[x_1,\ldots,x_s])@>\text{计值}>>M_n(\mathbb{Z}_p) \end{CD} \end{equation*}$ 为了从以上计算的结果得到有理数解或分式解,可采用模同态的中国剩余定理与多项式插值定理(本质上也相当于中国剩余定理)来"恢复"应有结果.这种思路可以判断线性方程组是否有解并求得其一般解.整体来讲,我们的算法包含如下的三个层次,相应于(2)中的三个映射:

(算法的最外层)利用模映射,将 $\mathbb{Z}$ 上的多元多项式系数的矩阵映射到 $\mathbb{Z}_p$ 上的多元多项式系数的矩阵,引用中层算法得到其标准雁阵形,再通过中国剩余定理得到整系数或有理系数多元多项式矩阵的RRE.据此得到线性方程组的通解;
(算法的中层,应用于多元多项式系数的线性方程组)对于 $\mathbb{Z}_p$ 上的多元多项式系数的矩阵,通过计值映射将其化为 $\mathbb{Z}_p$ 上的矩阵,引用内层算法将其化为RRE,再通过多元多项式的插值算法得到 $\mathbb{Z}_p$ 上的多元多项式系数的RRE矩阵;
(算法的内层)通过 $\mathbb{Z}_p$ 上的Gauss消去法,将 $\mathbb{Z}_p$ 上的矩阵化为RRE.

对这个整体思路可以打一个比方,为了减小"信息处理"的工作量,我们先对原问题中需要处理的 $\mathbb{Z}$ 上的多元多项式矩阵先做一个"压缩"映射,对这个"压缩包"进行应有的处理之后,再通过"解压缩"将其恢复为我们需要的处理之后的结果.这时,我们将面临一个重要问题,即怎样保证这样的压缩处理是"无损压缩",从而可以从压缩包中重新恢复我们需要的结果呢?也就是说,我们使用怎样的模映射,才能保证可以通过中国剩余定理与插值算法得到原矩阵的RRE呢?详细说来,可以提出如下的三个问题:

我们所谓"无损压缩",究竟保持什么"信息"稳定,从而可以据此恢复应有结果?
什么样的模映射构成我们所需的"无损压缩"?
我们需要多少这样的模映射才能得到真正的"无损压缩"?实际存在的这种模映射够不够用?

下面我们回答这些问题,并给出线性方程组的求解详细算法.

首先讨论模映射的交换性.引入如下记号:以 $\theta:A\longmapsto A^*$ 表示模映射,以 $\Gamma:A\longmapsto \bar{A}$ 表示计算REE形的行变换.若 $\theta$ 与 $\Gamma$ 可交换,则称模映射 $\theta$ 具有交换性.在这种情况下,我们对 $A^*$ 进行行变换得到 $\bar{A^*}=(\bar{A})^*$ ,从而可以对 $\bar{A^*}$ 应用中国剩余定理恢复原矩阵 $A$ 的REE形.下述定理指出了模映射的交换性的充分条件:

定理2 设 $A$ 为 $m\times n$ 阶矩阵,在模映射 $\theta$ 下的像为 $A^*=\theta(A)$ .若 $(J_A,I_A)=(J_{A^*},I_{A^*})$ ,则 $\bar{A^*}=(\bar{A})^*$ ,即 $\theta$ 可交换.

证明注意到在模映射下,求子矩阵的行列式与模映射可交换.利用 $J_A,I_A$ 的定义可知定理成立. □

我们称满足 $(J_A,I_A)=(J_{A^*},I_{A^*})$ 的模映射称为"可行的(accepted)",否则称为不可行的(rejected).注意到可行性并非模映射交换的必要条件,例如下述矩阵 $A=\begin{bmatrix}p&1&0\\0&p&1\\1&0&p\end{bmatrix}$ 在模 $p$ 映射下不满足上述条件,但模 $p$ 映射对于 $A$ 交换.尽管如此,我们仍然得到了一个有效的判别一个模映射是否可交换的充分条件.

下面的论述表明,不可行的模映射是有限的.在模映射下,矩阵子阵的行列式只可能从非0映为0.根据矩阵秩以及 $J_A,I_A$ 的定义可知,必有 $r(A^*)<r(A)$ 或 $(J_{A^*},I_{A^*})\ge (J_A,I_A)$ .模映射 $\theta$ 可行当且仅当对于 $1\le k\le r(A)$ , $A^*\begin{pmatrix}i_1&\cdots&i_k\\j_1&\cdots&j_k\end{pmatrix}\neq0.$ 换句话说,模映射 $\theta$ 不可行当且仅当对某个 $k$ , $A^*\begin{pmatrix}i_1&\cdots&i_k\\j_1&\cdots&j_k\end{pmatrix}=0,$ 这样的模映射在整数环与多项式环上都是有限的,分别由整数的素因子个数与多项式的次数给出限制.因此,我们可以期望得到充分多的可行的模映射完成计算.当然,这里的"充分"由中国剩余定理与多项式插值公式来决定.

然而,以上给出的条件依赖于事先知道 $A$ 的 $(J_A,I_A)$ ,而这是不可能的.我们在计算中已知的只能是模映射后的矩阵的RE序列等.因此,我们还需要将以上条件换成其他等价条件,用模映射后的矩阵表示出来.确切地说,有以下两个定理,它们分别处理了对 $J_A$ 和 $I_A$ 的检测.

定理3 设 $C$ 为 $m\times n$ 阶非零矩阵,秩 $r<n$ . $E$ 为 $m\times n$ 阶RRE矩阵,其秩 $r'< r$ 或 $r'=r,J_E\ge J_C$ .设 $Z$ 为由 $E$ 根据(1)构造出来的矩阵.若 $CZ=0$ ,则 $J_C=J_E$ .

证明若 $r'<r$ 或 $r'=r,J_E>J_C$ ,则很容易通过代入得到不可能有 $CZ=0$ 的结果. □

根据这一定理,我们在计算中,只要将每步得到的 $Z$ 代入 $CZ=0$ 进行验证即可.通过如下的构造,可以减少该步骤验证的计算.设RRE形矩阵 $E$ 的RE序列为 $J_E=(j_1,\cdots,j_r)$ ,其余行记为 $k_1,\cdots,k_{n-r}$ ,公共对角元为 $d$ ,定义 $E$ 的非对角部分 $\tilde{E}$ 如下: $\begin{equation*}\tag{3} \tilde{E}=E \begin{pmatrix} 1&\cdots&r\\ k_1&\cdots&k_{n-r} \end{pmatrix}, \end{equation*}$ >我们知道, $E$ 中其余部分除对角元 $d$ 外均为0元素.又定义 $C$ 的两个子阵如下: $\begin{equation*} C'= \begin{pmatrix} 1&\cdots&m\\ j_1&\cdots&j_r \end{pmatrix},\quad C$ 若 $Z$ 定义如上,则 $$CZ=C'\tilde{E}-dC$ ,从而与代入 $CZ=0$ 验证等价的条件是 $\begin{equation*}\tag{4} C\tilde{E}=dC$ 实际中我们将采用此式进行代入验证.

定理4 设 $C$ 为 $\mathbb{Z}_p[x_1,\cdots,x_s],s>1$ 上的矩阵, $b+1$ 为根据插值公式的次数要求给出的所需计值点数的上界(将在下面给出).设 $\Psi_{a_i},0\le i\le b$ 为不同的计值映射,且 $C^{(i)}=\Psi_{a_i}(C)$ 的RE序列 $J_{C^{(i)}}=J_C$ ,而 $I_{C^{(i)}}=H$ .则 $H=I_C$ .进一步,若 $E$ 是根据 $\bar{C^{(0)}},\cdots,\bar{C^{(b)}}$ 运用插值方法构造出的RRE矩阵,则 $E=\bar{C}$ .

证明反证法.记 $H=(h_1,\cdots,h_m),I_C=(i_1,\cdots,i_m)$ .若 $H\neq I_C$ ,不妨设最先有 $h_k\neq i_k$ ,也即 $h_k>i_k$ ,显然应有 $1\le k\le r$ .则存在 $b+1$ 个计值映射 $\Psi_{a_i}$ 使得 $\Psi_{a_i}(\delta_k(C))=\Psi_{a_i}\left(\mathrm{det}C\begin{pmatrix}i_1&\cdots&i_k\\j_1&\cdots&j_k\end{pmatrix}\right)=0.$ 但由于不可行映射数目的限制(也即上述方程解个数的限制),由 $\deg(\delta_k(C))\le b$ 知道这是不可能的.从而 $H=I_C$ .根据映射的交换性即得由插值方法构造的 $E=\bar{C}$ . □

上述定理给出了算法过程中的另一个检验步骤.

设 $C=(A,B)$ 为 $m\times n+q$ 阶增广矩阵,且 $\tilde{E}$ 与 $d$ 满足(4).则方程组 $AX=B$ 的一般解 $(d,Y,Z)$ 可以如下构造: $\begin{gather*} \tag{5} Z(j_i,j)=\tilde{E}(i,j),\quad Z(k_j,u)=\begin{cases}-d,&u=j,\\0,&u\neq j;\end{cases}\\ \tag{6} Y(j_i,j)=\tilde{E}(i,n-r+j),\quad Y(k_j,u)=0. \end{gather*}$

下面我们给出插值方法所需的模同态的个数,也即给出矩阵中某些子行列式或多项式次数的上界.以下命题的证明都不困难,只要注意到细节问题即可,故略去.读者可以参考[7].

定理5 设 $B$ 为 $I[x]$ 上的 $m\times m$ 阶矩阵,则对于任意整数 $i:1\le i\le m,$ $\begin{equation*} \deg(\det(B))\le\max\limits_{1\le k\le m}\deg(B(i,k))+\sum\limits_{u=1}^me_u-\min\limits_{1\le u\le m}e_u, \end{equation*}$ 其中, $e_u=\max\limits_{1\le t\le m,t\neq i}\deg (B(t,u)),1\le u\le m$ .

定理6 设 $C$ 为 $m\times n$ 阶非零矩阵, $J_C=(j_1,\cdots,j_r)$ .定义 $f=\max\limits_{j\neq j_t}\max\limits_{1\le i\le m}\deg (C(i,j))$ , $e_u=\max\limits_{1\le i\le m}\deg (C(i,j_u)),1\le u\le r$ ,以及 $e\sum\limits^r_{u=1}e_u,e_0=min e_u$ ,令 $g=f-e_0$ .定义 $b=\begin{cases}e,&g\le 0,\\e+g,&g>0.\end{cases}$ 则对任意 $H\in\mathcal{P}_C$ , $b$ 构成 $\bar{C}_H$ 中元素次数的上界.

以上两个定理提供了插值算法所需的对于多项式矩阵的元素次数的估计.一般来讲,中国剩余定理也要求相应的元素上界估计².然而,在本算法中由于有了(4)的代入检验,就不再需要上界估计了.

这样,我们可将算法整理叙述如下.

算法3(最外层算法PLES) 输入: $\mathbb{Z}[x_1,\cdots,x_s],s\le 0$ 上的 $m\times n'$ 阶矩阵 $C$ ,整数 $n:1\le n<n'$ .其中, $C=(A,B)$ 为线性方程组 $AX=B$ 的增广矩阵, $A$ 为 $m\times n$ 阶非零矩阵, $B$ 为 $m\times q$ 阶矩阵, $q=n'-n$ .³

输出:三元组 $(d,Y,Z)$ .若线性方程组有解,则 $(d,Y)$ 构成 $AX=B$ 的特解, $Z$ 为 $A$ 的零空间的基.若线性方程组无解,则 $d,Y,Z$ 全部设为0.

(初始化)置 $r=0$ .
(模 $p$ 映射)若前述计算已经穷尽素数库 $\mathcal{L}$ ,则算法失败⁴.否则,取出下一个素数 $p\in \mathcal{L}$ ,计算 $C^*\equiv C\pmod{p}$ .若 $C^*$ 为零矩阵,转到第2步.
(应用计值插值算法)运用CPRRE算法,计算 $d^*=\delta(C^*),J^*=J_{C^*},I^*=I_{C^*}$ 以及 $\bar{C^*}$ 的非对角部分 $W^*$ .
(进行模算法的可行性检测)置 $r^*=\mathrm{rank}(C^*)$ .若 $J^*(r^*)>n$ ,置 $d=0,Y=0,Z=0$ ,返回 $(d,Y,Z)$ ⁵.若 $r^*>r$ ,转至第5步.若 $r^*<r$ ,转至第2步.若 $(J^*,I^*)<(J,I)$ ,转至第5步;若 $(J^*,I^*>(j,i)$ ,转至第2步.其余情形,转至第6步.
(中国剩余定理算法初始化)置 $r=r^*,d=0,J=J^*,I=I^*,h=1$ ,置 $W$ 为 $r\times n'-r$ 阶零矩阵.
由中国剩余定理的算法迭代步骤,由 $d,d^*,h,p$ 计算 $d'$ .采用类似的步骤,由 $W,W^*,h,p$ 计算 $r\times n'-r$ 阶矩阵 $W'$ 的元素.置 $h$ 为 $p\cdot h$ .
(相等检验)⁶若 $d'=d$ , $W'=W$ ,转至第8步.否则,置 $d=d'$ , $W=W'$ ,转至第2步.
(代入检验)根据(3),由 $C$ 与 $J$ 构造 $C'$ 与 $$C$ .计算 $C'W$ 与 $$d\cdot C$ ⁷.若 $$C'W\neq d\cdot C$ ,转至第2步.
(构造一般解)若 $r<n$ ,根据(5),由 $J,d,W$ 构造 $r\times n-r$ 阶矩阵 $Z$ ;若 $r=n$ ,置 $Z$ 为零矩阵.根据(6)构造 $n\times q$ 阶矩阵 $Y$ .返回 $(d,Y,Z)$ .

这就给出了我们的算法的主要步骤.

下面的算法给出 $\mathbb{Z}_p[x_1,\cdots,x_s],s\ge 0$ 上矩阵RRE形的计算.注意到,其中大部分步骤与PLES算法是类似的.

算法4(中层算法CPRRE) 输入: $\mathbb{Z}_p[x_1,\cdots,x_s],s\ge 0$ 上的 $n\times n$ 阶非零矩阵 $C$ .

输出: $d=\delta(C),J=J_C,I=I_C$ ,以及 $\bar{C}$ 的非对角部分 $W$ .

(应用Gauss消去法)若 $C$ 为 $\mathbb{Z}_p$ 上的矩阵,由下面介绍的CRRE算法计算 $d=\delta(C),J=J_C,I=I_C$ 与 $W$ .
(算法初始化)置 $r=0,a=p,k=0.$
(计值同态)置 $a=a-1$ .若 $a<0$ ,则算法失败返回⁸.否则,置 $C^*=\Psi_{a}(C)$ .若 $C^*=0$ ,转至第3步.
递归调用CPRRE算法,计算关于 $C^*$ 的相应结果: $d^*=\delta(C^*),J^*=J_{C^*},I^*=I_{C^*}$ 以及 $\bar{C^*}$ 的非对角部分 $W^*$ .
(计值同态的可行性检测)置 $r^*=\mathrm{rank}(C^*)$ .若 $r^*>r$ ,转至第6步;若 $r^*<r$ ,转至第3步.若 $(J^*,I^*)<(J,I)$ ,转至第6步;若 $(J^*,I^*)>(J,I)$ ,转至第3步.其余情况,转至第7步.
(插值算法初始化)置 $r=r^*,d=0,J=J^*,I=I^*$ .若 $r<n$ ,构造 $r\times n-r$ 阶零矩阵 $W$ ,若 $r=n$ ,置 $W=0$ .置 $E(x)=1.$
通过插值算法的一个迭代步骤,由 $d,d^*,E(x),a$ 构造 $d'$ .若 $r<n$ ,通过类似的计算用 $W,W^*,E(x),a$ 构造 $W'$ .若 $r=n$ ,置 $W'=0$ .更新 $E(x)\leftarrow E(x)\cdot (x-a).$ 若 $k=1$ ,转至第11步.
(相等检验)若 $d'=d$ , $W'=W$ ,转至第9步.否则,置 $d=d',W=W'$ ,转至第3步.
(代入检验)若 $r=n$ ,转至第10步.根据(3),由 $C$ 和 $J$ 构造 $C'$ 和 $$C$ .代入 $C'W$ 与 $$d\cdot C$ 中进行检验,若 $$C'W\neq d\cdot C$ ,转至第3步，否则置 $k=1$ .
(次数上界)根据6计算DRRE形矩阵元素次数的上界 $b$ .
(次数检验)若 $\deg(E(x))\le b$ ,转到第3步.否则,置 $d=d',W=W'$ ,返回.

在算法的最内层,我们应用 $\mathbb{Z}_p$ 上的Gauss消去法给出矩阵的REE形.由于这是一个标准的算法,我们只要给出几个关键步骤的详细描述即可.

算法5(内层算法,CRRE) 输入: $\mathbb{Z}_p$ 上的 $m\times n$ 阶矩阵 $C$ .

输出:表征 $C$ 的RRE形的 $\delta(C),J_C,I_C$ ,以及非对角元素 $W$ .

(初始化)置 $J=\varnothing,I=(1,\cdots,m),d=1.$ 记 $C^{(0)}=C$ .
(向前消去步骤)设在向前消去的前步,我们已经得到了,并各有个元素添加到与中.在向前消去的第步,执行如下步骤:
- (寻找主元)在第 $k,k+1,\cdots,m$ 行中,按自左向右逐列扫描的顺序,找到第一个非零元素 $C^{(k-1)}(t,s)$ .将 $s$ 添加到 $J$ 中,即令 $j_k=J_C(k)=s$ ;置 $d\leftarrow d\cdot C^{(k-1)}(t,s)$ ;交换第 $k$ 行至第 $m$ 行的排列顺序得到 $C^{(k)}$ 的一个"预备"形式:将第 $t$ 行排到第 $k$ 行,其余各行依原序排列,也即进行如下置换: $\tau_k(h)=\begin{cases}h,&1\le h<k\text{或}t<h\le m,\\k,&h=t,\\h+1,&k\le h<t.\end{cases}$ 随后将 $I$ 中的元素进行类似的重排:将 $I(\tau_k(h))$ 换成 $I(h)$ .
- 令 $e=(C^{(k)}(t,s))^{-1}$ , $C^{(k)}(k,j)\leftarrow C^{(k)}(k,j)\cdot e,s\le j\le n.$
- (向前消去)对于 $h>k$ 行中第 $s$ 列元素非零者, $C^{(k)}(h,j)\leftarrow C^{(k)}(h,j)-C^{(k)}(h,s)\cdot C^{(k)}(k,j),s\le j\le n$ .
(向后消去)经过向前消去后,我们已经得到了 $C$ 的阶梯形式.设 $\mathrm{rank}(C)=r$ ,则须执行 $r-1$ 个向后消去步骤以得到 $C$ 的RRE形式.在向后消去的第 $k:1\le k\le r-1$ 步,利用第 $r-k+1$ 行将前 $r-k$ 行的第 $j_{r-k+1}$ 列元素消去.可以证明,由此得到的RRE矩阵记为 $\bar{C}$ ,据此可以得到其非对角元素 $W$ .

至此,我们已完成了全部算法的叙述.McClellan在[7]中给出了详细的算法复杂性分析,我们不再给出具体结果.我们仅指出,这种算法的实际执行效率高度依赖于素数库的选择与插值映射的选择.其平均时间效率与在此之前的算法(即非模算法,如借助整数严格除法执行的Exact Division算法)相比,大致提高了 $m$ 倍.

$p$ -adic算法

求解非奇异整系数线性方程组 $Ax=b$ 的模算法的另一种思路是对解进行 $p$ 进制( $p$ -adic, $p$ 为素数)展开,然后"恢复"精确的有理数解.这一算法相当于多项式计算中介绍的Hensel提升方法,是Dixon (1982)[8]提出的.需要注意的是,与上节提到的对称表示不同,本节中 $\mathbb{Z}_p$ 中的元素都位于 $[0,p-1]$ 中.

该算法包括三个主要步骤:

合理选定素数 $p$ ,要求 $p\nmid \mathrm{det}A$ .在 $\mathbb{Z}_p$ 上计算 $A$ 的逆 $C$ ,即 $AC\equiv I\pmod{p}$ , $C$ 的存在性由 $C\equiv A^*/\mathrm{det}A\pmod{p}$ 保证.这一步可以利用 $\mathbb{Z}_p$ 上的Gauss消元法实现.
对于选定的充分大的 $m$ ,计算 $\bar{x}$ ,使得 $A\bar{x}\equiv b\pmod{p}$ . $\bar{x}$ 称为 $x$ 的 $p$ 进展式.
由 $x$ 的 $p$ 进展式得到其有理解.

首先考察如何得到 $x$ 的 $p$ 进展式.执行如下步骤: $\begin{align*} b_0&\leftarrow b,\\ x_i&\leftarrow Cb_i\pmod{p},\\ b_{i+1}&\leftarrow p^{-1}(b_i-Ax_i),i=0,1,2,\ldots. \end{align*}$ 注意到由于 $b_i-Ax_i=b_i-ACb_i\equiv 0\pmod{p}$ ，故上式最后一步中的除法是严格的,所得到的 $b_{i+1}$ 为整数.这样的循环步骤将在计算得到 $x_{m-1}$ 与 $b_m$ 之后结束( $m$ 的定义后面给出).这时,令 $\begin{equation*} \bar{x}=\sum\limits^{m-1}_{i=0}x_ip^i, \end{equation*}$ 我们有 $\begin{equation*} A\bar{x}=\sum\limits^{m-1}_{i=0}p^iAx_i=\sum\limits^{m-1}_{i=0}p^i(b_i-pb_{i+1})=b_0-p^mb_m, \end{equation*}$ 从而得到 $\begin{equation*} A\bar{x}\equiv b\pmod{p^m}. \end{equation*}$

为了从 $\bar{x}$ 得到 $x$ ,我们注意到 $x$ 一般为有理系数向量,可以表达为 $x=f/g$ ,其中 $f$ 为整系数向量, $g|\mathrm{det}A$ .由此可得 $g\bar{x}\equiv f\pmod{p^m}$ .有如下定理:

定理7 设 $s,h>1$ 为整数,存在整数 $f,g$ 满足 $gs\equiv f\pmod{h},\text{且}|f|,|g|\le\lambda h^{\frac{1}{2}},$ 其中 $\lambda=0.618\ldots$ 为方程 $\lambda^2+\lambda-1=0$ 的一个根.设既约分数 $w_i/v_i(i=1,2,\ldots)$ 为 $s/h$ 的连分式展式序列,并记 $u_i=v_is-w_ih$ .若该序列中 $k$ 第一个满足 $|u_k|<h^{\frac{1}{2}}$ ,则 $f/g=u_k/v_k$ .

证明根据连分式展式的性质,我们知道序列 $\{w_i\}$ 与 $\{v_i\}$ 递增而 $\{u_i\}$ 则正负交替而绝对值递降, $w_i/v_i$ 收敛于 $s/h$ .关于连分式展开的基本性质,可以参考[9].

记 $f=gs-th$ ,则 $\left|\frac{s}{h}-\frac{t}{g}=\frac{fg}{hg^2}<\frac{1}{2g^2},\right|$ 从而对任意 $t',g'<g$ ,由 $\begin{equation*} \begin{split} \left|\frac{s}{h}-\frac{t'}{g'}\right|&=\left|\frac{s}{h}-\frac{t}{g}+\frac{t}{g}-\frac{t'}{g'}\right|\\ &\ge\left|\frac{t}{g}-{t'}{g'}\right|-\left|\frac{s}{h}-\frac{t}{g}\right|\\ &\ge\frac{1}{gg'}-\frac{1}{2g^2}\\ &\ge \frac{1}{2g^2} \end{split} \end{equation*}$ 得到, $t/g$ 为 $s/h$ 的一个最佳逼近.由[9]第二章定理1得到, $t/g$ 等于 $s/h$ 的一个连分式展式,记为 $w_j/v_j$ .由于 $w_j$ 与 $v_j$ 互素, $|u_j|\le|f|\le\lambda h^{\frac{1}{2}}$ ,由 $k$ 的定义知 $j\ge k$ .另一方面,由 $u_j=v_js-w_jh$ 及 $u_k=v_ks-w_kh$ 得到 $u_jv_k-u_kv_j\equiv 0\pmod{h}$ .由于 $j\ge k$ , $|ujv_k-u_kv_j|\le(|u_j|+|u_k|)|v_j|<\lambda(\lambda+1)h=h$ ,从而 $u_jv_k=u_kv_j$ ,即 $j=k$ .从而 $f/g=u_k/v_k$ .证毕. □

在上述定理中取 $h=p^m$ ,即可给出 $f/g$ 的一种算法.其中 $m$ 由定理中要求的不等式定义: $\delta\le \lambda p^{m/2}$ ,其中 $\delta$ 为 $g$ 与 $f$ 的元素绝对值的上界,可由Hadamard不等式(参考[6]278页48题)给出:对任意实系数 $n$ 阶可逆阵 $B=(b_{ij})$ 有 $|\mathrm{det}B|\le \prod\limits^{n}_{i=1}\left(\sum\limits^{n}_{k=1}b_{ki}^2\right)^{1/2}=\prod\limits_{i=1}^nl_i,$ 其中 $l_i$ 为列向量的2-范数.由Cramer法则知道, $g$ 与 $f$ 中的元素均为增广矩阵的 $n$ 阶子式,从而可以选取 $n+1$ 个列向量中最长的 $n$ 个向量的长度求得Hadamard界 $\delta$ ,并计算得到 $m=\lceil 2\log(\delta\lambda^{-1})/\log p\rceil$ .随后通过连分式展开的步骤,执行如下过程:

$u_{-1}\leftarrow h,u_0\leftarrow s,v_{-1}\leftarrow 0,v_0\leftarrow 1$ ;
对 $i=0,1,\ldots$ 执行 $q_i=\lfloor u_i/u_{i-1},u_{i+1}=u_{i-1}-q_iu_i,v_{i+1}=v_{i-1}+q_iv_i,$ 直到 $u_k\le h^{1/2}$ .则 $f/g=(-1)^ku_k/v_k.$

在实用中,对于连分式展式部分可以采用Lehmer加速算法,可以参考[10].经过详细的分析,Dixon指出这种算法的复杂度为 $O(n^3(\log n)^2)$ ,接近数值算法的渐近复杂度,而优于采用中国剩余定理的算法复杂度.

进入21世纪以来,利用Hensel提升求解整系数线性方程组的思想由Storjohnan等人进行了系统的发展[11][12][13][14].这种发展包括两方面:一是对原来的算法进行了改进,应用J.G.Dumas等人提出的有限域上运算的方法[15]以及改进的扩展Euclid算法[16]大大地提高了算法的效率;另一方面是将该算法与计算矩阵行列式结合起来[17],在一系列问题中得到应用,如求解不定方程组的Diophantine解[18]等.这些算法大多已经被应用于精确线性代数库LinBox中([19]及其索引文献).(这些算法发展待整理)

数值算法求精确解

如引言中所述,将精确线性代数计算过程压缩的另一种思路是将其化归为浮点数的运算.这是基于两方面的考虑:一是现代计算机往往具有超强的浮点数运算能力,通过将高精度的整数通过一定的计算手段化归为机器精度的浮点数运算,可以大大提高计算的效率;二是在精确线性代数快速发展之前,数值线性代数已经经历了长期而系统的发展,建立了基本线性代数子程序(Basic Linear Algebra Subroutines,BLAS)系统,著名的数值线性代数库LAPACK就是它的一个实现.相应的,数值线性代数对精确线性代数的启发有三个个层次:一是在进行 $\mathbb{Z}$ 或 $\mathbb{Z}_p$ 上的精确计算时,借助机器精度较长的浮点数来实现,而无需借助高精度整数环境;二是类似数值线性代数,也在精确线性代数领域引入BLAS技术,为更"高级"一些的算法提供基本的算法支持;三是在"高级"算法领域,如引言中所述,通过某种"压缩"步骤,将算法中大部分计算化归为浮点数的计算,在得到相应的数值结果之后,通过合理的手段恢复为整数或有理数.我们下面将要介绍的这一算法[20],就是这一方面的典型例子.相较Wan给出的原始证明,我们在下面给出的证明充分利用了连分式展开的性质,因而更简洁.

这一算法具有如下性质:它对于适当良态的(即数值稳定性好的,或者说条件数较小的)非奇异矩阵,能很快给出其精确解;而对于病态的矩阵,则很快判断并退出执行.它首先借助多次迭代,求得线性方程组的一个近似有理数解,当该有理数解构成方程组严格解的一个最佳逼近时⁹,用上节提到的连分式展式将精确解求出¹⁰.下面分别讨论这两个步骤.

对于良态矩阵解的迭代逼近,数值线性代数部分也已经有一些介绍,见"Gauss消元法"相关部分.然而,这种方法并没有将全部计算化归为单精度(或者说,机器精度)的计算,特别在将历次修正累加时,需要利用高精度的计算.因此,对于我们目前的要求来说还不够.我们的迭代计算过程如下.对于输入非奇异整系数矩阵 $A$ 与整数向量 $b$ ,将解初始化为 $x=n/d$ ,其中 $n=0$ , $d=1$ ,将余量 $r$ 置为 $b$ .随后在每次迭代中,执行如下的步骤:在机器精度运算中找到近似解 $\Delta x$ 满足 $A\Delta x=r$ ,选择合适的量 $\alpha$ ,对方程进行放大,即令 $\Delta d=\alpha$ , $\Delta n=(\approx\Delta \alpha\cdot x_1,\cdots,\approx\Delta\alpha\cdot x_n)$ ,将方程解修正为 $n=\alpha n+\Delta n$ , $d=\Delta d\cdot d$ ,余量修正为 $r=\alpha r-A\cdot \Delta n$ .这样,在每步迭代步骤之后,出现在结果中的数,包括近似解的分子与分母以及(扩大后的)余量,都是不超过 $\|A\|_\infty+\|b\|_\infty$ 的整数.从而,迭代过程中我们不需要高精度运算.当近似解达到足够精度之后(判别标准如下),我们进行如下算法中叙述的恢复过程.

算法6(数值方法求稠密线性方程组有理解) 输入: $A$ 为 $m\times m$ 非奇异整系数矩阵, $b$ 为整系数向量.

输出:当 $A$ 为良态矩阵时,快速输出方程 $Ax=b$ 的有理解 $x$ ;否则,快速退出并输出"数值精度不足"的信息.

使用机器精度的浮点运算得到 $A$ 的 $LUP$ 分解(其他分解也可使用).
置解向量的公分母 $d^{(0)}=1.$
置余向量 $r^{(0)}=0$ .
置计数器 $i=0$ .
计算 $A$ 的Hadamard界(定义见" $p-$ adic算法"一节).
重复执行以下步骤,直到:
- i:=i+1.
- 利用第1步得到的 $LUP$ 分解,用浮点运算计算 $\bar{x}^{(i)}=A^{-1}r^{(i-1)}$ .
- 计利用浮点运算算 $\alpha^{(i)}:=\min \left(2^{30},2^{\lfloor log_2(\frac{\|r^{(i-1)}\|_{\infty}}{\|r^{(i-1)}-A\bar{x}^{(i)}\|_{\infty}})-1\rfloor}\right)$ .¹¹
- 若 $\alpha^{(i)}<2$ ,则退出并提示信息"数值精度不足".¹²
- 严格计算整系数向量 $x^{(i)}:=(\approx \alpha^{(i)}\cdot \bar{x}_1^{(i)},\cdots,\approx\alpha^{(i)}\cdot\bar{x}_m^{(i)})$ ,也即使得 $x^{(i)}$ 为满足 $\|x^{(i)}-\alpha^{(i)}\cdot \bar{x}^{(i)}\|_{\infty}\le 0.5$ 的向量.
- 严格计算整数 $d^{(i)}:=d^{(i-1)}\cdot\Delta \alpha^{(i)}$ .
- 严格计算整系数向量 $r^{(i)}:=\alpha^{(i)}\cdot r^{(i-1)}-Ax^{(i)}$ ,这时余量扩大了 $d^{(i)}$ 倍.
置 $k:=i$ .
计算整系数向量 $n^{(k)}:=\sum\limits_{1\le i\le k}\frac{d^{(k)}}{d^{(i)}}\cdot x^{(i)}$ ,注意到 $\frac{d^{(k)}}{d^{(i)}}=\prod\limits_{i<j\le k}\alpha^{(j)}$ .
此时,精确解 $x$ 构成 $\frac{1}{d^{(k)}}\cdot n^{(k)}$ 的一个最佳逼近,其分母有上界 $B$ .可利用连分式展开的方法得到 $x$ .
返回 $x$ .

首先我们来估计每次迭代过后余量的上界,并由此导出以上算法的正确性.

定理8(算法的正确性) 若以上算法中每步计算得到的 $\alpha^{(i)}$ 都有 $\begin{equation*}\tag{7} \alpha^{(i)}\le \frac{\|r^{(i-1)}\|_{\infty}}{2\|r^{i-1}-A\bar{x}^{(i)}\|_{\infty}}, \end{equation*}$ 则以上算法将正确执行,即正常地由于矩阵病态而退出,或者得到正确的结果.且在第 $i$ 步迭代中 $\|r^{(i)}\|_{\infty}=\|d^{(i)}(b-A\frac{1}{d^{(i)}}\cdot n^{(i)})\|_{\infty}\le 2^{-i}\|b\|_{\infty}+\|A\|_{\infty}.$

证明对于输入 $A$ 与 $b$ ,我们只需要证明若每次计算出的 $\alpha^{(i)}\ge2$ ,则算法能够顺利执行完毕并得到正确结果.此时,由 $d^{(i）}=\prod\limits_{1\le j\le i}\alpha^{(j)}\ge 2^i$ 可知循环步骤只能执行有限次而最终退出.

记 $n^{(i)}=\sum\limits_{1\le j\le i}\frac{d^{(i)}}{d^{(j)}}\cdot x^{(j)}$ 为 $i$ 步迭代之后解的分子.我们要估计绝对误差 $\|e^{(i)}\|_{\infty}=\|\frac{1}{d^{(i)}}\cdot n^{(i)}-A^{-1}b\|_{\infty}$ ,由归纳法容易得到 $\begin{gather*} r^{(i)}=d^{(i)}(b-A\frac{1}{d^{(i)}}\cdot n^{(i)}),\\ e^{(i)}=\|\frac{1}{d^{(i)}}\cdot n^{(i)}\|_{\infty}=\frac{1}{d^{(i)}}\|A^{-1}r^{(i)}\|_{\infty}. \end{gather*}$ 在每步迭代中,根据定理的假设, $\|A\bar{x}^{(i)}-r^{(i-1)}\|_{\infty}\le \frac{1}{2\alpha^{(i)}}\cdot \|r^{(i-1)}\|_{\infty}$ .根据 $x^{(i)}$ 的定义,我们有 $\|x^{(i)}-\alpha^{(i)}\bar{x}^{(i)}\|_{\infty}\le 0.5$ .从而, $\begin{equation*} \begin{split} \|r^{(i)}\|_{\infty}&=\|Ax^{(i)}-\alpha^{(i)}\cdot r^{(i-1)}\|_{\infty}\\ &\le \|\alpha^{(i)}\cdot A\bar{x}^{(i)}-\alpha^{(i)}\cdot r^{(i-1)}\|_{\infty}+\|Ax^{(i)}-\alpha^{(i)}\cdot A\bar{x}^{(i)}\|_{\infty}\\ &\le \frac{1}{2}\|r^{(i-1)}\|_{\infty}+\frac{1}{2}\|A\|_{\infty}, \end{split} \end{equation*}$ 据此,可以得到 $\|r^{(i)}\|_{\infty}\le \frac{1}{2^i}\|b\|_{\infty}+\|A\|_{\infty}.$ 误差 $e^{(i)}=\frac{1}{d^{(i)}}\|A^{-1}r^{(i)}\|_{\infty}\le \frac{1}{d^{(i)}}\|A^{-1}\|_{\infty}\left(\frac{1}{2^i}\|b\|_{\infty}+\|A\|_{\infty}\right),\forall i\le 1.$ 设 $k$ 为循环停止时计数值,则 $2B\det (A)e^{k}<\frac{2}{d^{(k)}}\|B\det(A)\cdot A^{-1}\|_{\infty}(2^{-k}\|b\|_{\infty}+\|A\|_{\infty}.$ 根据Cramer法则(可参见[6]), $\det(A)A^{-1}$ 正是 $A$ 的古典伴随矩阵,它的每个元素都是 $A$ 的 $m-1$ 阶子式,从而, $2B\det(A)e^{(k)}\le \frac{2mB^2(2^{-k}\|b\|_{\infty}+\|A\|_{\infty})}{d^{(k)}}<1.$ 这样就得到 $e^{(k)}<\frac{1}{2B\det(A)}$ ,也即 $\|\frac{n^{(k)}}{d^{(k)}}-A^{-1}b\|_{\infty}<\frac{1}{2B\det(A)}$ .同样根据Cramer法则,我们知道 $\det(A)A^{-1}b$ 为整系数向量.于是根据定理7,知道 $A^{-1}b$ 构成 $\frac{n^{(k)}}{d^{(k)}}$ 的一个最佳逼近,从而可以由 $\frac{n^{(k)}}{d^{(k)}}$ 的连分式展开得到 $A^{-1}b$ .这就完成了我们的证明. □

在利用近似解恢复精确有理解的阶段,可以采用[21]中建议的一种方法,可以概括为:首先利用随机算法得到矩阵 $A$ 的最大的不变因子(定义见下,也可见[6])或它的一个因子,随后可以更快地计算上述恢复过程.

定理9(PID上矩阵的Smith标准型) 设 $A$ 为主理想整环(PID)上的矩阵,秩 $\mathrm{rank}(A)=r$ .存在PID上的可逆方阵 $P$ 和 $Q$ ,使得 $\begin{equation*} B=PAQ= \begin{bmatrix} s_1&&&&&\\ &\ddots &&&&\\ &&s_r&&&\\ &&&0&&\\ &&&&\ddots &\\ &&&&&0 \end{bmatrix}, \end{equation*}$ 其中 $s_i|s_{i+1}$ , $i=1,\cdots,r-1$ . $(s_1,\cdots,s_r)$ 称为 $A$ 的不变因子组, $B$ 称为 $A$ 的相抵标准型,也称为Smith标准型.

证明证明可参见[6]定理7.13及345页16题. □

注1 由于整数环 $\mathbb{Z}$ 和域 $F$ 上的多项式环 $F[X]$ 是PID,故以上定理适用.特别地,为了使Smith标准型唯一,有时规定整系数矩阵的不变因子 $s_i\ge 0,1\le i\le r$ ,而多项式系数矩阵的不变因子 $s_i(X)$ 为首一多项式, $1\le i\le r$ .以下我们将采用这一规定.

定理10 对于 $m\times m$ 阶非奇异非奇异方阵 $A$ ,设其最大不变因子为 $s_n$ ,则 $s_nA^{-1}$ 为整系数矩阵.

证明设 $B=PAQ$ 为 $A$ 的Smith标准型,则 $s_nB^{-1}$ 为整系数矩阵,而 $P$ , $Q$ 为整系数可逆方阵,从而 $s_nA^{-1}=Q(s_nB^{-1})P$ 为整系数方阵. □

根据该定理,若 $s_n$ 已知,则算法6中由近似解恢复精确解成为平凡的,因为这时 $s_nA^{-1}b$ 为整系数向量,从而在计算足够精度之后将近似解截断即可.即使仅能得到 $s_n$ 的一个非平凡因子 $s_n^*$ ,将其乘到近似解 $\frac{n^{(k)}}{d^{(k)}}$ 上,就只要计算 $s_n^*\frac{n^{(k)}}{d^{(k)}}$ 的一个分母不超过 $B/s^*_n$ 的最佳逼近即可.下面给出随机性地计算 $A$ 最大不变因子的算法.

算法7(最大不变因子) 输入: $n$ 阶非奇异整系数矩阵 $A$ .

输出:正整数 $s^*_n|s_n$ .

(初始化) $M:=\max\{\lceil\sqrt{n\log \|A\|_{\infty}}\rceil,4000\}$ , $p$ 为约为 $n\log A$ 的素数, $b^{(1)},b^{(2)},c^{(1)},c^{(2)}$ 为随机生成的 $n$ 维整系数列向量. $h:=\lceil 2\log_p(\|A\|_{\infty}^{2n-1}M)\rceil$ ,从而 $q:=p^h\ge \|A\|_{\infty}^{2n-1}M$ .
对,执行如下计算步骤:
- $x^{(k)}=(x^{(k)}_i)^n_{i=1}:=A^{-1}b^{(k)}\in \mathbb{Z}^n_q.$
- $y^{(k)}:=c^{(k)T}x^{(k)}=\sum\limits_{i=1}^nc^{(k)}_ix^{(k)}_i\in \mathbb{Z}_q.$
- 取 $t_n^{(k)}$ 为 $y^{(k)}$ 的分母, $t_n^{(k)}:=\delta(y^{(k)})$ ,从而 $0\le t^{(k)}_n\le \|A\|_{\infty}^n$ .
输出 $s_n^*:=\mathrm{lcm}(t^{(1)}_n,t^{(2)}_2)$ .

定理11(算法的有效性) 算法7中, $s_n^*|s_n$ ,且能以不小于 $1/2$ 的概率给出 $s_n^*=s_n$ 的结果.

为了证明7算法的有效性,首先证明如下引理,它给出了以上随机化过程带来的结果.为了叙述方便,我们引入如下记号:对整数 $x$ 与素数 $p$ ,定义 $g=\mathrm{ord}_px$ 满足 $p^g|x$ , $p^{g+1}\nmid x$ ,称为 $x$ 关于 $p$ 的阶.用 $P$ 表示概率.

引理1 记 $\delta^{(k)}:=\mathrm{lcm}(\delta(x^{(k)}_1),\cdots,\delta(x^{(k)}_n))$ ,显然 $t^{(k)}_n|\delta^{(k)}$ .则对任意素数 $\bar{p}$ ,有结论如下:

$P(\mathrm{ord}_ps_n\neq \mathrm{ord}_p\delta^{(k)})\le \max\{1/M,1/p\}$ ;
$P(\mathrm{ord}_pt^{(k)}_n\neq \mathrm{ord}_p\delta^{(k)})\le \max\{1/M,1/p\}$ .

证明证明参见[21]. □

证明(算法有效性的证明) $s_n^*|s_n$ 根据计算过程容易根据归纳法得出.只需证明后一论断. $\begin{equation*} \begin{split} P(s_n \neq s_n^*)&\le \sum\limits_{\bar{p}|s_n,\bar{p}~\text{prime}}P(\mathrm{ord}_{\bar{p}}s_n^*<\mathrm{ord}_{\bar{p}}s_n)\\ &= \sum\limits_{\bar{p}|s_n,\bar{p}~\text{prime}}P(\mathrm{ord}_{\bar{p}}t_n^{(1)}<\mathrm{ord}_{\bar{p}}s_n\wedge \mathrm{ord}_{\bar{p}}t_n^{(2)}<\mathrm{ord}_{\bar{p}}s_n)\\ &\le \sum\limits_{\bar{p}|s_n,\bar{p}~\text{prime}}(\max\{1/M,1/\bar{p}\})^2\\ &<\frac12. \end{split} \end{equation*}$ □

在算法7中,我们仍然需要求解线性方程组,这是不是陷入了一种循环呢?注意到,在算法中,我们只需在 $\mathbb{Z}_q$ 上求解两个线性方程组,而 $y^{(k)},k=1,2$ 的恢复只需要少量运算,而不像之前算法中需要 $n$ 次恢复运算.这样,当 $n$ 很大时,我们可以利用这种方法进行加速.

注2 由定义可知, $\alpha^{(i)}$ 大致只与矩阵 $A$ 的条件数有关.因此,我们可以只计算一次 $\alpha$ ,而将之后的计算中取 $\alpha$ 相同.这样带来的一种可能性是这样得到的 $\alpha$ 不满足定理8中 (7)式的要求¹³.然而,我们可以通过检查余量的范数是否达到理论预测的值来发现这种情形.每当这种情形发生时,我们将 $\alpha$ 减半,再次代入运算即可.

注3 算法6的数值求解算法不限于 $LUP$ 分解.任何具有一定数值稳定性的算法都可以用来求近似解,这种算法可以包括稠密线性方程组的多种分解方法以及稀疏线性方程组的迭代算法.在[20]一文中就利用这种思想构造了一个稀疏线性方程组的求解算法.

注4 事实上,如果不要求精确解的话,算法6不需要执行完毕,而仅作为一个可以有效地估计精度的数值型求解算法.

注5 [20]指出,这一算法的渐近复杂度与Dixon $p$ -adic算法相同,均为 $O^{\sim}(m^3)$ 阶的.但其实际执行效率则比 $p$ -adic算法好很多.

Wiedemann算法与线性代数中的黑箱方法

黑箱算法与随机性算法概述

在数值线性代数中,黑箱型(black box)算法已经广为人知.这一算法的基本特点是,在算法过程中,矩阵仅通过对向量的作用表现它的性质,而不考虑其内部结构.当矩阵是稀疏或有结构的,从而容易用较少的计算量得到其线性作用时,对于数值计算,迭代算法容易快速收敛到足够精度,从而有效减少计算量,同时不会因为对矩阵的显式操作而破坏矩阵的稀疏性和结构性质.自1986年Wiedemann[22]将数值线性代数中已经广为人知的黑箱型(Black box)算法¹⁴用于计算稀疏线性方程组的精确解以来,精确线性代数中的黑箱型算法已经得到了很大的发展.这一方面表现为越来越多的基于黑箱思想的算法的出现,另一方面表现为预处理等适用于数值型算法的随机化思想在精确计算中得到一席之地.在本节中,我们将着重介绍求解有限域上(事实上也不限于有限域)稀疏线性方程组¹⁵的Wiedemann算法([22],也可参考[1])以及求解 $F_2$ 上齐次线性方程,即寻找向量之间线性相关系数的分块Lanczos算法[23].

在介绍算法之前,先引入随机性算法的一些术语.我们知道,在进行精确计算时,由于问题的复杂性,有时会引入与传统计算不同的做法,即随机性算法.随机性算法往往通过引入完全不同的算法思想,达到经典意义上的"决定性"算法不能达到的速度.在数论算法与多项式算法中有很多问题需要考虑随机性算法.我们通常所说的"随机性"往往包含两层意思,一是算法总能快速地给出结果但不能保证结果的正确性.这种类型算法的典型例子是素性判定中的Solovay-Strassen算法,Rabin-Miller算法及APRCL方法等,可见"素数判定"一章中"概率性的检测方法"一节.或者程序快速返回,但可能带回一个警告信息而没有得到正确结果.我们在前面介绍的线性方程组精确求解的数值型算法正是这种情况.这种算法的执行过程中,往往包含随机化的步骤,如产生随机数,随机矩阵或对已知向量(矩阵)进行随机的线性变换等.特别地,在线性代数领域,这种随机化往往通过随机性的预处理器(pre-conditioner)[24]来实现.提出这种这种算法必须同时给出出错的概率的上界,从而提供这种算法的可靠性,这往往通过有关多项式零点个数的Schwartz-Zippel定理[25]进行估计.这种类型的算法往往称为Monte Carlo型算法.第二种情形是算法总能给出正确的结果,但不能保证每次都能快速执行,而只能在平均意义下,或在一些好的情形能快速给出正确结果,这称为Las Vegas型算法.实际上,我们遇到的很多算法都或多或少带有这种性质,即能够在一些好多情形特别快速地给出正确结果,例如以上给出的基于中国剩余定理的模算法.很容易看到,两种类型的随机性算法的区别不是绝对的,对于同一种算法甚至可以转换为Las Vegas型的或Monte Carlo型的.例如,如果很容易判断计算结果的正确性,(例如求解线性方程组的情形,求解过程一般需要 $O(n^3)$ 的算法复杂度,而判断解的正确性却只要 $O(n^2)$ 次计算.)这时,我们只要反复执行算法,直到我们得到一个正确结果.这时我们可以说算法被化归为Las Vegas型的.对于一个Las Vegas算法,如果我们发现程序执行了过长的时间,例如远远长于平均执行时间时,我们就令算法中止退出,同时返回一个警告.我们可以说,这时算法就是一个Monte Carlo型的算法¹⁶.我们下面介绍的几个算法都带有随机性质.

在本小节的最后,我们介绍上面提到的Schwartz-Zippel定理[25].

定理12(Schwartz-Zippel) 设 $Q\in F[x_1,\cdots,x_n]$ , $Q\neq 0$ 是 $n$ 元多项式.设 $Q_1$ 是 $Q$ 的最简形式, $d_1$ 为 $x_1$ 在 $Q_1$ 中的次数, $Q_2$ 为 $x_1^{d_1}$ 在 $Q_1$ 中的系数多项式.归纳地说,令 $d_i$ 为 $Q_i$ 中 $x_i$ 的次数, $Q_{i+1}$ 为 $Q_i$ 中 $x_i^{d_i}$ 的系数, $1\le i\le n$ .对 $1\le i\le n$ ,令 $x_i$ 在 $I_i\subseteq F$ 中取值,则在 $I_1\times \cdots \times I_n$ 中 $Q$ 至多有 $|I_1\times \cdots \times I_n|\left(\frac{d_1}{|I_1|}+\cdots +\frac{d_n}{|I_n|}\right)$ 个零点.

这一定理给出了多元多项式在其定义集合上零点个数的上界估计.在我们采取的随机化步骤中,往往需要在某个集合中随机取值,构造向量或矩阵,前者如线性递推列特征多项式的Berlekamp-Massey算法,而后者则是随机预处理器的核心步骤.在此过程中,我们常常希望避免某种奇异情形出现,而这种奇异情形往往表现为随机化算法中不小心取到了某个多元多项式的零点.根据Schwartz-Zippel定理,我们可以将随机取值的集合取得充分大,就可使算法由此失败的概率减低.

证明归纳法. $n=1$ 的情形由代数基本定理可得.

设以上命题对 $Q_2\in F[x_2,\cdots,x_n]$ 成立,即其零点集 $\mathrm{Null}(Q_2)$ 有 $|\mathrm{Null}(Q_2)|\le |I_2\times\cdots \times I_n|(\frac{d_2}{|I_2|}+\cdots+\frac{d_n}{|I_n|})$ .若 $(z_2,\cdots ,z_n)\in \mathrm{Null}(Q_2)$ ,则 $\forall x_1\in I_1$ , $(x_1,z_2,\cdots,z_n)\in I_1\times \cdots\times I_n$ 都可能是 $Q_1\in F[x_1,\cdots,x_n]$ 的零点.若 $(z_2,\cdots,z_n)\notin \mathrm{Null}(Q_2)$ ,则至多有 $d_1$ 个 $x_1\in I_1$ 使得 $(x_1,z_2,\cdots,z_n)\in \mathrm{Null}(Q_1)$ .从而 $Q_1$ 在 $I_1\times\cdots I_n$ 中至多有 $\begin{equation*} \begin{split} &|I_1||I_2\times\cdots\times I_n|\left(\frac{d_2}{|I_2|}+\cdots+\frac{d_n}{|I_n|}\right)+d_1|I_2\times \cdots\times I_n|\\ &|I_1\times I_n|\left(\frac{d_1}{|I_1|}+\cdots+\frac{d_n}{|I_n|}\right). \end{split} \end{equation*}$ 个零点. □

Schwartz-Zippel定理有如下推论,其证明很直接,但其结论对于概率性算法的分析有重要意义.

推论1 设 $Q\in F[x_1,\cdots,x_n]$ , $Q\neq 0$ ,且 $x_1,\cdots,x_n$ 均在 $U\subseteq F$ 中随机选取. $Q$ 结果为零的概率不超过 $\frac{D}{\sharp U}$ ,其中 $D$ 为 $Q$ 的次数.

推论2 若 $F=\mathbb{R}$ 或 $F=\mathbb{C}$ , $Q\in F[x_1,\cdots,x_n]$ , $Q\neq 0$ .若 $(x_1,\cdots,x_n)$ 在非零测集 $U\subseteq F$ 上选取,则 $Q(x_1,\cdots,x_n)$ 在概率1的意义下非0.

推论3 若 $F=F_q$ , $Q\in F[x_1,\cdots,x_n]$ , $Q\neq 0$ .若 $(x_1,\cdots,x_n)$ 在 $F$ 中随机选取,则 $Q(x_1,\cdots,x_n)$ 为零的概率不超过 $\frac{D}{\sharp F}$ , $D$ 为 $Q$ 的次数.

为了显示以上定理的作用,我们在随机预处理器中选取一个简单例子来看一看.

所谓预处理步骤,就是我们希望通过一个确定的或随机性的变换 $\mathcal{P}$ ,将给定矩阵 $A$ 映射为 $A'$ ,且 $A'$ 满足某种好的性质.根据要求结果满足的性质不同,预处理器分为许多种[24],例如使主子式非奇异,或保证前 $r,r=\mathrm{rank}A$ 列线性无关等.

我们知道,在经典的Gauss消元法中,对于一般矩阵求解必须附加选取主元的过程.该过程的目的就是为了保证主子式的可逆性,从而能够顺利进行以后的消元步骤.理论上,我们还可不限于对矩阵每个元素的操作,考虑如下的块消元过程: $\begin{equation*} \begin{bmatrix} A&B\\ C&D \end{bmatrix} \begin{bmatrix} I&-A^{-1}B\\ O&I \end{bmatrix}= \begin{bmatrix} A&O\\ C&D-CA^{-1}B \end{bmatrix}. \end{equation*}$ 这样的算法至少有如下的几点好处:首先它借助分块矩阵计算将消元化为了一个递归的过程,这有利于算法复杂度的降低,再者经过化归后,该算法变得容易并行化.然而,这样的算法应用于一般矩阵时将遭遇重大的困难:算法步骤中要求 $A$ 必须可逆,而这并非普遍得到满足的.

如果我们能够通过如下的线性变换 $UMV\rightarrow M'$ 对 $M$ 进行预处理,其中 $U$ 和 $V$ 为预先选取的矩阵,使得 $UMV$ 满足块消元过程中要求的性质,则该算法可以进行,从而充分利用这一算法的优点.实际中,我们常常根据要求,在一定范围内进行随机选取预处理器,使之满足我们所需要的性质.

例1(随机四分变换) 若 $M=\begin{bmatrix}A&B\\C&D\end{bmatrix}$ 为 $(2n)\times(2n)$ 阶非奇异复系数矩阵, $R$ , $S$ 为 $n\times n$ 阶对角阵,其元素在 $U\subseteq \mathbb{C}$ 中随机选取.则 $\begin{equation*} \tilde{M}= \begin{bmatrix} \tilde{A}&\tilde{B}\\ \tilde{C}&\tilde{D} \end{bmatrix}= \begin{bmatrix} I&R\\ I&-R \end{bmatrix} \begin{bmatrix} A&B\\ C&D \end{bmatrix} \begin{bmatrix} I&I\\ S&-S \end{bmatrix} \end{equation*}$ 的每一分块均以概率1非奇异. ◇

证明注意到以上得到的 $\tilde{M}$ 的每一分块的行列式均为 $(r_1,\cdots,r_n)$ 和 $(s_1,\cdots,s_n)$ 的非零多项式. □

以上的随机四分变换是一类更普遍的随机蝶形变换(random butterfly transformation)的特例.这类变换的思想源起于一类复杂网络问题,它提供了一个很方便的分块矩阵预处理器,可以参考[26].从以上例子中我们可以归纳出预处理器的构造思路:

预处理过程即对矩阵 $A$ 进行如下的线性变换 $A\underrightarrow{\mathcal{P}}PAQ$ , $P\in F[x_1,\cdots,x_s]$ , $Q\in F[y_1,\cdots,y_t]$ .
其中 $x_1,\cdots,x_s,y_1,\cdots,y_t$ 在 $F$ (或 $U\subseteq F$ )上随机选取.
由于要证明的性质对应于 $F[x_1,\cdots,y_r]$ 上非零多项式,因此变换后的矩阵在一定概率下满足我们要求的性质.

在构造 $P$ , $Q$ 的过程中,有一些原则是要遵守的:

$P$ , $Q$ 的选取要尽量简单.
$P$ , $Q$ 首先具有很好的性质,如求逆操作易于进行等.
若我们需要进行数值计算,则 $P$ , $Q$ 应具有好的条件数,能够保证数值计算的精度.

线性递推列

在本节,我们首先介绍线性递推序列(linear recurrent sequences)及其极小多项式的概念以及相关算法.

定义1(线性递推列) 设 $F$ 为域, $V$ 为 $F$ 上的线性空间,记 $V$ 上的无穷序列 $(a_i)_{i\in \mathbb{N}},a_i\in V,i\in\mathbb{N}$ 全体构成的线性空间为 $V^{\mathbb{N}}$ .称 $a=(a_i)\in V^{\mathbb{N}}$ 为线性递推列,若存在 $n\in\mathbb{N}$ 及 $f_0,\cdots,f_n\in F,f_n\neq 0$ ,使得 $\sum\limits_{0\le j\le n}f_ja_{i+j}=f_na_{i+n}+\cdots+f_1a_{i+1}+f_0a_i=0$ 对于任意的 $i\in \mathbb{N}$ 都成立. $n$ 阶多项式 $f\equiv \sum_{0\le j\le n}f_jx^j\in F[X]$ 称为 $a$ 的特征多项式,也称为零化多项式或生成多项式.

根据以上定义,将 $F$ 和 $V$ 均取为 $\mathbb{Q}$ ,则我们熟悉的Fibonacci序列是线性递推列, $x^2-x-1$ 构成它的一个特征多项式.根据我们的经验,对于一个线性递推列,在已知其特征多项式及最初的几个取值时,我们可以通过递推的方法快速计算其后记序列.下面我们举一些线性代数中涉及的线性递推列的例子,它们在我们下面的算法中将扮演重要角色.

例2 本例给出线性代数中一些重要的线性递推列.这里最关键的是利用了线性代数中的Caley-Hamilton定理[6].以下所取的基域均为一般域 $F$ .

取 $V=F^{n\times n}$ , $A\in V$ 为任意方阵.则 $a=(A^i),i\in \mathbb{N}$ 为线性递推列, $A$ 的最小多项式与特征多项式都是 $a$ 的特征多项式.
设 $V=F^n$ , $A\in F^{n\times n}$ , $b\in F^n$ 为任意向量.则注意到 $(A^i)$ 的任意特征多项式也将 $(A^ib)$ 化为零,故 $(A^jb)$ 也为线性递推列.由 $a$ 的元素生成的 $V$ 中的线性子空间称为 $A$ 与 $b$ 的Krylov子空间,这一定义与我们在数值线性代数中"迭代法"一章中引入的定义是一致的.
设 $V=F^n$ , $A\in F^{n\times n}$ , $b,u\in F^n$ 为 $V$ 中的任意向量.则 $(A^ib)$ 的任意特征多项式同样将 $(u^TA^ib)$ 化为零,故 $(u^TA^i)b$ 也为线性递推列. ◇

我们看到,在引入线性递推列的特征多项式概念后,二者的结构产生了一定的关系.我们还可通过引入多项式对序列的作用,将无穷序列与多项式更紧密地联系起来,同时引出一些重要的概念.

定义2(多项式对序列的作用) 设 $f=\sum_{0\le j\le n}f_jx^j\in F[X]$ 以及序列 $a=(a_i)_{i\in \mathbb{N}}\in V^{\mathbb{N}}$ ,定义 $f$ 对 $a$ 的作用如下 $f\bullet a=\left(\sum\limits_{0\le j\le n}f_ja_{i+j}\right)_{i\in \mathbb{N}}\in V^{\mathbb{N}}.$ 注意到,在此定义下,常系数对序列的作用即数乘作用,而未定元 $x$ 对 $a$ 的作用相当于平移算子.在以上定义下, $V^{\mathbb{N}}$ 成为 $F[X]$ 上的模¹⁷.

固定 $a\in V^{\mathbb{N}}$ ,容易验证 $a$ 的特征多项式的集合与0一起,构成了 $F[X]$ 中的理想,我们称之为 $a$ 的零化子理想,记为 $\mathrm{Ann}(a)$ .由于 $F[X]$ 为主理想整环,存在 $\mathrm{Ann}(a)$ 的生成元 $m_a$ .按照定义, $m_a$ 的倍数全体构成了 $a$ 的特征多项式集合, $m_a$ 是其中次数最低的首一多项式,我们称之为 $a$ 的最小多项式. $m_a$ 的次数称为 $a$ 的递推阶数.结合2中的几个例子,我们可以对线性代数中的特征多项式,最小多项式等概念有更深刻的理解.

本节的主要目的是给出线性递推序列的最小多项式的构造性算法.以下定理给出了极小多项式的一个判别法则.

定义3(多项式的逆转) 设 $f(x)=f_dx^d+\cdots+f_0\in F[X]$ 为 $d$ 阶多项式,记 $\mathrm{rev}(f)\equiv x^df(x^{-1})=f_0x^d+\cdots+f_d\in F[X]$ ,称为 $f$ 的倒逆(reverse).

定理13 设 $a=(a_i)_{i\in \mathbb{N}}\in F^{\mathbb{N}}$ 为线性递推序列, $$h=\sum_{i\in \mathbb{N}}a_ix^i\in F<a href=$ X$" class="latex-inline" style="vertical-align: -7px">, $f\in F[X]\backslash 0$ , $\deg f=d$ , $r=\mathrm{rev}f$ .

以下命题等价:
- $f$ 为 $a$ 的特征多项式;
- $r$ , $h$ 为低于 $d$ 阶的多项式;
- $h=g/r$ ,其中 $g\in F[X]$ ,且 $\deg g<d$ .
若 $f$ 为 $a$ 的极小多项式,则 $d=\max\{1+\deg g,\deg r\},$ 且 $\mathrm{gcd}(g,r)=1$ .

证明

直接代入按照形式幂级数的乘法规则验证即可.
$d\ge \deg r$ ,故 $d\ge \max\{1+\deg g,\deg r\}$ .若 $d>\deg r$ ,则 $x|f$ ,且 $f/x$ 为 $a$ 的特征多项式,与 $f$ 极小矛盾.设 $u=\mathrm{gcd}(g,r)$ ,则 $f^*\equiv f/\mathrm{rev}(u)$ 为 $d-\deg r$ 次多项式, $r/u=\mathrm{rev}(f^*)$ ,且 $(r/u)\cdot h=g/u$ 为 $\deg g-\deg u$ 次多项式.于是 $r/u$ 为 $a$ 的特征多项式, $u\in F$ . □

注6 从 $h$ 的定义可以看出, $h=g/r$ 为序列 $a$ 的生成函数(generating function).

在域 $F$ 上,取线性空间 $V$ 即基域 $F$ ,若给定一个线性递推列 $a$ ,且我们能够预先知道其最小多项式 $m_a$ 的次数上界.例如在前面的例2的3中, $m_a$ 的最高次数不会超过矩阵 $A$ 的阶数.在这种情况下, $m_a$ 可以通过求解如下的Padé逼近问题得到解答: $\begin{equation*}\tag{8} h\equiv \frac{s}{t}\pmod x^{2n},x\nmid t,\deg s<n, \deg t\le n,\mathrm{gcd}(s,t)=1. \end{equation*}$ 我们看到, $(s,t)=(g,r)$ 正是如上问题的一组解.由Padé逼近问题解的唯一性,即得如下算法.

算法8(Berlekamp-Massey) 输入:域 $F$ 上线性递推列 $a\in F^{\mathbb{N}}$ ,给定 $a$ 的递推阶数的一个上界 $n$ 以及 $a$ 的前 $2n$ 个元素 $a_0,\cdots,a_{2n-1}\in F$ .

输出: $a$ 的最小多项式 $m_a\in F[X]$ .

$h:=a_{2n-1}x^{2n-1}+\cdots+a_1x+a_0$ ,使用扩展Euclid算法计算 $s,t\in F[X]$ 且 $t(0)=1$ ,使得(8)成立.
$d:=\max\{1+\deg s,\deg t\}$ .
返回 $\mathrm{rev}_d(t)$ .

经过分析,这一算法能在 $O(n^2)$ 的域 $F$ 运算步骤内给出 $a$ 的最小多项式.

线性方程组的Wiedemann算法

下面我们考虑将线性方程组的解与线性递推列的最小多项式联系起来.我们首先考虑线性方程组 $Ax=b$ ,系数矩阵 $A\in F^{n\times n}$ 非奇异,则该线性方程组有唯一解.设 $m=m_{b}=\sum_{0\le j\le d}m_jx^j\in F[X]$ 为线性递推序列 $a_b=(A^ib)_{i\in \mathbb{N}}$ 的最小多项式,则 $m\bullet a=m(A)b=\sum_{0\le j\le d}m_jA^jb=0$ ,从而我们得到 $A\cdot (-m_0^{-1})\sum\limits_{1\le j\le d}m_jA^{j-1}b=b,$ 即 $j=-m_0^{-1}\sum_{1\le j\le d}m_jA^{j-1}b$ 为方程组的解.据此分析,我们得到如下算法:

算法9 给定域 $F$ 上的非奇异 $n$ 阶方阵 $A$ 与任意 $n$ 维向量 $b$ ,本算法计算线性方程组 $Ax=b$ 的解 $x=A^{-1}b$ .

计算线性递推列 $(A^ib)_{i\in \mathbb{N}}$ 的最小多项式 $m\in F[X]$ .
返回 $y:=-m_0^{-1}\sum_{1\le j\le d}m_jA^{j-1}b.$

以上给出了Wiedemann算法的整体思路,但其中还有一个重要问题没有解决,这就是算法第一步中要求的 $n$ 维向量组成的线性递推列的最小多项式的计算.但由例2中3,我们知道 $a_b\equiv(A^jb)_{j\in \mathbb{N}}$ 与 $a_{ub}\equiv(u^TA^jb)_{\mathbb{N}},u\in F^{\mathbb{N}}$ 的特征多项式有很重要的联系.详言之,我们知道 $a_b$ 的每个特征多项式都将 $a_{ub}$ 零化,即 $\mathrm{Ann}(a_b)\subseteq \mathrm{Ann}(a_{ub})$ ,从而 $m_{ub}|m_b$ ;另一方面,若 $f\in \mathrm{Ann}(a_{ub}),\forall u\in F^{\mathbb{N}}$ ,则 $f\in \mathrm{Ann}(a_b)$ ,否则总可取 $u=f\bullet a_b$ 使得 $f\bullet a_{ub}\neq 0$ .这样,我们有如下的随机性算法,若该算法顺利执行完毕,则能给出正确的结果.

算法10 给定域 $F$ 上 $n$ 阶方阵 $A$ 与 $n$ 维向量 $b$ ,该算法计算线性递推序列 $a_b\equiv (A^jb)_{j\in \mathbb{N}}$ 的最小多项式 $m_{b}$ .

随机选取 $u\in F^{n}$ ,计算 $(u^TA^ib)_{i\in\mathbb{N}}\rightarrow a$ .
运用算法8计算 $m:=m_{ub}\in F[X]$ .
校验步骤:若 $m(A)b\neq 0$ ,进行步骤1.
返回 $m$ .

注7 如果计算结果在校验步骤中失败,意味着 $m_{ub}|m_b$ 但 $m_{ub}\neq m_b$ .在以上算法中我们采取的是丢弃策略,即认为这并没有提供给我们任何有价值的信息.事实上,我们看到, $m_{ub}$ 至少提供了 $m_b$ 的一个非平凡因子.因此我们也可以尝试如下策略来对以上算法进行加速,将第3步的校验步骤扩展为以下两种方式中的任意一种:

$m:=m\cdot m_{ub}$ ,若 $m(A)b\neq 0$ ,则令 $b:=m(A)b$ ,返回步骤1.
$m:=\mathrm{lcm}(m,m_{ub})$ ,若 $m(A)b\neq 0$ ,则返回步骤1.[22]

这样做法的正确性很容易证明.

在给出一个随机性算法后,我们急需解决的一个问题是算法的有效性.与确定性算法不同的是,仅仅知道一个随机性算法的正确性是远远不够的,¹⁸我们还需要知道,这一算法到底以多大概率命中正确结果,从而能够很快返回.为此,我们还需要对 $F^{\mathbb{N}}$ 作为 $F[X]$ 上模的结构有更清楚的认识.

设 $f\in F[X]$ 为 $d$ 次多项式,记 $\langle f\rangle=F[X]\cdot f$ 为由 $f$ 生成的多项式理想.定义 $M_f:=\{a\in F^{\mathbb{N}}|f\bullet a=0\}$ ,则 $M_f$ 为 $F[X]$ 上的循环模,且 $M_f=F[X]\bullet c$ ,其中 $c=(0,\cdots,0,1,c_d,c_{d+1},\cdots)$ , $c_d,c_{d+1},\cdots$ 为根据前 $d$ 个元素递推得到的序列元素.这时 $c,x\bullet c,\cdots,x^{d-1}\bullet c$ 构成了 $M_f$ 的一组基.我们看到,多项式 $F[X]$ 在 $M_f$ 上的作用给出了 $F[X]$ 到 $M_f$ 的满同态,从而决定了如下的模同构: $M_f\cong F[X]/\mathrm{Ann}(M_f)=F[X]/\langle f\rangle.$ 这样的结构给出了如下的引理,它是我们分析算法有效性的基础.

引理2 设 $A\in F^{n\times n}$ , $b\in F^n\backslash 0$ . $f\in F[X]$ 为 $(A^ib)_{i\in \mathbb{N}}$ 的最小多项式.存在 $F-$ 线性的满射 $\psi:F^n\rightarrow F[X]/\langle f\rangle$ ,使得 $\forall u\in F^n$ ,以下两个命题等价:

$f$ 为 $(u^TA^ib)_{i\in \mathbb{N}}\in F^{\mathbb{N}}$ 的最小多项式;
$\psi(u)$ 为单位,即 $(\psi(u),f)=1$ .

证明由前面的讨论,存在 $F[X]/\langle f\rangle$ 到 $M_f$ 的同构 $\phi:F[X]/\langle f\rangle\rightarrow M_f$ .我们还可定义 $F-$ 线性映射 $\psi^*:F^n\ni u\longmapsto \psi^*(u)=(u^TA^ib)_{i\in \mathbb{N}}\in M_f$ ,并由 $f$ 为 $M_f$ 的极小多项式得到 $\psi^*$ 为满射.因此可以构造 $F-$ 线性满映射 $\psi=\phi^{-1}\circ \psi^*:F^n\rightarrow M_f,$ 用交换图表示如下: $\begin{equation*} \xymatrix{ F[X]/\langle f\rangle\ar[dr]_{\phi}&&F^n\ar[dl]^{\psi^*}\ar[ll]_{\psi=\phi^{-1}\circ \psi^*}\\ &M_f } \end{equation*}.$

这时,我们有 $\begin{equation*} \begin{split} f=m_{\psi^*(u)}&\Longleftrightarrow \forall g\in F[X](g\bullet \psi^*(u)=0\Leftrightarrow f|g)\\ &\Longleftrightarrow \forall g\in F[X]((g \bmod f)\bullet \psi(u)=0\Leftrightarrow g\bmod f=0)\\ &\Longleftrightarrow \forall h\in F[X]/\langle f\rangle (h\bullet\psi(u)=0\leftrightarrow h=0)\\ &\Longleftrightarrow \psi(u)\text{为单位}. \end{split} \end{equation*}$ □

定理14 设 $U\subseteq F$ , $A\in F^{n\times n}$ , $b\in F^n\backslash 0$ , $f$ 为 $(A^ib)_{i\in\mathbb{N}}$ 的极小多项式, $d=\deg f$ ,则 $u$ 从 $U^n$ 中随机选取时, $f$ 为 $(u^TA^ib)_{i\in\mathbb{N}}$ 的最小多项式的概率 $p\ge 1-d/\sharp U$ .

这一定理给出了算法有效性的估计,并可据此给出Wiedemann算法的平均复杂度估计.

证明记 $u=(u_1,\cdots,u_n)^T=u_1e_1+\cdots+u_ne_n\in F^n$ 为 $F^n$ 中的任意向量,由 $\psi$ 线性得到 $\begin{equation*} \begin{split} \psi(u)&=u_1\psi(e_1)+\cdots+u_n\psi(e_n)\\&=(u_1h_1+\cdots+u_nh_n)\bmod f, \end{split} \end{equation*}$ 其中 $h_1,\cdots,h_n\in F[X]$ ,且次数低于 $d$ . $\psi(e_j)=h_j\bmod h_j$ .令 $y_1,\cdots,y_n$ 为不定元, $r=\mathrm{res}_x(y_1h_1+\cdots+y_nh_n,f)\in F[y_1,\cdots,y_n]$ ,则 $r$ 的总次数不超过 $d$ .由于 $\psi(u)\text{单位}\Leftrightarrow r(u_1,\cdots,u_n)\neq 0$ ,但 $\psi$ 满,故 $r\neq 0$ .从而 $u_i\in U$ 随机选择时, $p\ge 1-d/\sharp U$ . □

根据该定理,可以给出Wiedemann算法9的平均复杂度估计,证明可参考[1].

定理15 若 $F$ 中含有不少于 $2n$ 个元素,则算法9的复杂度为 $4nc(A)+O(n^2)$ ,其中 $c(A)$ 为矩阵 $A$ 左乘向量所需的计算量,与矩阵 $A$ 的结构有关.

对于元素很少的域 $F_q$ ,上面的随机性算法不适合于直接应用.这时,我们可以采取两种方式:一是对 $F_q$ 进行适当的域扩张,从中进行元素的随机选择;此时算法复杂度大致要乘上一个 $O(M(\log_qn))$ 的因子,其中 $M(\log_qn)$ 是进行 $F_q$ 上的 $n$ 次多项式乘法的运算量.第二种方式是利用注记7中类似的方法,即使某次计算失败了,仍然从中提取出我们想要的信息.

1. 注意到,以下的算法步骤中并不显式包含整数的除法,这是我们构造该交换图的基础.然而,必须注意,这里出现的三个"行变换"是执行完全相同的操作,在这个意义下,该交换图总是成立的.而后面将会讨论的模同态的交换性是在如下意义下:对于一个已知的矩阵,要通过行变换求得其REE形,这一变换对于模映射前后的两个矩阵可能是不同的,例如矩阵 $\begin{bmatrix}0&p\\p&0\end{bmatrix}$ 在模 $p$ 映射下得到零矩阵,二者的RE序列显然不同,从而化为REE形所需行变换也是不同的.如果模映射前后的两个矩阵需要执行的行变换完全相同,则模映射自然是交换的,这构成"交换"的充分条件.对此,后面会给出更为详细的分析.

2. 这种估计往往由Hadamard不等式提供,见"p进展开方法"一节.

3. 在该算法中,我们还需要一个相当规模的素数库 $\mathcal{L}$ .通常取该库中的素数 $p$ 满足 $p^2\lessapprox \ell$ ,其中 $\ell$ 为计算机硬件计算的字长.这样的素数库可以预先存储在程序库中,也可以动态生成.该库的生成并非本算法讨论的内容.

4. 若素数库采取动态生成的方式,则不应有此问题.

5. 注意到此时,线性方程组右端向量组 $B$ 包含列独立与系数矩阵列空间的向量,从而方程组无解.

6. 这一检验成立是中国剩余定理算法结束的必要条件(不是充分条件).此处引入该检验,是为了让随后的"代入检验"成功的可能性更大,从而减少代入检验的运算量.

7. 此处要运用矩阵乘法.对于一般矩阵的乘法快速算法可以参考"快速矩阵乘法"一节.对于多元多项式系数的矩阵乘法(在CPRRE算法中要用到),也可采用计值插值的方法进行.

8. 注意!这里包含一个与PLES算法很重要的不同之处.由于PLES算法中被模掉的素数可以在全体整数中选取,因此其算法主体总能成功(只要它调用的CPRRE算法成功);然而,由于CPRRE算法的计值点只能在 $\mathbb{Z}_p$ 中选取,而这是有限的,因此可能不成功.鉴于此,通常我们选取模同态时,要求模掉的素数 $p$ 充分地大:只要满足 $p^2\lessapprox \ell$ 的条件,从而 $\mathbb{Z}_p$ 上的运算可以高效地执行.

9. 这种说法看上去很自然,但并不准确.严格来讲,应该是精确解构成近似解的一个最佳逼近.见下所述.

10. 据[20]的说法,这种思想来源于分母有限大的有理数是离散的.

11. 上述选择是基于如下考虑:选择2的幂次使得以下的计算变为简单的移位的过程,从而更为搞笑;若 $A$ 为充分良态的,则余量可能已经非常之小,需要加一个限制以保证以下的乘法不会溢出,这一限制被选择为 $2^{30}$ 同样是为了使计算结果在机器精度之内.

12. 注意到,这样的情形仅发生在 $A$ 为病态的时候.

13. 事实上,即使按照原算法中建议的去做,也可能由于浮点数运算的舍入误差出现这种情况(尽管已经尽量做了避免).

14. 可参考数值线性代数部分"迭代法"一章,也可参考其中所列出的参考文献.

15. 这里的"稀疏"条件只是为了保证矩阵对向量的线性作用能以较少的计算量得到估计.对于某些有结构的矩阵,如Vandermonde阵,由于也满足这一要求,因而同样适合应用黑箱型算法.

16. 它至少还返回了一个警告.从这种意义上说,它还优于通常意义下的Monte Carlo型算法.

17. 直观地讲,模就是"环上的线性空间".关于其严格定义以及一些一般性质,可参考[6].

18. 如果我们这里的正确性仅仅指,"如果按照正常步骤计算并正常退出,则算法返回正确结果"的话.

19. Lanczos算法最先是由Lanczos提出用于数值求解稀疏线性方程组的,参见"数值线性代数"部分.准确地说,Wiedemann算法正是类似Lanczos算法的黑箱技术在精确线性代数领域的对应.

参考文献

[1]Joachim von zur Gathen and Jürgen Gerhard, Modern Computer Algebra, Cambridge University Press, 2002.

[2]Strassen V., Gaussian Elimination is not Optimal, Numer. Math. 13 (1969), 354-356.

[3]Coppersmith D., Winograd S., Matrix multiplication via arithmetic progressions, proceedings of the nineteenth annual acm symposium on theory of computing, 1987.

[4]Pan V., How Can We Speed Up Matrix Multiplication?, SIAM Review 26 (1984), no.3.

[5]Winograd S., A New Algorithm for Inner Product, IEEE Trans. Comp. 17 (1968), 693-694.

[6]张贤科，许甫华, 高等代数学, 清华大学出版社, 北京, 2004.

[7]McClellan M., The exact solution of systems of linear equations with polynomial coefficients, Journal of the Association for Computing Machinery 20 (1973), no.4, 563-588.

[8]Dixon J., Exact solution of linear equations using p-adic expansions, Numer. Math. 40 (1982), 137-141.

[9]Rockett A., Szusz P., Continued fractions, World Scientific, 1992.

[10]Donald E. Knuth, The art of computer programming, volume 2 (3rd ed.): seminumerical algorithms, Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1997.

[11]Chen Z., Storjohann A., An implementation of certified linear system solving for integer matrices, 2004, Poster at ISSAC.

[12]Chen Z., Storjohann A., A fast implementation of system solving for integer matrices, 2004, Poster at ECCAD'2004.

[13]Storjohann A., High-order lifting and integrality certification, Journal of symbolic computation 36 (2003), 613-648.

[14]Storjohann A., The shifted number system for fast linear algebra on integer matrices, Journal of complexity 21 (2005), 609-650.

[15]Dumas J., Giogi P., Pernet C., FFPACK: finite field linear algebra package, ISSAC'04, 2004.

[16]Wang X., Pan V., Acceleration of euclidean algorithm and rational number reconstruction, SIAM J. Comput. 32 (2003), no.2, 548-556.

[17]Pan V., Complexity of parallel matrix computations, Theor. Comp. Sc. 54 (1987), 65-85.

[18]Giesbrecht M., Effient parallel solution of sparse systems of linear diophantine equations, PASCO'97, 1997, 1-10.

[19]Villard G., Some recent progress in exact linear algebra and related questions, ISSAC'07, 2007.

[20]Wan Z., An algorithm to solve linear systems exactly using numerical methods, 2005, preprint submitted to elsevier science.

[21]Wang X., Pan V., Acceleration of euclidean algorithms and extensions, ISSAC'02, 2002.

[22]Wiedemann D., Solving sparse linear equations over finite fields, IEEE transactions on infomation theory 32 (1986), no.1, 54-62.

[23]Montegomery P., A block lanczos algorithm for finding dependencis over gf(2), Advances in cryptology, EUROCRYPT'95, LNCS 921, Springer-Verlag, 1995, 106-121.

[24]Chen L., Eberly W., Kaltofen E. etc., Efficient matrix preconditioners for black box linear algebra, linear algebra and its applications 343-344 (2000), 119-146.

[25]Schwartz J., Fast probabilistic algorithms for verification of polynomial identities, Journal of association for computing machinery 27 (1980), no.4, 701-717.

[26]Parker S., A randomizing butterfly transformation useful in block matrix computations, (1995).