QR方法特征分解
- 1. 实对称矩阵的QR方法特征分解
- 2. 代码实现与验证
- 3. 参考资料
1. 实对称矩阵的QR方法特征分解
使用QR分解的方法求取实对称矩阵特征值、特征向量，首先通过householder变换将实对称矩阵转化为三对角矩阵，对三对角矩阵进行QR分解迭代
[ × × × × × × × × × × × × × × × × ] → H = I − 2 w ∗ w T [ × × 0 0 × × × 0 0 × × × 0 0 × × ] \begin{bmatrix}\times&\times&\times&\times\\\times&\times&\times&\times\\\times&\times&\times&\times\\\times&\times&\times&\times\end{bmatrix}\xrightarrow{H=I-2w*w^T}\begin{bmatrix}\times&\times&0&0\\\times&\times&\times&0\\0&\times&\times&\times\\0&0&\times&\times\end{bmatrix} ×××××××××××××××× H=I−2w∗wT ××00×××00×××00××
针对Hessenberg矩阵多采用带位移的QR算法1(如: Rayleigh quotient shift，Wilkinson shift)
{ A k − σ ∗ I = Q k ∗ R k A k + 1 = R k ∗ Q k + σ ∗ I \left\{ \begin{array}{l} {A}_{k} - \sigma * I = {Q}_{k} * {R}_{k} \\ {A}_{k + 1} = {R}_{k} * {Q}_{k} + \sigma * I \end{array}\right. {Ak−σ∗I=Qk∗RkAk+1=Rk∗Qk+σ∗I
记 Q ^ k = Q 0 ∗ Q 1 ∗ … ∗ Q k {\widehat{Q}}_{k} = {Q}_{0} * {Q}_{1} * \ldots * {Q}_{k} Q k=Q0∗Q1∗…∗Qk Rayleigh商加速取 σ = H n n \sigma = {H}_{nn} σ=Hnn ,(如下 H H H 是QR迭代中的某一步 A k {A}_{k} Ak )
H = [ × × × × × × × × 0 × × × 0 0 α λ n ] \begin{array}{r} H = \left\lbrack \begin{matrix} \times & \times & \times & \times \\ \times & \times & \times & \times \\ 0 & \times & \times & \times \\ 0 & 0 & \alpha & {\lambda }_{n} \end{matrix}\right\rbrack \end{array} H= ××00×××0×××α×××λn
由特征多项式知道,当 α = 0 \alpha = 0 α=0 时, λ n {\lambda }_{n} λn 是 H H H 的一个特征值,实际计算时,当 α < ϵ , λ n \alpha < \epsilon ,{\lambda }_{n} α<ϵ,λn 即所求特征值。进而目标矩阵变为 ( n − 1 ) \left( {n - 1}\right) (n−1) 阶的。
假设 σ = H n n = λ n \sigma = {H}_{nn} = {\lambda }_{n} σ=Hnn=λn 是(近似)特征向量,则针对(近似)奇异矩阵 ( A k − σ ∗ I ) \left( {{A}_{k} - \sigma * I}\right) (Ak−σ∗I) 的 QR分解得到的 Q , R Q,R Q,R 形状如下 R n n = ≈ 0 {R}_{nn} = \approx 0 Rnn=≈0 ;,进而新的 A k + 1 {A}_{k + 1} Ak+1
H n , n − 1 = R ( n , : ) ∗ Q ( : , n − 1 ) = ≈ 0 {H}_{n,n - 1} = R\left( {n, : }\right) * Q\left( { : ,n - 1}\right) = \approx 0 Hn,n−1=R(n,:)∗Q(:,n−1)=≈0
A n e w = R ∗ Q = [ × × × × 0 × × × 0 0 × × 0 0 0 ε ≈ 0 ] ∗ [ × × × × × × × × 0 × × × 0 0 × × ] = [ × × × × × × × × 0 × × × 0 0 α ≈ 0 λ n ] {A}_{new} = R * Q = \left\lbrack \begin{matrix} \times & \times & \times & \times \\ 0 & \times & \times & \times \\ 0 & 0 & \times & \times \\ 0 & 0 & 0 & \varepsilon \approx 0 \end{matrix}\right\rbrack * \left\lbrack \begin{matrix} \times & \times & \times & \times \\ \times & \times & \times & \times \\ 0 & \times & \times & \times \\ 0 & 0 & \times & \times \end{matrix}\right\rbrack = \left\lbrack \begin{matrix} \times & \times & \times & \times \\ \times & \times & \times & \times \\ 0 & \times & \times & \times \\ 0 & 0 & \alpha \approx 0 & {\lambda }_{n} \end{matrix}\right\rbrack Anew=R∗Q= ×000××00×××0×××ε≈0 ∗ ××00×××0×××××××× = ××00×××0×××α≈0×××λn
待求解矩阵阶数降低一个。逐次降低，直至2阶、1阶; 完毕。
算法伪代码2：
图1. 算法伪代码 2. 代码实现与验证
使用Matlab代码实现，其中涉及到householder变换到hessenberg矩阵(三对角矩阵)，针对hessenberg矩阵进行单位移迭代
function [hatQ,D,iter]=QR_eig(A) %实对称矩阵QR分解 n = size(A, 1); %化为Hessenberg矩阵 [hatQ,A]=rsmhessenb(A); k=n; iter=0; D = eye(n); %保存特征值 while k>1 iter=iter+1; sigma=A(k,k); [Q,R]=QR_givens(A-sigma*eye(k));%Givens变换实现QR分解 A=R*Q+sigma*eye(k);%移位 D(1:k,k+1:end) = Q.'*D(1:k,k+1:end); hatQ(:,1:k) = hatQ(:,1:k)*Q;%记录变换矩阵 if(abs(A(k,k-1))<1e-11)%误差 D(1:k,k) = A(1:k,k);%记录右下角特征值 k=k-1;%收缩矩阵 A=A(1:k,1:k); end end D(1,1) = A(1,1); end function [Q,A]=rsmhessenb(A) %实对称矩阵化为Hessenberg矩阵 n=length(A(:,1)); Q=eye(n); for i=2:n sigma=sign(A(i,i-1))*sqrt(sum(A(i:n,i-1).^2)); rou=sigma*(A(i,i-1)+sigma); A(i,i-1)=A(i,i-1)+sigma; R=eye(n-i+1)-1.0/rou*A(i:n,i-1)*A(i:n,i-1)'; A(i,i-1)=-sigma; A(i+1:n,i-1)=0; A(i-1,i:n)=A(i:n,i-1);%AT=A %A(1:i-1,i:n)=A(1:i-1,i:n)*R;%AT\=A A(i:n,i:n)=R*A(i:n,i:n)*R; m=length(R(:,1)); U=[eye(n-m,n-m),zeros(n-m,m);zeros(m,n-m),R]; Q=Q*U; end end function [Q, R] = QR_givens(A) %%输入的为三对角矩阵 [m, n] = size(A); % 获取矩阵 A 的尺寸 Q = eye(n); % 初始化 Q 为 m x m 的单位矩阵 R = A; % 初始化 R 为 A 的拷贝 for k = 1:n-1 i = k + 1; % 计算 Givens 旋转参数 c 和 s % [c, s] = givens(R(k, k), R(i, k)); c=R(k,k)/sqrt(R(k,k)^2+R(i,k)^2); s=-R(i,k)/sqrt(R(k,k)^2+R(i,k)^2); % 构造 Givens 旋转矩阵 G G = [c s; -s c]; T=Q(:, [k i]); % 应用 Givens 旋转到 R 的第 k 行和第 i 行 R([k i], k:n) = G' * R([k i], k:n); % 应用 Givens 旋转到 Q 的第 k 列和第 i 列 Q(:, [k i]) = Q(:, [k i]) * G; end end同样将其改为C语言实现
void QR_givens(double A[], double Q[], int n) { //对三对角矩阵进行Givens-QR分解 int i,j,k; for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { Q[i * n + j] = (i == j) ? 1.0 : 0.0; } } double c, s,tau,temp; for (k = 0; k < n - 1; k++) { i = k + 1; tau=sqrt(A[k*n+k]*A[k*n+k]+A[i*n+k]*A[i*n+k]); c=A[k*n+k]/tau; s=-A[i*n+k]/tau; //应用Givens旋转到R的第k行和第i行 for (j= k; j < n; j++) { temp=A[k*n+j]; A[k*n+j]= c* A[k*n+j] -s* A[i*n+j]; A[i*n+j]= s* temp +c* A[i*n+j]; } A[i*n+k]=0.0; // 应用Givens旋转到Q的第k列和第i列 for (j = 0; j < n; j++) { temp=Q[j*n+k]; Q[j*n+k] =c* Q[j*n+k]-s* Q[j*n+i]; Q[j*n+i] =s*temp+c* Q[j*n+i]; } } } // 计算特征值和特征向量的QR法 // A：输入n*n的矩阵，输出没有意义 // D：输入n*n的矩阵用于中间计算，输出第一行为特征值 // p：输入n*n的矩阵用于中间计算，输出没有意义 // Q: 输入n*n的矩阵为hessenberg变换矩阵，输出特征向量 // eps：精度 void QR_eig(double A[], double D[], double p[],double Q[], int n, double eps) { int i, j, m, k = n; double sigma; for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { p[i * n + j] = (i == j) ? 1.0 : 0.0; } } while (k > 1) { sigma = A[k*k-1];//选取A[K][K] for (j = 0; j < k; j++) { A[j * k + j] = A[j * k + j] - sigma; }//A-sigma*I QR_givens(A, p, k);//Givens变换实现QR分解，其中A为R,p为Q //记录R*Q+sigma*I for (i = 0; i < k; i++) { for (j = 0; j < k; j++) { D[i * k + j] = 0.0; for (m = 0; m < k; m++) { D[i * k + j] += A[i * k + m] * p[m * k + j]; } if (i == j) D[i * k + j] += sigma; } } //重新赋值给A_new，即A=R*Q+sigma*I for (i = 0; i < k; i++) { for (j = 0; j < k; j++) { A[i * k + j]=D[i*k+j]; } } //记录变换矩阵Q'=Q*p for (i = 0; i < n; i++) { // 对于hatQ的每一行 for (j = 0; j < k; j++) { // 对于前k列的每一个元素 // 清零当前元素 D[i*k+j] = 0.0; // 计算当前元素的值，它是第i行和第j列的Q矩阵的和乘以对应hatQ元素的总和 for (m = 0; m < k; m++) { D[i*k+j] += Q[i*n+m] * p[m*k+j]; } } } for (i = 0; i < n; i++) { for (j = 0; j < k; j++) { Q[i * n + j]=D[i*k+j]; } } //达到误差要求，更改矩阵 if (fabs(A[k*k-2]) < eps){ k-=1; for(i=0;i测试代码
int main() { int n = 4,i,j,k; // 矩阵的大小 double A[16] = { // 初始对称矩阵的元素 1, 3, 1, 4, 3, 2, 0, 1, 1, 0, 2, 3, 4, 1, 3, 2 }; double Q[16],D[16],p[16]; double u[4]; double t[4]={0.0}; // 调用函数 rsmhessenb(A, n, Q, u, t); QR_eig(A,D,p,Q,n,0.0000001); printf("eigenvalue:\n"); for (int i = 0; i < n; i++) { printf("%f ", D[i]); } printf("\neigen vector:\n"); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { printf("%f ", Q[i * n + j]); } printf("\n"); } return 0; }测试结果
图2. QR分解Matlab结果 图3. C代码结果 3. 参考资料
数值分析8(带位移的QR算法和一般矩阵的处理)苏州大学_哔哩哔哩_bilibili ↩︎
喻文健主编. 数值分析与算法 第3版[M]. 北京：清华大学出版社, 2020.03 ↩︎