Jacobi方法特征分解学习笔记

fwfwfwfwfwfwfwfwfw

已于 2024-05-25 12:49:44 修改

阅读量735

点赞数 26

文章标签：矩阵笔记

于 2024-05-25 12:48:51 首次发布

本文链接：https://blog.csdn.net/qq_55009067/article/details/139194986

版权

Jacobi方法求解特征值、特征向量

1. Givens变换

Jacobi 方法是用来计算实对称矩阵的全部特征值和特征向量的一种方法, 它的基本思想是: 因为任何一个实对称矩阵都与一个对角矩阵相似, 所以通过正交相似变换把实对称矩阵化为对角矩阵, 对角矩阵的对角线上的元素即为所求特征值. 由线性代数的知识得, 若 $\mathbf{A} {\in} {\mathbf{R}}^{n {\times} n}$ 为对称矩阵, 则存在一个正交矩阵 $\mathbf{P}$ , 使得
$\mathbf{PA}{\mathbf{P}}^{{-}1} = \operatorname{diag}\left( {\lambda}_{1},{\lambda}_{2},{\cdots},{\lambda}_{n} \right) = \mathbf{D}$

${\lambda}_{1},{\lambda}_{2},{\cdots},{\lambda}_{n}$ 即为矩阵 $\mathbf{A}$ 的特征值, ${\mathbf{P}}^{\mathrm{T}}$ 的 $n$ 个列向量 $v_{1},v_{2},{\cdots},v_{n}$ 即为所对应的特征向量
$\mathbf{P} = \begin{bmatrix} 1 & & & & & & & \\ & {\ddots} & & & & & & \\ & & \cos\theta & & & & \sin\theta & \\ & & & 1 & & & & \\ & & & & {\ddots} & & & \\ & & & & & 1 & & \\ & & {-} \sin\theta & & & & \cos\theta & \\ & & & & & & & 1 \\ & & & & & & & & {\ddots} \\ & & & & & & & & & 1 \end{bmatrix} {\equiv} \mathbf{P}(i, j)$
称 $\mathbf{P}$ 为平面旋转矩阵， $P A$ 只改变 $A$ 的第 $i$ 行与第 $j$ 行的元素; $AP^{\mathrm{T}}$ 只改变 $A$ 的第 $i$ 列与第 $j$ 列的元素; $PAP^{T}$ 只改变 $A$ 的第 $i$ 行和第 $j$ 行, 第 $i$ 列与第 $j$ 列的元素，考虑 $n = 2$ 的情况, 设 $\mathbf{A} = \left\lbrack \begin{array}{ll} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array} \right\rbrack$ 为对称矩阵, 今取 $\mathbf{P} = \begin{bmatrix} \cos\theta & \sin\theta \\ {-} \sin\theta & \cos\theta \end{bmatrix}$ , 显然 $\mathbf{P}$ 为正交矩阵, 令 $\mathbf{PA}{\mathbf{P}}^{\mathrm{T}} = \mathbf{C} = \left\lbrack \begin{array}{ll} c_{11} & c_{12} \\ c_{21} & c_{22} \end{array} \right\rbrack$ , 由矩阵乘法得

$c_{11} = a_{11}{\cos}^{2}\theta + a_{21}\sin 2\theta + a_{22}{\sin}^{2}\theta,$
$c_{22} = {-} a_{12}\sin 2\theta + a_{11}{\sin}^{2}\theta + a_{22}{\cos}^{2}\theta,$
$c_{21} = c_{12} = \frac{1}{2}\left( a_{22} {-} a_{11} \right)\sin 2\theta + a_{21}\cos 2\theta.$

为使得 $\mathbf{P}\mathbf{A}{\mathbf{P}}^{\mathrm{T}} = \mathbf{C}$ 成为对角矩阵, 应选择 $\theta$ 使 $c_{21} = c_{12} = 0$ , 即 $\frac{1}{2}\left( a_{22} {-} a_{11} \right)\sin 2\theta + a_{21}\cos 2\theta = 0,$ 得 $\tan 2\theta = \frac{2a_{21}}{a_{11} {-} a_{22}}$ , 由此可得 $\theta$ 的值. 若 $a_{11} = a_{22}$ 时, 取 $|\theta| = \frac{\pi}{4};a_{11} > 0$ 时, $\theta = \frac{\pi}{4};a_{11} < 0$ 时, $\theta = {-} \frac{\pi}{4}$ , 结果就使得 $\mathbf{P}\mathbf{A}{\mathbf{P}}^{\mathrm{T}} = \operatorname{diag}\left( {\lambda}_{1},{\lambda}_{2} \right)$ ，则 $\mathbf{C} = \mathbf{P}\mathbf{A}{\mathbf{P}}^{\mathrm{T}}$ 的元素 $c_{ij}$ 的计算公式为:

$c_{ii} = a_{ii}{\cos}^{2}\theta + a_{ij}\sin 2\theta + a_{jj}{\sin}^{2}\theta$ , $c_{jj} = {-} a_{ij}\sin 2\theta + a_{ii}{\sin}^{2}\theta + a_{jj}{\cos}^{2}\theta.$
$c_{ij} = c_{ji} = \frac{1}{2}\left( a_{jj} {-} a_{ii} \right)\sin 2\theta + a_{ij}\cos 2\theta$ .
第 $i$ 行元素 $c_{ik} = c_{ki} = a_{ik}\cos\theta + a_{jk}\sin\theta,k {\neq} i,j$ .
第 $j$ 行元素 $c_{jk} = c_{kj} = a_{jk}\cos\theta {-} a_{ik}\sin\theta,k {\neq} i,j$ .
第 $i$ 列元素 $c_{ki} = a_{ki}\cos\theta + a_{kj}\sin\theta,k {\neq} i,j$ .
第 $j$ 列元素 $c_{kj} = a_{kj}\cos\theta {-} a_{ki}\sin\theta,k {\neq} i,j$ .
其他元素不变, $c_{lk} = a_{lk},l,k {\neq} i,j$ .

由此可见, 若矩阵 $\mathbf{A}$ 的非对角元素 $a_{ij} {\neq} 0$ , 我们就可以选择一个正交矩阵 $\mathbf{P}(i,j)$ 使得 $\mathbf{C} = \mathbf{P}\mathbf{A}{\mathbf{P}}^{\mathrm{T}}$ 的元素 $c_{ij} = c_{ji} = 0$ , 即选择 $\theta$ 满足

$\tan 2\theta = \frac{2a_{ij}}{a_{ii} {-} a_{jj}},|\theta| {\leq} \frac{\pi}{4}.$

通过不断左乘右乘旋转矩阵，就能完成对角化，对角矩阵即为特征值，旋转矩阵的乘积为特征向量。

2. 经典Jacobi

算法步骤：

在 $\mathbf{A}$ 的非对角元素中选取一个绝对值最大的元素 (称为主元素), 设 $\left| {\mathbf{a}}_{i_{1},j_{1}} \right| = \mathop{\max}\limits_{l {\neq} k}$ $\left| a_{lk} \right|$ , 可设 $a_{i_{1},j_{1}} >tol$ 即大于精度误差, 否则, 认为 $A$ 已对角化迭代结束;
构造旋转矩阵 ${\mathbf{P}}_{1}\left( i_{1},j_{1} \right)$ 使得 ${\mathbf{A}}_{1} = {\mathbf{P}}_{1}\mathbf{A}{\mathbf{P}}_{1}^{{\top}}$ 的非对角元 $a_{i_{1}j_{1}}^{(1)} = a_{j_{1}i_{1}}^{(1)} = 0$
更新元素 $c_{ij}$ 进行下一轮迭代

Matlab代码如下：

function [V,D,iter]=Jocobi_classical(A,maxIter,tol)
n = size(A, 1); % 矩阵的大小
V = eye(n); % 初始化特征向量矩阵为单位矩阵
iter = 0; % 初始化迭代次数
% 设置最大迭代次数和误差精度
if nargin < 3 || isempty(tol)
    tol = 1e-9; % 默认误差精度
end
if nargin < 2 || isempty(maxIter)
    maxIter = 1000; % 默认最大迭代次数
end
while(iter < maxIter)
    iter=iter+1;
    D=A;
    n=size(D,1);
    p=1;q=2;
    for i=1:n
        for j=i+1:n
            if(abs(D(i,j))>abs(D(p,q)))%找到对称矩阵的上三角矩阵中最大的元素的下标
                p=i;q=j;
            end
        end
    end
    if(abs(D(p,q))<tol)
        break;
    end
    if(A(p,q)~=0)
        d=(A(q,q)-A(p,p))/(2*A(p,q));
        if(d>0)
            t=1/(d+sqrt(d^2+1));
        else
            t=-1/(-d+sqrt(d^2+1));
        end
        c=1/sqrt(t^2+1);s=c*t;
    else
        c=1;s=0;
    end
    R=[c s;-s c];
    A([p,q],:)=R'*A([p,q],:);
    A(:,[p,q])=A(:,[p,q])*R;
    V(:, [p, q]) = V(:,[p,q])*R;
end

在这里插入图片描述

图1. 测试矩阵

在这里插入图片描述

图2. eig函数结果

在这里插入图片描述

图3. 经典Jacobi结果 最后将该算法移植为C语言实现

// 计算特征值和特征向量的雅可比迭代法
// a：输入n*n的矩阵，计算后对角元素为特征值
// n：矩阵a的阶数
// v：用于存储特征向量的数组
// eps：精度
// jt：最大迭代次数
// 返回值：1表示成功，-1表示迭代次数超过最大迭代次数
int jcbi(double a[], int n, double v[], double eps, int jt) {
  int i, j, p, q, u, w, t, s, l;
  double fm, cn, sn, omega, x, y, d;
  l = 1;
  // 初始化v为单位矩阵
  for (i = 0; i <= n - 1; i++) {
    v[i * n + i] = 1.0;
    for (j = 0; j <= n - 1; j++)
      if (i != j)
        v[i * n + j] = 0.0;
  }
  while (1) {
    fm = 0.0;
    // 下三角找最大值
    for (i = 1; i <= n - 1; i++)
      for (j = 0; j <= i - 1; j++) {
        d = fabs(a[i * n + j]);
        if ((i != j) && (d > fm)) {
          fm = d;
          p = i;
          q = j;
        }
      }
    if (fm < eps)
      return (1);
    if (l > jt)
      return (-1);
    l = l + 1;     // 迭代次数
    u = p * n + q; // a[p][q]
    w = p * n + p; // a[p][p]
    t = q * n + p; // a[q][p]
    s = q * n + q; // a[q][q]
    x = -a[u];
    y = (a[s] - a[w]) / 2.0;
    omega = x / sqrt(x * x + y * y);
    if (y < 0.0)
      omega = -omega;
    sn = 1.0 + sqrt(1.0 - omega * omega);
    sn = omega / sqrt(2.0 * sn);
    cn = sqrt(1.0 - sn * sn);
    fm = a[w];
    // 更新a[p][p]，a[q][q]，a[p][q]，a[q][p]
    a[w] = fm * cn * cn + a[s] * sn * sn + a[u] * omega;
    a[s] = fm * sn * sn + a[s] * cn * cn - a[u] * omega;
    a[u] = 0.0;
    a[t] = 0.0;
    // 更新a[p][*],a[q][*]
    for (j = 0; j <= n - 1; j++)
      if ((j != p) && (j != q)) {
        u = p * n + j;
        w = q * n + j;
        fm = a[u];
        a[u] = fm * cn + a[w] * sn;
        a[w] = -fm * sn + a[w] * cn;
      }
    // 更新a[*][p],a[*][q]
    for (i = 0; i <= n - 1; i++)
      if ((i != p) && (i != q)) {
        u = i * n + p;
        w = i * n + q;
        fm = a[u];
        a[u] = fm * cn + a[w] * sn;
        a[w] = -fm * sn + a[w] * cn;
      }
    // 更新v[*][*]
    for (i = 0; i <= n - 1; i++) {
      u = i * n + p;
      w = i * n + q;
      fm = v[u];
      v[u] = fm * cn + v[w] * sn;
      v[w] = -fm * sn + v[w] * cn;
    }
  }
  return (1);
}