如何求对称矩阵的特征值和特征向量

jocabi方法解决对称矩阵特征值问题

1.特征向量相关公式

对于给定方阵$A$ 存在非零向量$x$和标量$\lambda$

满足 $Ax = \lambda x$,则称 x 为 A 的特征向量，$\lambda$为$A$ 的特征值

从空间角度表示在该特定线性变换下特征向量方向的输入保持方向不变，只进行的伸缩

对于得到的n组特征值我们可以表示为：

$\Sigma_{n \times n} = \begin{pmatrix} \lambda_1 & & & \\ &\lambda_2 & & \\ & & ... & \\ & & & \lambda_n\end{pmatrix}$

n组特征向量(列向量)表示为

$X_{n \times n} = \begin{pmatrix} v_1, v_2, ..., v_n\end{pmatrix}$

则$X\Sigma = \begin{pmatrix}\lambda_1v_1, \lambda_2 v_2, ..., \lambda_n v_n \end{pmatrix}$

得到公式$AX = X\Sigma$ ，若$X$为可逆矩阵则$X^{-1}AX = \Sigma$

经过适当的变换就可以得到特征值构成的对角矩阵$\Sigma$

2. Jocabi方法迭代

对于对称矩阵$A$ ，提取部分元素$A' = \begin{pmatrix} a_{ii}&a_{ij} \\ a_{ij}&a_{jj}\end{pmatrix}$ 旋转矩阵$X' = \begin{pmatrix} cos\theta & -sin\theta \\ sin \theta & cos\theta\end{pmatrix}$

可得到 $\Sigma' = \begin{pmatrix} cos\theta & sin\theta \\ -sin \theta & cos\theta\end{pmatrix} \begin{pmatrix} a_{ii}&a_{ij} \\ a_{ij}&a_{jj}\end{pmatrix} \begin{pmatrix} cos\theta & -sin\theta \\ sin \theta & cos\theta\end{pmatrix}$

接下来公式密集，可能有算错的地方

$\sigma_{ii} = a_{ii}cos^{2}\theta+a_{jj}sin^{2}\theta + a_{ij}sin2\theta$

$\sigma_{jj} = a_{ii}sin^{2}\theta+a_{jj}cos^{2}\theta - a_{ij}sin2\theta$

$\sigma_{ij} = \sigma_{ji} = \frac{1}{2}sin2\theta(a_{jj} - a_{ii}) + a_{ij}cos2\theta$

实际上旋转矩阵是 $\begin{pmatrix}1&&&&&\\&...&&&&&\\&&cos\theta &&sin\theta \\&&&...&&& \\&&-sin\theta&&cos\theta \\&&&&&...&\\&&&&&&1 \end{pmatrix}$

则对于$\Sigma$ 整体其他行列还有

$\sigma_{ki}=\sigma_{ik} = a_{ik}cos \theta + a_{jk}sin\theta$

$\sigma_{kj}=\sigma_{jk} = -a_{ik}sin\theta + a_{jk}cos\theta$

上述两式$(k!=i,j)$

要消去 $\sigma_{ij},\sigma_{ji}$ 这类非对角线的元素 就需要巧妙的设置$\theta$（注意分子为0的情况讨论，如果使用cmath 的atan2(y,x)大致不需要）

$2\theta = arctan(\frac{2a_{ij}}{ a_{ii} - a_{jj}})$

$\theta = \frac{1}{2}arctan(\frac{2a_{ij}}{ a_{ii} - a_{jj}})$

然后将$A$更新为$\Sigma$即赋值 $A =\Sigma$通过不断迭代就可以消去除对角线外的元素

总体上说如果将$N$轮迭代合并得到式子$X^{-1}AX = \Sigma$

可得到特征向量$X = X_1X_2...X_n$ ，计算中我们可以初始化$X_0$ 为单位矩阵

得到$X = X_0X_1X_2...X_n$ 便于编程

可得到特征向量$X = X_1X_2...X_n$的证明如下

比如$X_{1}^{-1}A_1X_1 = A_2$

$X_{2}^{-1}A_2X_2 = A_3$

得到$X_{2}^{-1}X_{1}^{-1}A_{1}X_{1}X_{2}= A_3$ 由此不难推得

$X_n^{-1}X_2^{-1}...X_1^{-1}A_1X_1X_2...X_n = A_{N+1}$ 合并得到式子$X^{-1}AX = \Sigma$

$X$的更新方式

$x_{ki}= x_{ki}cos \theta + x_{kj}sin\theta$

$x_{kj}= -x_{ki}sin \theta + x_{kj}cos\theta$

3. 代码设计

主体思路

    while (epoch < MaxEpoch) {
        //循环次不超过最大次数的时候
        MaxElem(row, col, E);//找到最大元素的行和列，E为其绝对值
        if (E < Eps)break; //误差在一定范围内跳出
        Givens(row, col);  //利用上述旋转法旋转
        epoch++;
    }

完全代码如下

#pragma once
#include<vector>
#include<cmath>
​
/*
    jacobi方法求对称矩阵的特征值和特征向量
    Dim:        输入的维度
    MaxEpoch:   如果达不到指定误差范围，最多多少次退出
    Eps:        设定误差，如果达到误差则退出
    Mat:        输入的方阵，经过分解后对角线上是特征值
    Vec:        特征向量
*/
​
class Jacobi {
private:
    int Dim;
    int MaxEpoch;
    double Eps;
    std::vector<std::vector<double>> Mat;
    std::vector<std::vector<double>> Vec;
private:
    void MaxElem(int&, int&, double&);
    void Givens(int, int);
public:
    Jacobi() = default;
    Jacobi(std::vector <std::vector<double>>&, double, int);
    ~Jacobi() = default;
    std::vector<std::vector<double>> EigenVec();
    std::vector<std::vector<double>> EigenVal();
};
​
Jacobi::Jacobi(std::vector<std::vector<double>>& mat, double eps, int MaxN) {
    if (mat.size() != mat[0].size())throw("Dim");
    Eps = eps;
    MaxEpoch = MaxN;
    Dim = int(mat.size());
    Mat = mat;
    for (int i = 0; i < Dim; i++) {
        for (int j = 0; j < Dim; j++) {
            //确保矩阵是对称矩阵
            if (i != j && mat[i][j] != mat[j][i])throw("symmtry");
        }
    }
    Vec = std::vector<std::vector<double>>
        (Dim, std::vector<double>(Dim, 0));
    for (int i = 0; i < Dim; i++)Vec[i][i] = 1;
​
    int epoch = 0;
    int row = 0;
    int col = 0;
    double E = 0;
    while (epoch < MaxEpoch) {
        MaxElem(row, col, E);
        if (E < Eps)break;
        Givens(row, col);
        epoch++;
    }
​
}
​
void Jacobi::MaxElem(int& row, int& col, double& E) {
    //非对角线最大的元素
    float MaxE = fabs(Mat[0][1]);
    row = 0; col = 1;
    for (int i = 0; i < Dim; i++) {
        for (int j = 0; j < Dim; j++) {
            float e = fabs(Mat[i][j]);
            if (i != j && e > MaxE) {
                MaxE = e;
                row = i;
                col = j;
            }
        }
    }
    E = MaxE;
}
void Jacobi::Givens(int p, int q) {
    double y = 2 * Mat[p][q];
    double x = Mat[p][p] - Mat[q][q];
    double theta = 0.5 * atan2(y, x);
    double sint = sin(theta);
    double cost = cos(theta);
    double sin2t = sin(2.0 * theta);
    double cos2t = cos(2.0 * theta);
​
    double pp = Mat[p][p];
    double qq = Mat[q][q];
    double pq = Mat[p][q];
​
    Mat[p][p] = pp * (cost * cost) + qq * (sint * sint) + pq * sin2t;
    Mat[q][q] = pp * (sint * sint) + qq * (cost * cost) - pq * sin2t;
    Mat[q][p] = Mat[p][q] = 0.5 * (qq - pp) * sin2t + pq * cos2t;
​
    for (int i = 0; i < Dim; i++) {
        if (i != p && i != q) {
            double pi = Mat[p][i];
            double qi = Mat[q][i];
​
            Mat[p][i] = cost * pi + sint * qi;
            Mat[q][i] = - sint * pi + cost * qi;
            Mat[i][p] = Mat[p][i];
            Mat[i][q] = Mat[q][i];
        }
    }
    for (int k = 0; k < Dim; k++) {
        double Ekp = Vec[k][p];
        double Ekq = Vec[k][q];
​
        Vec[k][p] = Ekp * cost + Ekq * sint;
        Vec[k][q] = -Ekp * sint + Ekq * cost;
    }
}
​
std::vector<std::vector<double>> 
Jacobi::EigenVec() {
    return Vec;
}
​
std::vector<std::vector<double>> 
Jacobi::EigenVal(){
    return Mat;
}