【c++/c】幂法求解绝对值最大的特征值、最大和最小的特征值

【part 1：幂法的迭代格式】

　　　　

 

 【part 2：计算绝对值最大特征值的步骤】

　　1.明确要求解特征值的矩阵A和初始非零向量u0

　　2.将u0单位化得到y_k-1

　　3.通过矩阵乘法计算得到uk

　　4.通过矩阵乘法计算得到β_k（上图中的β_k-1应为β_k，为书写时的笔误），判断误差

　　　　如果误差小于允许误差，即为计算得到的绝对值最大的特征值

　　　　反之，则继续迭代计算

　　5.如果迭代次数大于设定的最大次数，计算失败。

示例：

step1：确定向量/矩阵

　　1.需要求解特征值的矩阵A：

#include <iostream>
#include <vector>
#include <math.h>

using vec = std::vector<std::vector<double>>;
using vecRow = std::vector<double>;

int g_time{ 5000 };
double g_err{ 10e-12 };

vec createA(double b = 0.16, double c = -0.064,int size = 501) {
    vec A(size, vecRow(size));
    double e{ exp(1) };
    for (int i = 0;i != size;++i) {
        for (int j = 0;j != size;++j) {
            if (i == j) {
                A[i][j] = (1.64 - 0.024 * (i+1)) * sin(0.2 * (i+1)) - 0.64 * exp(0.1/(i+1));
            }
            else if (fabs(i - j) == 1) {
                A[i][j] = b;
            }
            else if (fabs(i - j) == 2) {
                A[i][j] = c;
            }
            else {
                A[i][j] = 0;
            }
        }
    }
    return A;
}

　　2.初始非零向量u0。当然可以是别的形式。

vecRow createU(int num = 1,int size=501) {
    vecRow u(size);
    for (int i = 0;i < size;++i) {
        u[i] = num;
    }
    return u;
}

Step2  矩阵单位化

vecRow normalize(vecRow u) {
    double mid{};
    mid = multi(u, u);
    double eta{ sqrt(mid) };
    int row{ static_cast<int> (u.size()) };
    vecRow y(row);
    for (int i = 0;i < row;++i) {
        y[i] = u[i] / eta;
    }
    return y;
}

Step3 矩阵乘法

vec multi(vec A, vec B) {
    int row1{ static_cast<int> (A.size()) };
    int col1{ static_cast<int>(A[0].size()) };
    int row2{ static_cast<int>(B.size()) };
    int col2{ static_cast<int>(B[0].size()) };
    vec res(row1, vecRow(col2));
    double temp{};
    assert(col1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    for (int i = 0;i != row1;++i) {
        for (int j = 0;j != col2;++j) {
            temp = 0;
            for (int k = 0;k != col1;++k) {
                temp += A[i][k] * B[k][j];
            }
            res[i][j] = temp;
        }
    }
    return res;
}

double multi(vecRow A, vecRow B) {
    int row1{ static_cast<int> (A.size()) };
    int row2{ static_cast<int>(B.size()) };
    assert(row1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    double sum{0.0};
    for (int i = 0;i != row1;++i) {
        sum += A[i] * B[i];
    }
    return sum;
}

vecRow multi(vec A, vecRow y) {
    int row1{ static_cast<int> (A.size()) };
    int col1{ static_cast<int>(A[0].size()) };
    int row2{ static_cast<int>(y.size()) };
    vecRow result(row2);
    double sum{};
    assert(col1 == row2 && "第一个矩阵的列和第二个矩阵的行不相等");
    for (int i = 0;i != row1;++i) {
        sum = 0;
        for (int j = 0;j != col1;++j) {
            sum += A[i][j] * y[j];
        }
        result[i] = sum;
    }
    return result;
}

Step4 迭代计算

double mi(vec A, vecRow u0) {
    int k{ 1 };
    int row{ static_cast<int> (u0.size()) };
    vecRow y(row);
    vecRow uk(row);
    //double betaSet[static_cast<int>(g_time)]{};
    //vecRow betaSet(g_time);
   // betaSet[0]= 0;
    double beta0{};
    double betaStep{};
    double distance{};
    double beta{};
    while (k) {
        y = normalize(u0);
       // print(y);
        uk = multi(A, y);
        //betaSet[k] = multi(y, uk);
        //std::cout << betaSet[k];
        beta0 = betaStep;
        betaStep = multi(y, uk);
        distance = fabs((betaStep - beta0) / betaStep);
        if(distance<=g_err){
            beta = betaStep;
            //std::cout <<beta<< " is found\n";
            return beta;
        }
        if (k > g_time) {
            std::cout << "out of g_time" << '\n';
            return 0;
        }
        u0 = uk;
        //std::cout << "step " << k << '\n';
        ++k;
    }
    return 0;
}

 【part 3：计算最小和最大特征值】

 设绝对值最大的特征值为β

　　（1）beta > 0，则最大特征值为β

　　（2）beta < 0，则最小特征值为β

如何求解另一个最小/最大的特征值？ 令矩阵B为A-βI，利用幂法计算得出B的绝对值最大的特征值。特征值再加上减去的β，即为我们想要的另一个特征值。

step 1 创建B矩阵

vec createB(vec A, double beta) {
    int row{ static_cast<int> (A.size()) };
    int col{ static_cast<int>(A[0].size()) };
    for (int i = 0;i < row;++i) {
        for (int j = 0;j < col;++j) {
            if (i == j) {
                A[i][j] -= beta * 1;
            }
            else {
                continue;
            }
        }
    }
    return A;
}

Step 2 判断矩阵A的绝对值最大的特征值的正负

vecRow chargeMi(vec A,double beta,vecRow u0) {
    vec B{ createB(A, beta) };
    double lambda{ mi(B, u0) + beta };
    vecRow lambdaMinMax(2);
    if (beta > 0) {
        lambdaMinMax[0] = lambda;
        lambdaMinMax[1] = beta;
    }
    else {
        lambdaMinMax[0] = beta;
        lambdaMinMax[1] = lambda;
    }
    return lambdaMinMax;
}