特殊线性方程组

Yule-Walker问题的Durbin算法
一般右端向量的Levinson算法
正定对称Toeplitz阵求逆Trench算法

数值分析的一个基本原则是:求解任何问题都应充分利用它的特殊的结构特性.在数值线性代数中,这意味着当问题中出现诸如对称性,定型性和稀疏性等特性时,需要将适用于求解一般矩阵的算法修改,使其效率更高.这将是本章的主题.我们的主要目的是设计一些特殊的 $LU$ 分解的专用算法.

$LDM^T$ 和 $LDL^T$ 分解

由前述 $LU$ 分解容易证明以下的 $LDM^T$ 分解定理:

定理1如果 $A\in\mathbb{R}^{n\times n}$ 的所有顺序主子阵都是非奇异的,则存在唯一的单位下三角阵 $L$ 和 $M$ 和唯一的对角阵 $D=\mathrm{diag}\{d_1,\cdots,d_n\}$ 满足 $A=LDM^T$ .

特别地,在 $A$ 为对称方阵的情形,成立以下 $LDL^T$ 分解:

定理2 如果 $A=LDM^T$ 是非奇异对称阵 $A\in\mathbb{R}^{n\times n}$ 的 $LDM^T$ 分解,则 $L=M$ .

由以上定理可知,对对称阵施行 $LU$ 或 $LDL^T$ 分解时,工作量可减半.在第 $j$ 步,由于 $M=L$ 且假定 $L$ 的前 $j-1$ 列已知,则 $L(j,1:j-1)$ 也为已知.定义 $v(1:j)$ 为 $DL^Te_j$ 的前 $j$ 个分量,则 $\begin{equation*} v(1:j)= \begin{bmatrix} d(1)L(j,1)\\ \vdots\\ d(j-1)L(j,j-1)\\ d(j) \end{bmatrix}. \end{equation*}$ 于是,向量 $v(1:j-1)$ 可通过对 $L$ 的第 $j$ 行做简单数乘得到.由 $L(1:j,1:j)v=A(1:j,j)$ 的第 $j$ 个方程,有关系式 $v(j)=A(j,j)-L(j,1:j-1)v(1:j-1)$ .

相应地有以下算法:

算法1( $LDL^T$ ) 如果 $A\in\mathbb{R}^{n\times n}$ 对称且存在 $LDL^T$ 分解,则本算法计算单位下三角阵 $L$ 和对角阵 $D$ ,其元素覆盖 $A$ . $\begin{align*} \label{eq:LDLT} &\text{在算法的第}j\text{步,计算}v(1:j).\\ &\text{由}v(j)=A(j,j)-L(j,1:j-1)v(1:j-1)\text{计算}v(1:j).\\ &A(j,j)\leftarrow v(j).\\ &A(j+1:n,j)\leftarrow (A(j+1:n,j)-A(j+1:n,1:j-1)v(1:j-1))/v(j). \end{align*}$

本算法大约减少了一半的工作量,需要 $n^3/2$ 个flop.可以证明,当 $A$ 对称且正定时,以上算法不但能够顺利执行完,而且非常稳定[1].如果 $A$ 对称但非正定,则需要结合对称地选主元,即在算法的第 $j$ 步比较 $A(j:n,j:n)$ 得到其中绝对值最大的元素并通过置换阵 $P$ 的对称作用 $A\leftarrow P^TA^TP$ 将其换到 $A(j,j)$ 位置,使之后的算法稳定执行.

Cholesky分解

定义1(正定矩阵) 若对所有非零向量 $x\in\mathbb{R}^n$ 有 $x^TAx>0$ ,则称矩阵 $a\in\mathbb{R}^{n\times n}$ 是正定的.

可以证明,若 $A$ 是正定的,则在 $A$ 的分解 $A=LDM^T$ 中, $D$ 的对角元均大于 $0$ .然而,仅仅存在 $LDM^T$ 分解仍不足以保证分解算法的稳定性.

以下着重介绍对称正定阵的Cholesky分解.

定理3(Cholesky分解) 如果 $A\in \mathbb{R}^{n\times n}$ 是对称正定的,则存在唯一的一个对角元大于 $0$ 的下三角阵 $G\in \mathbb{R}^{n\times n}$ ,满足 $A=GG^T$ .

证明由 $LDL^T$ 分解及 $A$ 的正定性容易得到 $D$ 对角元大于 $0$ ,从而得到本定理. □

我们下面导出一个含有大量Gaxpy运算的Cholesky分解算法.比较等式 $A=GG^T$ 的第 $j$ 列得到 $\begin{equation*} A(:,j)=\sum\limits_{k=1}^nG(j,k)G(:,k), \end{equation*}$ 也就是说, $\begin{equation*} G(j,j)G(:,j)=A(:,j)-\sum\limits_{k=1}^{j-1}G(j,k)G(:,k)\equiv v. \end{equation*}$ 如果 $G$ 的第 $j-1$ 列已知,则可计算出 $v$ .由上式中各元素间的相等关系推出 $\begin{equation*} G(j:n,j)=v(j:n)/\sqrt{v(j)}. \end{equation*}$ 于是得到以下算法:

算法2(Cholesky分解:基于Gaxpy运算) 给出对称正定阵 $A\in \mathbb{R}^{n\times n}$ ,本算法计算出一个下三角阵 $G\in \mathbb{R}^{n\times n}$ 满足 $A=GG^T$ .对所有 $i\geqslant j$ , $G(i,j)$ 覆盖 $A(i,j)$ .

$\begin{align*} &A(1:n,1)\leftarrow A(1:n,1)/\sqrt{A(1,1)}.\\ &\text{对}j>1,\text{通过以下步骤计算}G(j:n,j)\text{并覆盖}A(j:n,j):\\ &A(j:n,j)\leftarrow A(j:n,j)-A(j:n,1:j-1)A(j,1:j-1)^T. \end{align*}$

本算法共需 $n^3/2$ 个flop.这一算法在 $A$ 正定时相当稳定.

Vandemonde方程组

Vandemonde方程组与下一节将要讲到的正定对称的Toeplitz阵问题是两类重要的能在 $O(n^2)$ 的计算量内稳定求解的问题.

设 $x(0:n)\in \mathbb{R}^{n+1}$ ,形如 $\begin{equation*} V=V(x_0,\cdots,x_n)= \begin{bmatrix} 1& 1& \cdots& 1\\ x_0& x_1& \cdots& x_n\\ \vdots& \vdots& & \vdots\\ x_0^n& x_1^n& \cdots & x_n^n \end{bmatrix} \end{equation*}$ 的矩阵 $V\in \mathbb{R}^{(n+1)\times(n+1)}$ 被称作Vandemonde矩阵.下面我们讨论线性方程组 $V^Ta=f$ 及 $Vx=b$ 的求解.

Vandemonde矩阵与多项式插值有密切的联系.首先注意到,若有 $V^Ta=f$ ,并定义 $p(x)=\sum\limits_{j=0}^na_jx^j$ ,则对于 $i=0:n$ , $p(x_i)=f_i$ .故若 $x_i$ 是互异的,则可通过构造插值多项式 $p(x)$ 求解.

第一步是计算插值多项式 $p$ 的Newton表达式: $\begin{equation*} p(x)=\sum\limits_{k=0}^nc_k \left(\prod\limits_{i=0}^{k-1} (x-x_i)\right), \end{equation*}$ 常数 $c_k$ 是差商,可按以下步骤确定: $\begin{align*} &c(0:n)\leftarrow f(0:n).\\ &\text{for } k=0:n-1\\ &~~~~\text{for } i=n:-1:k+1\\ &~~~~~~~~c_i\leftarrow (c_i-c_{i-1})/(x_i-x_{i-k-1}).\\ &~~~~\text{end}\\ &\text{end} \end{align*}$ 下一步是由 $c(0:n)$ 产生 $a(0:n)$ .通过迭代 $\begin{align*} &p_n(x)\leftarrow c_n,\\ &\text{for } k=n-1:-1:0\\ &~~~~p_k(x)\leftarrow c_k+(x-x_k)p_{k+1}(x).\\ &\text{end} \end{align*}$ 来定义多项式 $p_k(x) (x=n:-1:0)$ .可求得 $p(x)\equiv p_0(x)$ .于是系数 $a_i\equiv a_i^{(0)}$ 可如下求得 $\begin{align*} &a(0:n)\leftarrow c(0:n),\\ &\text{for } k=n-1:-1:0\\ &~~~~\text{for } i=k:n-1\\ &~~~~~~~~a_i\leftarrow a_i-x_ka_{i+1}.\\ &~~~~\text{end}\\ &\text{end} \end{align*}$ 以上两步骤合并即得到求解 $V^Ta=f$ 的算法,共需 $5n^2/2$ 个flop.

该算法事实上即多项式插值问题的一般算法.它事实上利用了如下分解 $\begin{equation*} (V^T)^{-1}=L^TU^T= (L_0(x_0)^T \cdots L_{n-1}(x_{n-1})^T) (D_{n-1}^{-1}L_{n-1} (1)\cdots D_0^{-1}L_0 (1)), \end{equation*}$ 其中 $\begin{displaymath} L_k (\alpha)= \left[\begin{array}{c|c} I_k& O\\ \hline\\ O& \begin{matrix} 1& \cdots& & O\\ -\alpha& \ddots\\ \vdots& \ddots& \ddots\\ 0& \cdots&\ -\alpha& 1 \end{matrix} \end{array}\right], \end{displaymath} \begin{equation*} D_k=\mathrm{diag} (\underbrace{1,\cdots. 1}_{k+1},x_{k+1}-x_0,\cdots,x_n-x_{n-k-1}). \end{equation*}$ 于是对于方程组 $Vz=b$ ,只要利用 $V^{-1}=UL$ 即可得到其三角阵分解,从而求解方程组.

对称正定Toeplitz阵及相关问题

设 $T\in \mathbb{R}^{n\times n}$ ,存在 $r_{-n+1},\cdots,r_0,\cdots,r_{n-1}$ 对所有 $i$ 和 $j$ 满足 $a_{ij}=r_{j-i}$ ,则 $T$ 称为Toeplitz阵.Toeplitz阵属于一个更大的反向对称(persymmetric)矩阵类 (即沿反对角线进行"转置"不改变).容易验证,Toeplitz阵及非奇异Toeplitz阵的逆也是反向对称的.本节将利用这一特性给出对称正定Toeplitz阵的 $O(n^2)$ 的算法.不失一般性,以下均设 $r_0=1$ .

Yule-Walker问题的Durbin算法

Yule-Walker问题是系数矩阵为Toeplitz阵 $T_n$ ,非齐次项 $-r=- (r_1,\cdots,r_n)^T$ 的线性方程组 $T_ny=-r$ .可以看到,这种问题具有很大的特殊性,即它的非齐次项与系数矩阵有密切关系,但它是进一步讨论的基础.设我们已经求解了 $k$ 阶的Yule-Walker方程组 $T_ky=-r=- (r_1,\cdots,r_k)$ .下面给出如何在 $O(k)$ 内求解 $(k+1)$ 阶Yule-Walker方程组,其中 $E_k$ 为 $k$ 阶反向置换阵: $\begin{equation*} E_k= \begin{bmatrix} & & 1\\ & \adots& \\ 1& & \end{bmatrix} \end{equation*}$ $\begin{equation*} \begin{bmatrix} T_k& E_kr\\ r^TE_k& 1 \end{bmatrix} \begin{bmatrix} z\\ \alpha \end{bmatrix} =- \begin{bmatrix} r\\ r_{k+1} \end{bmatrix}. \end{equation*}$ 注意到 $\begin{gather*} z=T_k^{-1} (-r-\alpha E_kr)=y-\alpha T_k^{-1}E_kr,\\ \alpha=-r_{k+1}-r^TE_kz, \end{gather*}$ 由于 $T_kE_k=E_kT_k$ ,故有 $\begin{equation*} z=y-\alpha E_kT_k^{-1}r=y+\alpha E_ky. \end{equation*}$ 代入 $\alpha$ 表达式 $\begin{equation*} \alpha=-r_k-r^T (y+\alpha E_ky)=- (r_{k+1}+r^TE_ky)/(1+r^Ty), \end{equation*}$ 其中 $\beta\equiv 1+r^Ty$ 是大于 $0$ 的,因Toeplitz阵 $T_{k+1}$ 为正定的: $\begin{equation*} \begin{bmatrix} I& E_ky\\ 0& 1 \end{bmatrix} ^T \begin{bmatrix} T_k& E_kr\\ r^TE_k& 1 \end{bmatrix} \begin{bmatrix} I& E_ky\\ 0& 1 \end{bmatrix} = \begin{bmatrix} T_k& \\ & 1+r^Ty \end{bmatrix}. \end{equation*}$ 进一步利用 $\beta^{(k)}=1+ (r^{(k)})^Ty^{(k)}= (1-\alpha_{k-1}^2)\beta_{k-1}$ ,即得到完整的Durbin算法,需 $2n^2$ 个flop.

一般右端向量的Levinson算法

可以用类似方法求解右端为任意向量的Toeplitz方程组 $T_nx=b$ .

假设已解出 $k$ 阶Toeplitz方程组 $T_kx=b= (b_1,\cdots,b_k)^T$ ,现需求解 $\begin{equation*} \begin{bmatrix} T_k& E_kr\\ r^TE_k& 1 \end{bmatrix} \begin{bmatrix} v\\ \mu \end{bmatrix} = \begin{bmatrix} b\\ b_{k+1} \end{bmatrix}, \end{equation*}$ 同时假定 $k$ 阶Yule-Walker方程组 $T_ky=-r$ 的解也已得到.由 $T_kv+\mu I_kr=b$ 可知 $\begin{gather*} v=T_k^{-1} (b-\mu E_kr)=x-\mu T_k^{-1}E_kr=x+\mu E_ky,\\ \mu=b_{k+1}-r^TE_kv= (b_{k+1}-r^TE_kx)/ (1+r^Ty). \end{gather*}$ 这样,我们可在 $O(k)$ 个flop内完成 $k+1$ 阶方程组的求解.全部算法共需 $4n^2$ 个flop.

正定对称Toeplitz阵求逆Trench算法

对称正定的Toeplitz阵 $T_n$ 的最令人吃惊的性质之一是它的逆可以在 $O(n^2)$ 个flop内算出.为了得到这个算法,我们将 $T_n^{-1}$ 做如下划分: $\begin{equation*} T_n^{-1}= \begin{bmatrix} A& Er\\ r^TE& 1 \end{bmatrix} ^{-1}= \begin{bmatrix} B& v\\ v^T& \gamma \end{bmatrix}, \end{equation*}$ 由 $\begin{equation*} \begin{bmatrix} A& Er\\ r^TE 1 \end{bmatrix} \begin{bmatrix} v\\ \gamma \end{bmatrix} = \begin{bmatrix} 0\\ 1 \end{bmatrix} \end{equation*}$ 得到 $\begin{gather*} Av=-\gamma Er,\\ \gamma=1-r^TEv. \end{gather*}$ 若 $y$ 是 $n-1$ 阶Yule-Walker方程组 $Ay=-r$ 的解,则 $\begin{gather*} \gamma=1/(1+r^Ty),\\ v=\gamma Ey. \end{gather*}$ 于是我们得到了 $T_n^{-1}$ 的最后一列.

因为 $AB+Erv^T=I_{n-1}$ ,故 $B=A^{-1}- (A^{-1}Er)v^T=A^{-1}+vv^T/\gamma$ .由于 $A=T_{n-1}$ 是非奇异的Toeplitz阵,其逆是反向对称的,于是 $\begin{equation*} \begin{split} b_{ij}&= (A^{-1})_{ij}+\frac{v_iv_j}{\gamma}\\ &= (A^{-1})_{n-j,n-i}+\frac{v_iv_j}{\gamma}\\ &= b_{n-j,n-i}-\frac{v_{n-j}v_{n-i}}{\gamma}+\frac{v_iv_j}{\gamma}\\ &=b_{n-j,n-i}-\frac{1}{\gamma} (v_iv_j-v_{n-j}v_{n-i}), \end{split} \end{equation*}$ 于是我们可以"从外向里"地确定 $B$ ,以 $6$ 阶矩阵为例形象地描述如下:

假定 $T_n^{-1}$ 的最后一行和最后一列已知 $\begin{equation*} T_n^{-1}= \begin{bmatrix} *& *& *& *& *& k\\ *& *& *& *& *& k\\ *& *& *& *& *& k\\ *& *& *& *& *& k\\ *& *& *& *& *& k\\ k& k& k& k& k& k \end{bmatrix}, \end{equation*}$ 这里 $*$ 和 $k$ 分别代表未知元素和已知元素.交替利用 $T_{n-1}^{-1}$ 的反向对称性和前述递推式,可求出 $(n-1)\times(n-1)$ 的顺序子块 $B$ 如下:

$\begin{gather*} \underrightarrow{\text{反向对称}} \begin{bmatrix} k& k& k& k& k& k\\ k& *& *& *& *& k\\ k& *& *& *& *& k\\ k& *& *& *& *& k\\ k& *& *& *& *& k\\ k& k& k& k& k& k \end{bmatrix} \underrightarrow{\text{递推关系}} \begin{bmatrix} k& k& k& k& k& k\\ k& *& *& *& k& k\\ k& *& *& *& k& k\\ k& *& *& *& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k \end{bmatrix}\\ \underrightarrow{\text{反向对称}} \begin{bmatrix} k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& *& *& k& k\\ k& k& *& *& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k \end{bmatrix} \underrightarrow{\text{递推关系}} \begin{bmatrix} k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& *& k& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k \end{bmatrix}\\ \underrightarrow{\text{反向对称}} \begin{bmatrix} k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k\\ k& k& k& k& k& k \end{bmatrix}. \end{gather*}$ 当然,当求解一个既对称又反对称的矩阵时,只需求解矩阵的"上楔形"即可,即 $\begin{equation*} \begin{matrix} *& *& *& *& *& *\\ & *& *& *& *\\ & & *& * \end{matrix} (n=6). \end{equation*}$ 由以上分析,不难得到对于 $n$ 阶正定对称Toeplitz阵的求逆算法,称为Trench算法.它需要 $13n^2/4$ 个flop.