隐藏目录
对称$QR$算法
Jacobi方法
对称三对角矩阵的特征值方法
Sturm序列与二分法
分而治之方法
计算SVD

由于对称矩阵的结构性质,其特征值问题在数值线性代数中有许多漂亮的算法.除$QR$迭代可以应用于对称特征值问题外,还有Jacobi迭代以及与SVD直接相关的三对角方法.

回到顶部对称$QR$算法

将$QR$迭代的思想结合矩阵的对称性质可以得到对称$QR$算法.其导出思路可以简述如下:

用于求最大特征值及特征向量的幂法和用于求不变子空间的$QR$迭代分别以$\left|\frac{\lambda_2}{\lambda_1}\right|^{2k}$和$\left|\frac{\lambda_{r+1}}{\lambda_r}\right|^k$的方式收敛到相应特征向量和不变子空间.这引导我们考虑$QR$迭代,这相当于$r=n$的正交迭代.

我们考察对称$QR$迭代的加速算法.首先通过一系列Householder变换$U$使得$U^TAU=T$是三对角阵,从而将之后每步迭代的运算量降至$O (n)$个flop,而三对角化这一步骤约需$4n^3/3$个flop.其次,应用位移思想,可使收敛到对角形的速率为立方收敛.我们着重分析位移思想.

若$s$是一个近似特征值,则我们期望用位移$s$进行一次$QR$迭代后,$(n,n-1)$元会很小.若 \begin{equation*} T= \begin{bmatrix} a_1& b_1\\ b_1& a_2& \ddots\\ & \ddots& \ddots& b_{n-1}\\ & & b_{n-1}& a_n \end{bmatrix}, \end{equation*} 由Gerschgorin定理,位移的一个合理选择是$\mu=a_n$,而一个更有效的选择是用$\left[\begin{smallmatrix}a_{n-1}& b_{n-1}\\ b_{n-1}& a_n \end{smallmatrix}\right]$靠近$a_n$的特征值作位移,这称为Wilkinson位移,它由 \begin{equation*} \mu=a_n+d-\mathrm{sign} (d)\sqrt{d^2+b_{n-1}^2} \end{equation*} 给出,其中$d\equiv (a_{n-1}-a_n)/2$.Wilkinson[1]证明了以上两种位移策略都是立方收敛,并给出了理由为什么更偏好后者.

应用隐式Q定理,不必显式形成矩阵$T-\mu I$就能实现从$T$到$T^+\equiv RU+\mu I=U^TTU$的变换,这在$\mu$远大于某个$a_j$时有优点.

首先取 \begin{equation*} G_1=G(1,2,\theta)=\left[ \begin{array}{cc|c} c& s& \\ -s& c& \\\hline & & I_{n-2} \end{array}\right], \end{equation*} 使得 \begin{equation*} \begin{bmatrix} c& s\\ -s& c \end{bmatrix}^T \begin{bmatrix} a_1-\mu\\ b_1 \end{bmatrix}= \begin{bmatrix} x\\ 0 \end{bmatrix}, \end{equation*} 则 \begin{equation*} G_1e_1=Ue_1. \end{equation*} 再计算$G_2,\cdots,G_{n-1}$使得$Z\equiv G_1\cdots G_{n-1}$满足$Ze_1=G_1e_1=Ue_1$且$Z^TTZ$是三对角阵,就可以应用隐式Q定理.而实际中,我们使 \begin{equation*} G_i=G (i,i+1,\theta_i),(i=2:n-1), \end{equation*} 则$Z$和$Q$的第一列相等且$G_i$可用于将多余的非零元素逐出矩阵$G_1^TTG_1$.这样就完成了$T$到$T^+=Z^TTZ$的变换,从而完成了对称$QR$算法的一个步骤.

回到顶部Jacobi方法

求解对称特征值问题的Jacobi迭代法因易于并行化而引起人们的注意.它通过进行一系列正交相似变换 \begin{equation*} A\leftarrow Q^TAQ, \end{equation*} 逐步减小 \begin{equation*} \text{off} (A)\equiv \sqrt{\sum\limits_{i=1}^n\sum\limits_{j\neq i}a_{ij}^2}. \end{equation*} 我们可以形象地称$\text{off}(A)$为非对角元素的Frobenius范数.实现它的工具是Jacobi-Givens旋转变换.

Jacobi方法的基本步骤:

  1. 选择指标$(p,q) (1\leqslant p<q\leqslant n)$,在经典Jacobi算法中选择$(p,q)$使$a_{pq}^2$最大.
  2. 计算一个余弦-正弦对$(c,s)$,使得 \begin{equation*} \begin{bmatrix} c& s\\ -s& c \end{bmatrix}^T \begin{bmatrix} a_{pp}& a_{pq}\\ a_{qp}& a_{qq} \end{bmatrix} \begin{bmatrix} c& s\\ -s& c \end{bmatrix}= \begin{bmatrix} b_{pp}& 0\\ 0& b_{qq} \end{bmatrix} \end{equation*} 为对角阵.
  3. 令$J=J(p,q,\theta)$,计算$B=J^TAJ$并覆盖$A$.

容易得到,经过一次Jacobi步骤后, \begin{equation*} \text{off} (B)^2=\text{off} (A)^2-2a_{pq}^2. \end{equation*} 在此意义下,每个Jacobi步骤后$A$更"靠近"对角形.

经典Jacobi方法(选取$a_{pq}^2$最大)的每次校正要做$O (n)$次运算,但选取最大元素却需$O (n^2)$次比较.解决此矛盾的途径是将变换顺序固定下来,如逐行进行"扫描",这称为行循环方法,它是二次收敛的.作为对比，$QR$算法是立方收敛的.

如前所述,Jacobi方法的优势在于其并行本性.事实上,Jacobi旋转$J(p,q,\theta)$只作用于$p$,$q$行和$p$,$q$列.这样,每次将矩阵的全部行与列划分为若干互不冲突的组,每组的运算可以得到并行的处理.各族变换并行处理完成后,变换行与列的划分(类似一次循环赛的过程),直至所有$(p,q)$对都得到处理.这种并行的思想也同样适用于分块的Jacobi迭代.

回到顶部对称三对角矩阵的特征值方法

回到顶部Sturm序列与二分法

令$T_r$表示$T $的$r $阶顺序主子阵,定义其特征多项式$p_r (x)=\det (T_r-xI)$,则有以下递推公式(设$p_0 (x)\equiv 1$): \begin{equation*} p_r (x)= (a_r-x)p_{r-1} (x)-b_{r-1} ^2p_{r-2}(x). \end{equation*} 因通过$O (n)$次运算即可计算$p_n (x)$的值,故利用二分法找它的根是可行的.显然,二分迭代是线性收敛,即每步将误差减半.

对于不可约三对角阵,有如下经典结果,其证明可参见[2].

定理1(Sturm序列性质) 对于不可约三对角阵$T$,$T_{r-1}$的特征值严格隔离$T_r$的特征值: \begin{equation*} \lambda_r (T_r)<\lambda_{r-1} (T_{r-1})<\lambda_{r-1} (T_r)<\cdots<\lambda_2 (T_r)<\lambda_1 (T_{r-1})<\lambda_1 (T_r). \end{equation*}

此外,若$a(\lambda)$表示在序列$\{p_0 (\lambda),p_1 (\lambda),\cdots,p_n (\lambda)\}$中符号改变的个数,则$a(\lambda)$等于$T$的比$\lambda$小的特征值个数.其中$p_r (\lambda)$如上定义,且约定若$p_(\lambda)=0$,则$p_r(\lambda)$与$p_{r-1} (\lambda)$反号.

结合Gerschgorin定理可知$\lambda_k (T)\in [y,z]$.其中$y\equiv \min\limits_{1\leqslant i\leqslant n} (a_i-|b_i|-|b_{i-1}|)$,$z\equiv \max\limits_{1\leqslant i\leqslant n} (a_i+|b_i|+|b_{i-1}|)$.于是我们可以此为初始值,进行二分法迭代,每次迭代的判断条件为$a (x)\geqslant n-k$,则可求出$T$的第$k $大的特征值.

与$QR$迭代相比,对分法的优势在于无论特征值大小,计算值都具有小的相对误差.

回到顶部分而治之方法

适合于并行处理的三对角阵特征值算法是以以下观察为基础的:设$n=2m$,定义 \begin{equation*} v= \begin{bmatrix} e_m^{(m)}\\ \theta e_1^{(m)} \end{bmatrix} \in \mathbb{R}^n, \end{equation*} 则 \begin{equation*} \tilde{T}=T-\rho vv^T \end{equation*} 除了"中间四个"元素 \begin{equation*} \tilde{T} (m:m+1,m:m+1)\equiv \begin{bmatrix} a_m-\rho& b_m-\rho\theta\\ b_m-\rho\theta& a_{m+1}-\rho\theta^2 \end{bmatrix}, \end{equation*} $\tilde{T}$与$T $相等.若取$\rho\theta=b_m$,则 \begin{equation*} T= \begin{bmatrix} T_1\\ & T_2 \end{bmatrix}+\rho vv^T, \end{equation*} 其中, \begin{align*} T_1 &\equiv \begin{bmatrix} a_1& b_1\\ b_1& a_2& \ddots\\ & \ddots& \ddots& b_{m-1}\\ & & b_{m-1}& a_m-\rho \end{bmatrix},\\ T_2 &\equiv \begin{bmatrix} a_{m+1}-\rho\theta^2& b_{m+1}\\ b_{m+1}& a_{m+2}& \ddots\\ & \ddots& \ddots& b_{n-1}\\ & & b_{n-1}& a_n \end{bmatrix}. \end{align*} 这样,我们就把$T$"撕"成了两片,并加上一个修正项.分别解决$T_1$,$T_2$的特征值问题 \begin{align*} Q_1^TT_1Q_1&=D_1,\\ Q_2^TT_2Q_2&=D_2 \end{align*} 后,只要快速求解一个对角矩阵向量外积修正的特征值问题 \begin{equation*} \begin{bmatrix} D_1\\ & D_2 \end{bmatrix}+\rho (U^Tv) (U^Tv)^T \end{equation*} 就可以了.

对角矩阵向量外积修正的特征值问题的主要计算依赖于以下结果:

定理2 假定$D\in \mathrm{R}^{n\times n}$具有性质$d_1>\cdots>d_n$.又设$\rho\neq 0$且$z\in \mathrm{R}^n$没有零分量.若$(D+\rho zz^T)v=\lambda v,v\neq 0$,则$z^Tv\neq0$且$D-\lambda I$非奇异.
定理3 设$D=\mathrm{diag} (d_1,\cdots,d_n)\in \mathrm{R}^{n\times n}$且对角元满足$d_1>\cdots>d_n$.假定$\rho\neq 0$且$z\in \mathrm{R}^n$没有零分量,若正交阵$V\in \mathrm{R}^{n\times n}$使得 \begin{equation*} V^T (D+\rho zz^T)=\mathrm{diag} (\lambda_1,\cdots,\lambda_n), \end{equation*} 其中$\lambda_1\geqslant\cdots\geqslant\lambda_n$,且$V=[v_1,\cdots,v_n]$.则
  1. $\lambda_i$是 \begin{equation*} f (\lambda)=1+\rho z^T (D-\lambda I)^{-1}z \end{equation*} 的$n$个零解.
  2. 当$\rho>0$,则 \begin{equation*} \lambda_1>d_1>\lambda_2>\cdots>d_n. \end{equation*} 当$\rho<0$,则 \begin{equation*} d_1>\lambda_1>d_2>\cdots>\lambda_n. \end{equation*}
  3. 特殊向量$v_i$是$(D-v_i I)^{-1}z$的倍数.

上述定理表明要计算$V$,我们应首先用Newton型算法找出$f$的根$\lambda_1,\cdots,\lambda_n$,然后对$i=1:n$,通过正规化向量$(D-\lambda_iI)^{-1}z$来计算$V$的列向量.

以下考察有重复的$d_i$或零分量$z_i$的情形,我们对以下定理给出一个构造性的证明,从而解决问题:

定理4 若$D=\mathrm{diag} (d_1,\cdots,d_n)$且$z\in \mathrm{R}^n$,则存在一正交阵$V_1$使得 \begin{equation*} V_1^TDV_1=\mathrm{diag} (\mu_1,\cdots,\mu_n), \end{equation*} 且 \begin{equation*} w=V_Tz, \end{equation*} 其中 \begin{equation*} \mu_1>\mu_2>\cdots>\mu_r\geqslant\mu_{r+1}\geqslant\mu_n, \end{equation*} 其中,$w_i\neq 0,\forall i=1:r$而$w_i=0,\forall i=r+1:n$

证明
  1. 假设对某个$i<j$有$d_i=d_j$,令$J (i,j,\theta)$是$(i,j)$平面上的旋转变换,使$J (i,j,\theta)^Tz$的第$j $个分量为$0$.易知 \begin{equation*} J (i,j,\theta)^TDJ (i,j,\theta)=D. \end{equation*} 于是,若有一个重复的$d_i$,我们就能使$z$的一个分量化为$0$.
  2. 令$p $是$i $和$j $列互换的单位阵,可推出$P^TDP$是对角阵,这样我们就可将为零的$z_i$放在底部.重复以上步骤就可得到所需的标准结构.$V_1$是这些旋转阵之积. □

综合上述所述,我们就得到了分而治之算法.

回到顶部计算SVD

矩阵$A$的奇异值分解与对称阵$A^TA$,$AA^T$的Schur分解有着密切的联系.若$U^TAV=\mathrm{diag} (\sigma_1,\cdots,\sigma_n)$是$A\in \mathbb{R}^{m\times n},(m\geqslant n)$的SVD,则 \begin{align*} V^T (A^TA)V&=\mathrm{diag} (\sigma_1^2,\cdots,\sigma_n^2)\in \mathbb{R}^{n\times n},\\ U^T (AA^T)U&=\mathrm{diag} (\sigma_1^2,\cdots,\sigma_n^2;\underbrace{0,\cdots,0}_{m-n})\in \mathbb{R}^{n\times n}. \end{align*}

我们计算$A$的SVD的思路就是将其化为$A^TA$的Schur分解来进行.但如果显式形成$A^TA$,则可能造成信息的损失.Golub和Kahan (1965)[3]给出了基于隐式对称$QR$算法的SVD方法.

首先用Householder变化将$A$化为双对角形: \begin{equation*} U_B^AV_B= \begin{bmatrix} B\\ O \end{bmatrix}, \end{equation*} 其中 \begin{equation*} B= \begin{bmatrix} d_1& f_1& \cdots& O\\ & d_2& \ddots& \vdots\\ & & \ddots& f_{n-1}\\ & & & d_n \end{bmatrix}. \end{equation*} 接下来只要计算$B$的SVD.为此考虑三对角阵$T=B^TB$,应用$QR$算法:

  1. 计算矩阵 \begin{equation*} \begin{bmatrix} d_m^2+f_m^2& d_mf_m\\ d_mf_m& d_n^2+f_n^2 \end{bmatrix} \end{equation*} (其中$m\equiv n-1$)的靠近$d_n^2+f_n^2$的特征值$\lambda$.
  2. 计算$c_1=\cos\theta_1$和$s_1=\sin\theta_1$使得 \begin{equation*} \begin{bmatrix} c_1& s_1\\ -s_1& c_1 \end{bmatrix}^T \begin{bmatrix} d_1^2-\lambda\\ d_1f_1 \end{bmatrix}= \begin{bmatrix} *\\ 0 \end{bmatrix}. \end{equation*} 令$G_1=G(1,2,\theta)$.
  3. 计算Givens旋转$G_2,\cdots,G_{n-1}$,使得当$Q=G_1\cdots G_{n-1}$时,$Q^TTQ$为三对角阵,且$Qe_1=G_1e_1$.

正如之前指出的,以上算法需显式计算$B^TB$,这可能导致算法的不稳定.若我们把Givens变化$G_1$直接作用到$B $上,在决定Givens旋转$U_1,V_2,U_2,\cdots,V_{n-1},U_{n-1}$将$B$恢复为双对角阵,得到 \begin{equation*} \bar{B}= (U_{n-1}^T\cdots U_1^T)B (G_1V_2\cdots V_{n-1})=\bar{U}^TB\bar{V}. \end{equation*} 由$V_i$的形式可知 \begin{equation*} \bar{V}e_1=Qe_1. \end{equation*} 从而由隐式Q定理断言$\bar{V}$和$Q$本质相同.于是我们隐式地完成了$T=B^TB$到$\bar{T}=\bar{B}^T\bar{B}$的过程.

由于以上本质上是对称$QR$迭代的方法,所以其收敛是三次的.实际中我们通过检测$B$的小次对角元来决定何时将问题降阶或终止迭代.

Chan (1982)[4]对Golub和Kahan的算法进行了改进.通过合理安排运算顺序,在$m>n$时可以有效地减少运算量.

  1. 形成双对角矩阵时,先通过Householder变换,将$A$上三角化: \begin{equation*} L^TA= \begin{bmatrix} R\\ O \end{bmatrix} \begin{array}[c]{l} n\\ m-n \end{array}, \end{equation*} 其中$R$为上三角阵.随后利用$R$的上三角性质,利用Givens变换或快速Givens变换将其化为双对角形.由于R-双对角化时只需对较短的向量进行操作,而Givens变换较Householder变换可以减少加法的数目,特别在不需要求变换矩阵时可以有效地减少运算量.当$m\gg n$时尤其有效.
  2. 当求得$R=X\Sigma Y^T$之后,SVD的变换矩阵可如下求出: \begin{equation*} A=L \begin{bmatrix} X\\ O \end{bmatrix} \Sigma Y^T=L \begin{bmatrix} I\\ O \end{bmatrix} X\Sigma Y^T. \end{equation*}

Chan给出了算法复杂度的详细分析.他指出在$m\simeq 2n$时,R-双对角化几乎总是比Golub-Kahan算法更快速,特别在不需显式求出变换矩阵的情况下.

设我们已求得$m\times n$矩阵$A$的奇异值$\{\sigma_1,\cdots,\sigma_r\}$,来考虑如何求正交矩阵$U$,$V$,使$A=U^T\Sigma V$,其中$\Sigma=\mathrm{diag} (\sigma_1,\cdots,\sigma_r)$.

事实上,我们已求得$A^TA$的特征值$\sigma_1^2,\cdots,\sigma_r^2,0,\cdots,0$,要求正交矩阵$V$,其列向量是$A^TA$的特征向量,我们只要用逆迭代求解$A^TAy=\sigma^2y$ .实用中采用双对角矩阵$J$代替$A$,之后令$U$的相应列$x=\sigma^{-1}Jy$即可.这些算法都是稳定有效的.

参考文献

[1]Wilkinson J.H., Global Convergence of Tridiagonal QR Algorithm With Origin Shifts, Lin. Alg. and Its Applic. I (1968), 409-420.

[2]Wilkinson J. H., The Algebraic Eigenvalue Problem, Clarendon Press, Oxford, UK, 1965.

[3]Golub G.H. and Kahan W., Calculating the Singular Values and Pseudo-Inverse of a Matrix, SIAM J. Num. Anal. Ser. B 2 (1965), 205-224.

[4]Chan T., An Improved Algorithm for Computing the Singular Value Decomposition, ACM Trans. Math. Soft. 8 (1982), no.1, 72-83.

最后更新： 2009年01月28日 22:41:14 返回主页 / 更新记录 / 许可协议 Copyright © 2007,2008 maTHμ Project Group