把一个向量反复乘上同一矩阵，会发生什么？取
A = [ 2 1 1 2 ] , x 0 = ( 1 , 0 ) T . A=\begin{bmatrix}2&1\\1&2\end{bmatrix},\qquad x_0=(1,0)^T. A = [ 2 1 1 2 ] , x 0 = ( 1 , 0 ) T . 第一次得到 ( 2 , 1 ) T (2,1)^T ( 2 , 1 ) T ( 5 , 4 ) T (5,4)^T ( 5 , 4 ) T ( 1 , 1 ) T (1,1)^T ( 1 , 1 ) T
这个现象给出了一条求特征值的路：不展开高次行列式，而是观察矩阵怎样反复改变向量。只是，方向有没有收敛、找到的是哪个特征值、剩下的误差有多大，都需要各自检查。
幂法留下的是最大模方向 特征向量 v ≠ 0 v\ne0 v = 0 A v = λ v Av=\lambda v A v = λ v A A A v 1 , … , v n v_1,\ldots,v_n v 1 , … , v n
x 0 = ∑ j = 1 n c j v j . x_0=\sum_{j=1}^n c_jv_j. x 0 = j = 1 ∑ n c j v j . 连续乘 A A A
A k x 0 = ∑ j = 1 n c j λ j k v j . A^kx_0=\sum_{j=1}^nc_j\lambda_j^kv_j. A k x 0 = j = 1 ∑ n c j λ j k v j . 假设
∣ λ 1 ∣ > ∣ λ 2 ∣ ≥ ⋯ ≥ ∣ λ n ∣ , c 1 ≠ 0. |\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_n|,\qquad c_1\ne0. ∣ λ 1 ∣ > ∣ λ 2 ∣ ≥ ⋯ ≥ ∣ λ n ∣ , c 1 = 0. 把第一项的系数提出来：
A k x 0 = c 1 λ 1 k [ v 1 + ∑ j = 2 n c j c 1 ( λ j λ 1 ) k v j ] . A^kx_0=c_1\lambda_1^k
\left[v_1+\sum_{j=2}^n\frac{c_j}{c_1}
\left(\frac{\lambda_j}{\lambda_1}\right)^kv_j\right]. A k x 0 = c 1 λ 1 k [ v 1 + j = 2 ∑ n c 1 c j ( λ 1 λ j ) k v j ] . 括号里其余系数的模都趋于零。因此，忽略整体倍数后，方向趋向 v 1 v_1 v 1 A A A A k A^k A k
这个论证里的条件各有用途。可对角化保证上述展开成立，严格的模间隔让其他方向衰减，c 1 ≠ 0 c_1\ne0 c 1 = 0 “最大模”不等于“代数值最大”。 如果特征值为 − 5 , 2 , 1 -5,2,1 − 5 , 2 , 1 − 5 -5 − 5
方向逐渐靠向占优势的特征方向。这里用未归一化箭头示意方向变化，长度、夹角和背景椭圆不对应某个指定矩阵的实际迭代；数值过程在交互中逐步计算。
每轮归一化，避免数值膨胀 直接保存 A k x 0 A^kx_0 A k x 0
y k = A x k , x k + 1 = y k ∥ y k ∥ ∞ . y_k=Ax_k,\qquad x_{k+1}=\frac{y_k}{\|y_k\|_\infty}. y k = A x k , x k + 1 = ∥ y k ∥ ∞ y k . 这里采用无穷范数，是为了让手算的最大分量保持为 1；用 2-范数同样可以。除以一个非零标量不改变方向，所以前面的论证仍然适用。
如果 y k = 0 y_k=0 y k = 0 x k x_k x k
开头的例子变成
x 1 = ( 1 , 1 / 2 ) T , A x 1 = ( 5 / 2 , 2 ) T , x 2 = ( 1 , 4 / 5 ) T . x_1=(1,1/2)^T,\qquad
Ax_1=(5/2,2)^T,\qquad
x_2=(1,4/5)^T. x 1 = ( 1 , 1/2 ) T , A x 1 = ( 5/2 , 2 ) T , x 2 = ( 1 , 4/5 ) T . 两个精确特征方向分别为
v 1 = ( 1 , 1 ) T 2 , λ 1 = 3 ; v 2 = ( 1 , − 1 ) T 2 , λ 2 = 1. v_1=\frac{(1,1)^T}{\sqrt2},\quad \lambda_1=3;
\qquad
v_2=\frac{(1,-1)^T}{\sqrt2},\quad \lambda_2=1. v 1 = 2 ( 1 , 1 ) T , λ 1 = 3 ; v 2 = 2 ( 1 , − 1 ) T , λ 2 = 1. x 0 x_0 x 0 v 1 , v 2 v_1,v_2 v 1 , v 2 v 1 v_1 v 1 θ k \theta_k θ k
tan θ k = 3 − k . \tan\theta_k=3^{-k}. tan θ k = 3 − k . 所以这里每轮把角度的正切缩小到原来的三分之一。一般情况下，∣ λ 2 / λ 1 ∣ |\lambda_2/\lambda_1| ∣ λ 2 / λ 1 ∣
看起来不动，未必完成了目标 对于 diag ( 3 , 1 ) \operatorname{diag}(3,1) diag ( 3 , 1 ) ( 0 , 1 ) T (0,1)^T ( 0 , 1 ) T
对于 diag ( 1 , − 1 ) \operatorname{diag}(1,-1) diag ( 1 , − 1 ) ( 1 , 1 ) T (1,1)^T ( 1 , 1 ) T ( 1 , 1 ) T (1,1)^T ( 1 , 1 ) T ( 1 , − 1 ) T (1,-1)^T ( 1 , − 1 ) T u , v u,v u , v
min { ∥ u − v ∥ 2 , ∥ u + v ∥ 2 } , \min\{\|u-v\|_2,\ \|u+v\|_2\}, min { ∥ u − v ∥ 2 , ∥ u + v ∥ 2 } , 避免把同一直线的相反朝向误判为不收敛。这个方向差仍不能替代原方程残差。
在实验里比较谱比相差很大的两种矩阵。运行之前，试着判断哪一种会更快；再切换到缺少主方向的初值和等模例子，观察为什么“多迭代几次”解决不了它们。表中的单位向量、Rayleigh 商与残差必须来自同一个迭代时刻。
用 Rayleigh 商估计特征值 手里已有非零近似向量 x x x μ \mu μ
∥ A x − μ x ∥ 2 2 = ∥ A x ∥ 2 2 − 2 μ x T A x + μ 2 x T x \|Ax-\mu x\|_2^2
=\|Ax\|_2^2-2\mu x^TAx+\mu^2x^Tx ∥ A x − μx ∥ 2 2 = ∥ A x ∥ 2 2 − 2 μ x T A x + μ 2 x T x 看成关于实数 μ \mu μ
ρ ( x ) = x T A x x T x . \rho(x)=\frac{x^TAx}{x^Tx}. ρ ( x ) = x T x x T A x . 这就是 Rayleigh 商。它不随 x x x x x x
对于实对称矩阵，可以取一组标准正交特征向量。这个事实也有一个直观的证明路线：二次型 u T A u u^TAu u T A u u u u w w w 2 w T A u 2w^TAu 2 w T A u A u Au A u u u u A A A u T A w = ( A u ) T w = 0 u^TAw=(Au)^Tw=0 u T A w = ( A u ) T w = 0
将单位向量展开为 u = ∑ j a j v j u=\sum_j a_jv_j u = ∑ j a j v j
∑ j a j 2 = 1 , ρ ( u ) = ∑ j λ j a j 2 . \sum_j a_j^2=1,\qquad
\rho(u)=\sum_j\lambda_j a_j^2. j ∑ a j 2 = 1 , ρ ( u ) = j ∑ λ j a j 2 . 所以对称情形的 Rayleigh 商是特征值的加权平均，必定落在最小与最大特征值之间。若 θ \theta θ u u u v 1 v_1 v 1
ρ ( u ) − λ 1 = ∑ j ≠ 1 ( λ j − λ 1 ) a j 2 , \rho(u)-\lambda_1
=\sum_{j\ne1}(\lambda_j-\lambda_1)a_j^2, ρ ( u ) − λ 1 = j = 1 ∑ ( λ j − λ 1 ) a j 2 , 从而
∣ ρ ( u ) − λ 1 ∣ ≤ max j ≠ 1 ∣ λ j − λ 1 ∣ sin 2 θ . |\rho(u)-\lambda_1|
\le \max_{j\ne1}|\lambda_j-\lambda_1|\sin^2\theta. ∣ ρ ( u ) − λ 1 ∣ ≤ j = 1 max ∣ λ j − λ 1 ∣ sin 2 θ . 这解释了一个常见现象：特征值读数已经很准，向量方向还没有同样准确。此处的平方误差依赖对称性，不能直接照搬到任意非对称矩阵。
把第二轮结果代回去 开头第二轮的向量是 x 2 = ( 1 , 4 / 5 ) T x_2=(1,4/5)^T x 2 = ( 1 , 4/5 ) T
A x 2 = ( 14 / 5 , 13 / 5 ) T , x 2 T A x 2 = 122 25 , x 2 T x 2 = 41 25 . Ax_2=(14/5,13/5)^T,\qquad
x_2^TAx_2=\frac{122}{25},\qquad
x_2^Tx_2=\frac{41}{25}. A x 2 = ( 14/5 , 13/5 ) T , x 2 T A x 2 = 25 122 , x 2 T x 2 = 25 41 . 因此
ρ ( x 2 ) = 122 41 ≈ 2.97561. \rho(x_2)=\frac{122}{41}\approx2.97561. ρ ( x 2 ) = 41 122 ≈ 2.97561. 它小于 3，与加权平均的范围一致。为了让残差不受向量长度影响，将 x 2 x_2 x 2
u 2 = ( 5 , 4 ) T 41 . u_2=\frac{(5,4)^T}{\sqrt{41}}. u 2 = 41 ( 5 , 4 ) T . 本章的特征残差约定为 r = A u − ρ u r=Au-\rho u r = A u − ρ u
r = 9 41 41 ( − 4 , 5 ) T , ∥ r ∥ 2 = 9 41 . r=\frac{9}{41\sqrt{41}}(-4,5)^T,\qquad
\|r\|_2=\frac9{41}. r = 41 41 9 ( − 4 , 5 ) T , ∥ r ∥ 2 = 41 9 . 特征值误差是 3 − ρ = 1 / 41 3-\rho=1/41 3 − ρ = 1/41 tan θ 2 = 1 / 9 \tan\theta_2=1/9 tan θ 2 = 1/9 sin 2 θ 2 = 1 / 82 \sin^2\theta_2=1/82 sin 2 θ 2 = 1/82 3 − ρ = 2 sin 2 θ 2 3-\rho=2\sin^2\theta_2 3 − ρ = 2 sin 2 θ 2
同一个单位向量既给出 Rayleigh 商，也给出特征残差。主例第二轮的值误差为 1/41，残差长度为 9/41，两项检查不能互相替代。
残差能证明多少 下文都让 u u u ∥ u ∥ 2 \|u\|_2 ∥ u ∥ 2
令 r = A u − ρ u r=Au-\rho u r = A u − ρ u
E = − r u T . E=-ru^T. E = − r u T . 由于 u T u = 1 u^Tu=1 u T u = 1
( A + E ) u = A u − r = ρ u , ∥ E ∥ 2 = ∥ r ∥ 2 . (A+E)u=Au-r=\rho u,\qquad \|E\|_2=\|r\|_2. ( A + E ) u = A u − r = ρ u , ∥ E ∥ 2 = ∥ r ∥ 2 . 因此，小残差保证这对数值是一个邻近矩阵的精确特征对。这是后向解释，适用于这里的任意实矩阵。若 A A A ρ = ρ ( u ) \rho=\rho(u) ρ = ρ ( u ) u T r = 0 u^Tr=0 u T r = 0
E = − r u T − u r T . E=-ru^T-ur^T. E = − r u T − u r T . 它同样满足 E u = − r Eu=-r E u = − r u u u r / ∥ r ∥ 2 r/\|r\|_2 r /∥ r ∥ 2
[ 0 − ∥ r ∥ 2 − ∥ r ∥ 2 0 ] , \begin{bmatrix}0&-\|r\|_2\\-\|r\|_2&0\end{bmatrix}, [ 0 − ∥ r ∥ 2 − ∥ r ∥ 2 0 ] , 其他正交方向上为零，故其 2-范数仍为 ∥ r ∥ 2 \|r\|_2 ∥ r ∥ 2 r = 0 r=0 r = 0 E = 0 E=0 E = 0
对称矩阵：值误差与方向误差分开估计 使用刚才的正交特征基，
r = ∑ j ( λ j − ρ ) a j v j , ∥ r ∥ 2 2 = ∑ j ( λ j − ρ ) 2 a j 2 . r=\sum_j(\lambda_j-\rho)a_jv_j,\qquad
\|r\|_2^2=\sum_j(\lambda_j-\rho)^2a_j^2. r = j ∑ ( λ j − ρ ) a j v j , ∥ r ∥ 2 2 = j ∑ ( λ j − ρ ) 2 a j 2 . 令 d = min j ∣ λ j − ρ ∣ d=\min_j|\lambda_j-\rho| d = min j ∣ λ j − ρ ∣ d 2 a j 2 d^2a_j^2 d 2 a j 2
min j ∣ λ j − ρ ∣ ≤ ∥ r ∥ 2 . \boxed{\min_j|\lambda_j-\rho|\le\|r\|_2.} j min ∣ λ j − ρ ∣ ≤ ∥ r ∥ 2 . 这保证附近至少有一个特征值，还没有指定是哪一个。
若要判断与一个简单特征值 λ ℓ \lambda_\ell λ ℓ
g = min j ≠ ℓ ∣ λ j − ρ ∣ > 0. g=\min_{j\ne\ell}|\lambda_j-\rho|>0. g = j = ℓ min ∣ λ j − ρ ∣ > 0. 只保留残差展开中 j ≠ ℓ j\ne\ell j = ℓ
∥ r ∥ 2 2 ≥ g 2 ∑ j ≠ ℓ a j 2 = g 2 sin 2 θ , sin θ ≤ ∥ r ∥ 2 g . \|r\|_2^2\ge
g^2\sum_{j\ne\ell}a_j^2
=g^2\sin^2\theta,
\qquad
\boxed{\sin\theta\le\frac{\|r\|_2}{g}.} ∥ r ∥ 2 2 ≥ g 2 j = ℓ ∑ a j 2 = g 2 sin 2 θ , sin θ ≤ g ∥ r ∥ 2 . 这里的间隔是其他特征值到当前估计 ρ \rho ρ g g g
取
A = diag ( 1 , 1 + δ ) , u = ( 1 , 1 ) T 2 , δ > 0. A=\operatorname{diag}(1,1+\delta),\qquad
u=\frac{(1,1)^T}{\sqrt2},\quad \delta>0. A = diag ( 1 , 1 + δ ) , u = 2 ( 1 , 1 ) T , δ > 0. 有 ρ = 1 + δ / 2 \rho=1+\delta/2 ρ = 1 + δ /2 ∥ r ∥ 2 = δ / 2 \|r\|_2=\delta/2 ∥ r ∥ 2 = δ /2 δ \delta δ u u u 45 ∘ 45^\circ 4 5 ∘
两个坐标方向标示特征方向，青色箭头表示固定的 u。特征值更接近后，残差变小，但夹角仍为 45 度；图中箭头长度不用于比较向量的模。
非对称问题不能照搬这个界 看看
A M = [ 1 M 0 2 ] , u = ( M , 2 ) T M 2 + 4 , M > 0. A_M=\begin{bmatrix}1&M\\0&2\end{bmatrix},\qquad
u=\frac{(M,2)^T}{\sqrt{M^2+4}},\qquad M>0. A M = [ 1 0 M 2 ] , u = M 2 + 4 ( M , 2 ) T , M > 0. 它的特征值一直是 1 和 2，但
ρ ( u ) = 3 M 2 + 8 M 2 + 4 , ∥ A M u − ρ ( u ) u ∥ 2 = 2 M M 2 + 4 . \rho(u)=\frac{3M^2+8}{M^2+4},\qquad
\|A_Mu-\rho(u)u\|_2=\frac{2M}{M^2+4}. ρ ( u ) = M 2 + 4 3 M 2 + 8 , ∥ A M u − ρ ( u ) u ∥ 2 = M 2 + 4 2 M . 当 M M M
实际停止时，可以报告
η = ∥ A u − ρ u ∥ 2 ∥ A ∥ F + ∣ ρ ∣ , \eta=\frac{\|Au-\rho u\|_2}{\|A\|_F+|\rho|}, η = ∥ A ∥ F + ∣ ρ ∣ ∥ A u − ρ u ∥ 2 , 其中 ∥ A ∥ F \|A\|_F ∥ A ∥ F ρ \rho ρ
移位反迭代：把目标附近的方向放大 幂法只能优先留下最大模方向。如果目标在 10 附近，最大模特征值却是 1000，继续原来的乘法不会改变目标。
对移位 μ \mu μ A v j = λ j v j Av_j=\lambda_jv_j A v j = λ j v j
( A − μ I ) v j = ( λ j − μ ) v j . (A-\mu I)v_j=(\lambda_j-\mu)v_j. ( A − μ I ) v j = ( λ j − μ ) v j . 只要 μ \mu μ
( A − μ I ) − 1 v j = 1 λ j − μ v j . (A-\mu I)^{-1}v_j=\frac1{\lambda_j-\mu}v_j. ( A − μ I ) − 1 v j = λ j − μ 1 v j . 原来的特征方向没变，特征值经过“减去移位、再取倒数”变换。离 μ \mu μ
算法不形成逆矩阵。每轮解
( A − μ I ) y k = x k , x k + 1 = y k / ∥ y k ∥ ∞ , (A-\mu I)y_k=x_k,\qquad
x_{k+1}=y_k/\|y_k\|_\infty, ( A − μ I ) y k = x k , x k + 1 = y k /∥ y k ∥ ∞ , 再用原矩阵 A A A μ \mu μ A − μ I A-\mu I A − μ I
若 λ ∗ \lambda_\ast λ ∗ λ n e x t \lambda_{\mathrm{next}} λ next
∣ λ ∗ − μ λ n e x t − μ ∣ . \left|\frac{\lambda_\ast-\mu}
{\lambda_{\mathrm{next}}-\mu}\right|. λ next − μ λ ∗ − μ . 两个特征值与移位等距时，严格优势可能消失。移位精确等于特征值时，线性系统奇异，不能按普通求解步骤硬算。
同一组谱，换一个观察位置 令
A = diag ( 1 , 3 , 6 ) , μ = 14 5 , x 0 = ( 1 , 1 , 1 ) T . A=\operatorname{diag}(1,3,6),\qquad
\mu=\frac{14}{5},\qquad x_0=(1,1,1)^T. A = diag ( 1 , 3 , 6 ) , μ = 5 14 , x 0 = ( 1 , 1 , 1 ) T . 第一轮求解得到
y 0 = ( − 5 / 9 , 5 , 5 / 16 ) T , x 1 = ( − 1 / 9 , 1 , 1 / 16 ) T . y_0=(-5/9,\ 5,\ 5/16)^T,
\qquad
x_1=(-1/9,\ 1,\ 1/16)^T. y 0 = ( − 5/9 , 5 , 5/16 ) T , x 1 = ( − 1/9 , 1 , 1/16 ) T . 第二轮归一化后是
x 2 = ( 1 / 81 , 1 , 1 / 256 ) T . x_2=(1/81,\ 1,\ 1/256)^T. x 2 = ( 1/81 , 1 , 1/256 ) T . 目标方向是中间的坐标轴。与特征值 1 对应的分量每轮相对缩小 1 / 9 1/9 1/9 1 / 16 1/16 1/16
减去 14/5 再取倒数之后，中间特征值对应的放大倍数模最大。每一行保留同一个特征方向，改变的是该方向的伸缩倍数。
靠近目标通常加快方向收敛，同时也让线性系统接近奇异。这里不能只凭“病态”二字断言归一化后的方向一定很差：被放大的误差可能主要沿目标方向，归一化会消去一部分尺度影响。但必须用原矩阵重算残差，处理求解失败，并避免溢出。这个判断不能由迭代次数代替。
让移位跟着当前估计改变 把固定 μ \mu μ
在对称二阶问题中可以把原因算透。令 v 1 , v 2 v_1,v_2 v 1 , v 2 λ 1 ≠ λ 2 \lambda_1\ne\lambda_2 λ 1 = λ 2 v 1 + t v 2 v_1+tv_2 v 1 + t v 2 v 1 v_1 v 1
ρ = λ 1 + λ 2 t 2 1 + t 2 . \rho=\frac{\lambda_1+\lambda_2t^2}{1+t^2}. ρ = 1 + t 2 λ 1 + λ 2 t 2 . 用这个值移位并解一次系统后，两个方向的系数比变成
t n e w = t λ 1 − ρ λ 2 − ρ = − t 3 . t_{\mathrm{new}}
=t\frac{\lambda_1-\rho}{\lambda_2-\rho}
=-t^3. t new = t λ 2 − ρ λ 1 − ρ = − t 3 . 误差比例从 t t t − t 3 -t^3 − t 3 ∣ t ∣ = 1 |t|=1 ∣ t ∣ = 1 t = 0 t=0 t = 0
实验中先固定移位，查看同一份分解怎样被反复使用；然后改用当前 Rayleigh 商，比较分解次数与方向误差。把移位放在两个特征值中点，再精确放到某个特征值上，观察两种失败原因的区别。每次最终检查仍使用原来的 A A A
QR 迭代：把整组谱留在矩阵里 如果需要一组特征值，逐个试移位并不总是合适。第 7 章的 QR 分解提供了另一种组织方式。记 A 0 = A A_0=A A 0 = A
A k = Q k R k , A k + 1 = R k Q k , A_k=Q_kR_k,\qquad A_{k+1}=R_kQ_k, A k = Q k R k , A k + 1 = R k Q k , 其中 Q k Q_k Q k
Q k T A k Q k = Q k T Q k R k Q k = R k Q k , Q_k^TA_kQ_k
=Q_k^TQ_kR_kQ_k
=R_kQ_k, Q k T A k Q k = Q k T Q k R k Q k = R k Q k , 每一步都是正交相似变换。
相似为什么保留特征值？若 A k v = λ v A_kv=\lambda v A k v = λ v
A k + 1 ( Q k T v ) = Q k T A k v = λ ( Q k T v ) . A_{k+1}(Q_k^Tv)=Q_k^TA_kv
=\lambda(Q_k^Tv). A k + 1 ( Q k T v ) = Q k T A k v = λ ( Q k T v ) . Q k T v Q_k^Tv Q k T v
det ( λ I − A k + 1 ) = det ( Q k T ) det ( λ I − A k ) det ( Q k ) = det ( λ I − A k ) . \det(\lambda I-A_{k+1})
=\det(Q_k^T)\det(\lambda I-A_k)\det(Q_k)
=\det(\lambda I-A_k). det ( λ I − A k + 1 ) = det ( Q k T ) det ( λ I − A k ) det ( Q k ) = det ( λ I − A k ) . 只是坐标换了，特征向量的分量通常也跟着变。
对于主例，可以取
Q 0 = 1 5 [ 2 − 1 1 2 ] , R 0 = 1 5 [ 5 4 0 3 ] . Q_0=\frac1{\sqrt5}
\begin{bmatrix}2&-1\\1&2\end{bmatrix},
\qquad
R_0=\frac1{\sqrt5}
\begin{bmatrix}5&4\\0&3\end{bmatrix}. Q 0 = 5 1 [ 2 1 − 1 2 ] , R 0 = 5 1 [ 5 0 4 3 ] . 检查 Q 0 T Q 0 = I Q_0^TQ_0=I Q 0 T Q 0 = I Q 0 R 0 = A Q_0R_0=A Q 0 R 0 = A
A 1 = R 0 Q 0 = 1 5 [ 14 3 3 6 ] . A_1=R_0Q_0=\frac15
\begin{bmatrix}14&3\\3&6\end{bmatrix}. A 1 = R 0 Q 0 = 5 1 [ 14 3 3 6 ] . 非对角元从 1 变成 3 / 5 3/5 3/5 14 / 5 , 6 / 5 14/5,6/5 14/5 , 6/5
同一组 Q、R 按两种顺序相乘，得到正交相似的矩阵。右边的整体因子 1/5 作用于所有元素，因此非对角元实际为 3/5。
与幂法的联系，以及不能省掉的条件 累积正交变换
Z k = Q 0 Q 1 ⋯ Q k − 1 , Z 0 = I , Z_k=Q_0Q_1\cdots Q_{k-1},\qquad Z_0=I, Z k = Q 0 Q 1 ⋯ Q k − 1 , Z 0 = I , 便有 A k = Z k T A Z k A_k=Z_k^TAZ_k A k = Z k T A Z k A k = Q k R k A_k=Q_kR_k A k = Q k R k
A Z k = Z k + 1 R k . AZ_k=Z_{k+1}R_k. A Z k = Z k + 1 R k . 观察第一列，上三角矩阵的第一列只有首项，因此
A ( Z k e 1 ) = r 11 ( k ) Z k + 1 e 1 . A(Z_ke_1)=r_{11}^{(k)}Z_{k+1}e_1. A ( Z k e 1 ) = r 11 ( k ) Z k + 1 e 1 . 累积基的第一列实际上沿着幂法更新，只是把归一化信息放进了 R k R_k R k
主例是对称二阶矩阵，第一列向主特征方向靠近时，第二列因正交性也向另一个特征方向靠近。因此 Z k T A Z k Z_k^TAZ_k Z k T A Z k
比如
A = [ 0 1 1 0 ] A=\begin{bmatrix}0&1\\1&0\end{bmatrix} A = [ 0 1 1 0 ] 是正交矩阵。取 Q = A , R = I Q=A,R=I Q = A , R = I A A A 1 , − 1 1,-1 1 , − 1
移位必须加回来 带移位的一步写成
A k − μ k I = Q k R k , A k + 1 = R k Q k + μ k I . A_k-\mu_kI=Q_kR_k,\qquad
A_{k+1}=R_kQ_k+\mu_kI. A k − μ k I = Q k R k , A k + 1 = R k Q k + μ k I . 这样仍有 A k + 1 = Q k T A k Q k A_{k+1}=Q_k^TA_kQ_k A k + 1 = Q k T A k Q k
一种容易尝试的选择是末对角元 μ k = ( A k ) n n \mu_k=(A_k)_{nn} μ k = ( A k ) nn μ 0 = 2 \mu_0=2 μ 0 = 2 A 0 − 2 I A_0-2I A 0 − 2 I 2 I 2I 2 I
这里还有一个与反迭代不同的地方：移位矩阵奇异时，QR 分解本身仍然存在，只是 R R R R x = b R x=b R x = b
在实验里逐步查看 Q , R Q,R Q , R
结构、去耦与停止 对一个稠密大矩阵，每轮都完整做 QR 代价较高。常用准备步骤是通过左右配对的 Householder 反射，将它变成上 Hessenberg 矩阵：
H = Z T A Z , h i j = 0 ( i > j + 1 ) . H=Z^TAZ,\qquad h_{ij}=0\quad(i>j+1). H = Z T A Z , h ij = 0 ( i > j + 1 ) . 这里只消去第一条次对角线以下的元素。左边做一次反射改变行，右边乘上同一反射保证相似性；不能只做左变换。后续反射避开已经处理的坐标，便可逐列保留先前的零。
上 Hessenberg 矩阵的 QR 消元每次只需处理相邻两行，重新相乘也保持这种结构，所以一次结构化 QR 步只需 O ( n 2 ) O(n^2) O ( n 2 ) O ( n 3 ) O(n^3) O ( n 3 )
上 Hessenberg 结构只要求第一条次对角线以下为零；再加上对称性，第一条超对角线以外的上方元素也必须为零。叉号表示可非零，并不要求一定非零。
以对称二阶块说明“去耦”：
B = [ a b b d ] . B=\begin{bmatrix}a&b\\b&d\end{bmatrix}. B = [ a b b d ] . 当 b b b
Δ B = [ 0 − b − b 0 ] , ∥ Δ B ∥ 2 = ∣ b ∣ . \Delta B=\begin{bmatrix}0&-b\\-b&0\end{bmatrix},
\qquad \|\Delta B\|_2=|b|. Δ B = [ 0 − b − b 0 ] , ∥Δ B ∥ 2 = ∣ b ∣. 这个操作使两个坐标方向分开，每个对角元都可以作为一个特征值近似。若 B = Z T A Z B=Z^TAZ B = Z T A Z Z e 2 Ze_2 Z e 2
A ( Z e 2 ) − d ( Z e 2 ) = Z ( b , 0 ) T , A(Ze_2)-d(Ze_2)=Z(b,0)^T, A ( Z e 2 ) − d ( Z e 2 ) = Z ( b , 0 ) T , 残差长度恰好为 ∣ b ∣ |b| ∣ b ∣
对一般实非对称矩阵，目标常是准上三角形式：实特征值对应一阶块，共轭复特征值对应二阶块。实数算法不能保证把含复特征值的矩阵变成实对角矩阵。这里已经足够解释 QR 的基本选择；隐式双移位、复杂去耦准则和大型稀疏谱算法不在本章的实现范围内。
练习：同时检查计算和结论 选择方法时，先弄清需要哪部分谱，以及矩阵允许什么运算。只需最大模方向、又能快速算 A x Ax A x O ( n 3 ) O(n^3) O ( n 3 ) O ( n 2 ) O(n^2) O ( n 2 )
题 1 对 A = diag ( − 4 , 2 ) A=\operatorname{diag}(-4,2) A = diag ( − 4 , 2 ) x 0 = ( 1 , 1 ) T x_0=(1,1)^T x 0 = ( 1 , 1 ) T
查看解答 第一轮 x 1 = ( − 1 , 1 / 2 ) T x_1=(-1,1/2)^T x 1 = ( − 1 , 1/2 ) T x 2 = ( 1 , 1 / 4 ) T x_2=(1,1/4)^T x 2 = ( 1 , 1/4 ) T x k = ( ( − 1 ) k , 2 − k ) T x_k=((-1)^k,2^{-k})^T x k = (( − 1 ) k , 2 − k ) T e 1 e_1 e 1 − e 1 -e_1 − e 1 e 1 e_1 e 1 2 − k 2^{-k} 2 − k 1 / 2 1/2 1/2 − 4 -4 − 4
题 2 对 diag ( 5 , 2 , 1 ) \operatorname{diag}(5,2,1) diag ( 5 , 2 , 1 ) ( 0 , 1 , 1 ) T (0,1,1)^T ( 0 , 1 , 1 ) T diag ( 2 , − 2 ) \operatorname{diag}(2,-2) diag ( 2 , − 2 ) ( 1 , 1 ) T (1,1)^T ( 1 , 1 ) T
查看解答 第一个起点没有 e 1 e_1 e 1 e 2 e_2 e 2
第二个例子交替得到两个方向 ( 1 , 1 ) T (1,1)^T ( 1 , 1 ) T ( 1 , − 1 ) T (1,-1)^T ( 1 , − 1 ) T
题 3 令 A = diag ( 4 , 1 ) A=\operatorname{diag}(4,1) A = diag ( 4 , 1 ) u = ( 3 , 1 ) T / 10 u=(3,1)^T/\sqrt{10} u = ( 3 , 1 ) T / 10 ρ , r , ∥ r ∥ 2 \rho,r,\|r\|_2 ρ , r , ∥ r ∥ 2 e 1 e_1 e 1 u u u 10 − 8 u 10^{-8}u 1 0 − 8 u
查看解答 加权平均给出 ρ = 37 / 10 \rho=37/10 ρ = 37/10
r = 1 10 ( 9 / 10 , − 27 / 10 ) T , ∥ r ∥ 2 = 9 / 10. r=\frac1{\sqrt{10}}(9/10,-27/10)^T,\qquad \|r\|_2=9/10. r = 10 1 ( 9/10 , − 27/10 ) T , ∥ r ∥ 2 = 9/10. 至少有一个特征值距离 ρ \rho ρ 0.9 0.9 0.9 0.3 0.3 0.3 e 1 e_1 e 1 ρ \rho ρ g = 27 / 10 g=27/10 g = 27/10 sin θ ≤ 1 / 3 \sin\theta\le1/3 sin θ ≤ 1/3 sin θ = 1 / 10 \sin\theta=1/\sqrt{10} sin θ = 1/ 10
缩放向量后 Rayleigh 商和方向不变，未经归一化的残差长度缩小 10 8 10^8 1 0 8
题 4 把正文聚集谱例子的起点改为 u = ( 3 / 2 , 1 / 2 ) T u=(\sqrt3/2,1/2)^T u = ( 3 /2 , 1/2 ) T e 1 e_1 e 1 30 ∘ 30^\circ 3 0 ∘
查看解答 对 A = diag ( 1 , 1 + δ ) A=\operatorname{diag}(1,1+\delta) A = diag ( 1 , 1 + δ )
ρ = 1 + δ / 4 , r = ( − 3 δ / 8 , 3 δ / 8 ) T , ∥ r ∥ 2 = 3 δ / 4. \rho=1+\delta/4,\qquad
r=(-\sqrt3\,\delta/8,\ 3\delta/8)^T,\qquad
\|r\|_2=\sqrt3\,\delta/4. ρ = 1 + δ /4 , r = ( − 3 δ /8 , 3 δ /8 ) T , ∥ r ∥ 2 = 3 δ /4. δ → 0 \delta\to0 δ → 0 30 ∘ 30^\circ 3 0 ∘
题 5 对正文的 A M , u A_M,u A M , u M = 1000 M=1000 M = 1000
查看解答 代入得 ρ ≈ 2.999996 \rho\approx2.999996 ρ ≈ 2.999996 0.001999992 0.001999992 0.001999992 0.999996 0.999996 0.999996 E = − r u T E=-ru^T E = − r u T ρ \rho ρ A M + E A_M+E A M + E 0.002 0.002 0.002 ρ \rho ρ
题 6 令 A = diag ( − 2 , 1 , 5 ) A=\operatorname{diag}(-2,1,5) A = diag ( − 2 , 1 , 5 ) μ = 4 \mu=4 μ = 4 ( 1 , 1 , 1 ) T (1,1,1)^T ( 1 , 1 , 1 ) T
查看解答 线性系统的对角元为 − 6 , − 3 , 1 -6,-3,1 − 6 , − 3 , 1 ( − 1 / 6 , − 1 / 3 , 1 ) T (-1/6,-1/3,1)^T ( − 1/6 , − 1/3 , 1 ) T 1 / 3 1/3 1/3
移位 3 到 1 和 5 等距，这两个方向的放大模相同，通常不会选出唯一方向。移位 5 使系数矩阵奇异；若当前还没有通过特征残差检查，就应报告求解问题，而非继续除法。
题 7 在特征基中，A = diag ( 4 , 1 ) A=\operatorname{diag}(4,1) A = diag ( 4 , 1 ) ( 1 , 1 / 2 ) T (1,1/2)^T ( 1 , 1/2 ) T
查看解答 当前 ρ = 17 / 5 \rho=17/5 ρ = 17/5
diag ( 3 / 5 , − 12 / 5 ) y = ( 1 , 1 / 2 ) T \operatorname{diag}(3/5,-12/5)y=(1,1/2)^T diag ( 3/5 , − 12/5 ) y = ( 1 , 1/2 ) T 得到 y = ( 5 / 3 , − 5 / 24 ) T y=(5/3,-5/24)^T y = ( 5/3 , − 5/24 ) T ( 1 , − 1 / 8 ) T (1,-1/8)^T ( 1 , − 1/8 ) T t n e w = − ( 1 / 2 ) 3 t_{\mathrm{new}}=-(1/2)^3 t new = − ( 1/2 ) 3
4 + 1 / 64 1 + 1 / 64 = 257 65 , \frac{4+1/64}{1+1/64}=\frac{257}{65}, 1 + 1/64 4 + 1/64 = 65 257 , 离 4 还差 3 / 65 3/65 3/65
题 8 对 A = [ 3 1 1 1 ] A=\begin{bmatrix}3&1\\1&1\end{bmatrix} A = [ 3 1 1 1 ]
Q = 1 10 [ 3 − 1 1 3 ] , R = 1 10 [ 10 4 0 2 ] . Q=\frac1{\sqrt{10}}\begin{bmatrix}3&-1\\1&3\end{bmatrix},
\qquad
R=\frac1{\sqrt{10}}\begin{bmatrix}10&4\\0&2\end{bmatrix}. Q = 10 1 [ 3 1 − 1 3 ] , R = 10 1 [ 10 0 4 2 ] . 计算一次 QR 换序结果，并求原坐标中 Q e 2 Qe_2 Q e 2
查看解答 两列正交且单位长，乘回 Q R QR QR
R Q = 1 10 [ 34 2 2 6 ] . RQ=\frac1{10}\begin{bmatrix}34&2\\2&6\end{bmatrix}. R Q = 10 1 [ 34 2 2 6 ] . 迹为 4，行列式为 2，与原矩阵相同。取 d = 3 / 5 d=3/5 d = 3/5 u = Q e 2 = ( − 1 , 3 ) T / 10 u=Qe_2=(-1,3)^T/\sqrt{10} u = Q e 2 = ( − 1 , 3 ) T / 10
A u − d u = 1 10 ( 3 / 5 , 1 / 5 ) T , Au-du=\frac1{\sqrt{10}}(3/5,1/5)^T, A u − d u = 10 1 ( 3/5 , 1/5 ) T , 长度为 1 / 5 1/5 1/5
题 9 对 B = [ 2 b b 2 ] B=\begin{bmatrix}2&b\\b&2\end{bmatrix} B = [ 2 b b 2 ] b > 0 b>0 b > 0 b = 10 − 4 b=10^{-4} b = 1 0 − 4 10 − 6 10^{-6} 1 0 − 6
查看解答 移位后为 b [ 0 1 1 0 ] b\begin{bmatrix}0&1\\1&0\end{bmatrix} b [ 0 1 1 0 ] Q Q Q R = b I R=bI R = b I B B B
真实特征值为 2 + b , 2 − b 2+b,2-b 2 + b , 2 − b b = 10 − 4 b=10^{-4} b = 1 0 − 4 10 − 6 10^{-6} 1 0 − 6
10 一个单位向量的特征残差很小，就足以断定它接近最大模特征值对应的唯一方向。
