QR方法求解特征值、特征向量(附Matlab、C代码)

QR方法特征分解

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​×××0​0×××​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​∗Rk​Ak+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<k;i++){
            for(j=0;j<k;j++){
                A[i*k+j]=A[i*(k+1)+j];
            }
        }
    }
  }
  //变换后的A矩阵保留特征值，读取特征值保存到D矩阵的第一行
  for(i=1;i<=n;i++){
    D[i-1]=A[i*i-1];
  }
}

测试代码

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. 参考资料

  1. 数值分析8(带位移的QR算法和一般矩阵的处理)苏州大学_哔哩哔哩_bilibili ↩︎

  2. 喻文健主编. 数值分析与算法 第3版[M]. 北京：清华大学出版社, 2020.03 ↩︎