Homepage of Jinghua Xu

2025-04-10发表2025-04-29更新课程项目 / 工程数值分析1 小时读完 (大约6923个字)

研究背景意义

矩阵的特征值计算虽然有比较可靠的理论方法，但是，理论方法只适合于矩阵规模很小或者只是在理论证明中起作用，而实际问题的数据规模都比较大，不太可能采用常规的理论解法。计算机擅长处理大量的数值计算，所以通过适当的数值计算理论，写成程序，让计算机处理，是一种处理大规模矩阵的方法，而且是一种好的方法。乘幂法（又称幂法）是求矩阵按模最大特征值的常用方法，其结构简单、便于使用，在实际工程中应用广泛：

在数值分析和优化问题中，乘幂法可通过求解矩阵按模最大特征值来得到其谱半径，从而帮助分析迭代算法的收敛性；
在物理学和工程学中，特征值问题常用于分析系统的稳定性、振动频率和模态分析，乘幂法能够有效求解系统的最大特征值，从而帮助理解系统的动态行为；
在图像处理中，特征值分解用于图像压缩、特征提取和模式识别，乘幂法可用于求解图像矩阵的最大特征值，从而实现高效的图像处理算法；
……

但该方法也存在一定的局限性，其只适用于求矩阵按模最大的特征值，而在实际问题中往往需要求矩阵全部特征值，如在主成分分析中需要求相关矩阵的全部特征值以确定各主成分的贡献率，在求解线性微分方程组中通过求系数矩阵的全部特征值以确定其基础解系；同时该方法在求特定矩阵（如具有一对互为相反数的按模最大特征值的矩阵等）的按模最大特征值时是不收敛的。因此，我们希望以乘幂法基本思想为出发点，在此基础上提出一些改进算法以对其进行泛化，使其能够解决更为广泛且通用的应用场景下的实际问题。

算法基本思想

设n阶矩阵A具有n个线性无关的特征向量
$$
x_1,x_2,…,x_n
$$
相应的特征值
$$
\lambda_1,\lambda_2,…,\lambda_n
$$
满足：
$$
\left|\lambda_1\right|>\left|\lambda_2\right|\geq\left|\lambda_3\right|\geq…\geq\left|\lambda_n\right|
$$
现任取一非零向量u_0，作迭代
$$
u_k=Au_{k-1}, k=0,1,2,…
$$
得到向量序列{u_k}(k=0,1,2,…)。因各特征向量线性无关，故n维向量u_0必可由他们线性表示，即：
$$
u_0=\alpha_{1}x_{1}+\alpha_{n}x_{2}+…+\alpha_{n}x_{n}
$$
显然有：
$$
\begin{aligned}u_{k}&=A^{k}v_{0}=\alpha_{1}A^{k}x_{1}+\alpha_{2}A^{k}x_{2}+\cdots+\alpha_{n}A^{k}x_{n}=\alpha_{1}\lambda_{1}^{k}x_{1}+\alpha_{2}\lambda_{2}^{k}x_{2}+\cdots+\alpha_{n}\lambda_{n}^{k}x_{n}\&=\lambda_{1}^{k}\left[\alpha_{1}x_{1}+\alpha_{2}\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{k}x_{2}+\cdots+\alpha_{n}\left(\frac{\lambda_{n}}{\lambda_{1}}\right)^{k}x_{n}\right]\end{aligned}
$$
设
$$
\alpha_{1} \neq 0
$$
由
$$
\left|\frac{\lambda_{i}}{\lambda_{1}}\right|<1,i=2,3,\cdots,n
$$
可得：
$$
\operatorname*{lim}{k\to\infty}\frac{u{k}}{\lambda_{1}^{k}}=\operatorname*{lim}{k\to\infty}\frac{\lambda{1}^{k}\left[\alpha_{1}x_{1}+\sum_{i=2}^{n}\alpha_{i}\left(\frac{\lambda_{i}}{\lambda_{1}}\right)^{k}x_{i}\right]}{\lambda_{1}^{k}}=\operatorname*{lim}{k\to\infty}\left[\alpha{1}x_{1}+\sum_{i=2}^{n}\alpha_{i}\left(\frac{\lambda_{i}}{\lambda_{1}}\right)^{k}x_{i}\right]=\alpha_{1}x_{1}
$$

$$
\operatorname*{lim}{k\to\infty}\frac{\left(u{k+1}\right){m}}{\left(u{k}\right){m}}=\operatorname*{lim}{k\to\infty}\frac{\left{\lambda_{1}^{k+1}\left[\alpha_{1}x_{1}+\sum_{i=2}^{n}\alpha_{i}\left(\frac{\lambda_{i}}{\lambda_{1}}\right)^{k+1}x_{i}\right]\right}{m}}{\left{\lambda{1}^{k}\left[\alpha_{1}x_{1}+\sum_{i=2}^{n}\alpha_{i}\left(\frac{\lambda_{i}}{\lambda_{i}}\right)^{k}x_{i}\right]\right}{m}}=\lambda{1}\frac{\left(x_{1}\right){m}}{\left(x{1}\right){m}}=\lambda{1}
$$

这表明向量序列u_k具有收敛性，其收敛速度由比值
$$
\left|\frac{\lambda_{2}}{\lambda_{1}}\right|
$$
确定，该比值越小说明收敛速度越快，比值越接近于1收敛速度就越慢。

同时，该结论表明，当k取得足够大时（在实际应用时往往认为某次迭代前后误差小于设定的误差限即满足该条件），有：
$$
u_{k}≈\lambda_{1}^{k}\alpha_{1}x_{1}, \lambda_{1}≈\frac{\left(u_{k+1}\right)}{\left(u_{k}\right)}
$$
即u_ k可近似看作特征值λ_ 1对应的特征向量（较原特征向量x_ 1而言仅乘上了有限的常数系数，其结果仍为该特征值对应的特征向量），而对应的按模最大特征值λ_ 1的估计值可由前后两次迭代的特征向量u_ k相除得到。

当矩阵的按模最大特征值是重根时，上述结论仍然成立。设λ_ 1为r重根，即：
$$
\lambda_{1}=\lambda_{2}=\cdots=\lambda_{r}
$$
且满足条件
$$
\mid\lambda_{r}\mid>\mid\lambda_{r+1}\mid\geq\cdots\geq\mid\lambda_{n}\mid
$$
则有：
$$
u_{k}=A^{k}u_{0}=\lambda_{1}^{k}\left[\sum_{i=1}^{r}\alpha_{i}x_{i}+\sum_{j=r+1}^{n}\alpha_{j}\left(\frac{\lambda_{j}}{\lambda_{1}}\right)^{k}x_{j}\right]
$$
易推得（式中(u_ k)_ m表示向量u_ k的第m个分量）：
$$
\lim_{k\to\infty}\frac{u_{k}}{\lambda_{i}^{k}}=\sum_{i=1}^{^{r}}\alpha_{i}x_{i},\lim_{k\to\infty}\frac{\left(u_{k+1}\right){m}}{\left(u{k}\right){m}}=\lambda{1}
$$
与一般情况略有不同的是，此时
$$
u_{k}≈\lambda_{i}^{k}\sum_{i=1}^{^{r}}\alpha_{i}x_{i}
$$
但事实上我们并没有办法在求解之前事先预知一个矩阵是否具有r重按模最大特征值，因此在实际操作的过程中采用与一般情形相同的处理方式（即认为r=1），所得的向量序列u_k仍然收敛，不影响求解结果。

基于以上的算法原理，可以编写如下MATLAB程序，实现一般情况下矩阵按模最大特征值及其对应特征向量的计算：

function calculatePowerMethod(matrixEdit, iterEdit)
    % 获取用户输入的矩阵matrixEdit和最大迭代次数iterEdit
    A = str2num(matrixEdit.Value); % 将字符串转换为矩阵
    n = length(A);
    V = rand(n,1);  % 以随机方式初始化迭代向量u0
    max_iter = iterEdit.Value;    % 获取迭代次数

    Eps = 1E-4; % 迭代精度
    k = 0; % 初始迭代次数
    lambda0 = 0;

    % 乘幂法迭代，求得矩阵按模最大特征值lambda及其对应特征向量v
    while k <= max_iter - 1
        v = A * V;
        [vmax, i] = max(abs(v));
        lambda = v(i) / V(i);
        V = v;
        if abs(lambda - lambda0)<Eps
            k = k + 1;
            break;
        end
        lambda0 = lambda;
        k = k + 1;
    end
end

改进的乘幂法

在实际应用层面，上述的算法逻辑只在一小部分的常规矩阵上对于矩阵的按模最大特征值求解有较好效果。事实上，仅依托以上的代码逻辑进行求解往往会出现以下几大问题：

当矩阵A不具有n个线性无关的特征向量时（事先无法判断），乘幂法不适用。
当待求取的按模最大特征值abs(λ_ 1)>1时，迭代向量u_ k的各个分量可能会随着abs(λ_ 1)^k变得很大而使计算机“上溢”；而当待求取的按模最大特征值abs(λ_ 1)<1时，迭代向量u_ k的各个分量可能会随着abs(λ_ 1)^k变得很小使u_k成为零向量。
矩阵的按模最大特征值是一对相反数。
矩阵的按模最大特征值是一对共轭复数。

（上述算法中还给出了α_1≠0的假设，事实上如果不满足这一点也不影响乘幂法的成功使用。因为舍入误差的影响，在迭代某一步会产生u_k在x_1方向上的分量不为零，以后的迭代仍会收敛）

针对问题1，我们暂时无法提出具有针对性的解决方案，但后续的实际测试证明，当迭代次数k足够大时，我们针对问题3所提出的改进方案在该种情况下相比基本算法能够有更好的收敛性（收敛更快）；针对问题2，我们可以通过将向量序列u_k进行规范化处理，限制向量的模在特定小范围内，防止计算机运算的上下溢出；针对问题3，我们在原有的乘幂法算法思想的基础上进行了改进，改进后的算法能够有效对原算法无法求解（不收敛）的按模最大特征值互为相反数的矩阵的按模最大特征值进行有效求解；针对问题4，我们目前还没有找到有效的解决方案。

迭代向量的归一化

设迭代向量u为非零向量，将其归一化得到向量
$$
y=\frac{u}{max(u)}
$$
其中max(u)表示向量u的模最大的分量（即向量的无穷范数）。这样的规范化会保证每一次迭代后得到的新的迭代向量的模最大分量均为1，这有效起到了限制迭代向量分量范围的作用，能够很好地解决上述提到的问题2。值得注意的是，对向量规范化的方式并不只有归一化这一种，实际上可以使用任意一种范数对向量进行规范化处理，如2范数等。可以证明，对于基于不同范数的归一化都有相同的向量序列收敛结论，下面以无穷范数为例给出基本的证明。（显然这样的规范化只是在原有的向量基础上除以了某一个常数，这并不会改变向量序列的收敛性）

取初始向量
$$
u_{0}\neq0
$$
规范化得
$$
y_{0}=\frac{u_{0}}{\max(u_{0})}
$$
构造向量序列：
$$
\begin{aligned}
&u_{1}=Ay_{0}=\frac{Au_{0}}{\max(u_{0})},
&y_{1}=\frac{u_{1}}{\max(u_{1})}=\frac{Au_{0}}{\max(Au_{0})}\
&u_{2}=Ay_{1}=\frac{A^{2}u_{0}}{\max(Au_{0})},
&y_{2}=\frac{u_{2}}{\max(u_{2})}=\frac{A^{2}u_{0}}{\max(A^{2}u_{0})}\
&……\
&u_{k}=Ay{{k-1}}=\frac{A^{k}u{0}}{\max(A^{k-1}u_{0})},
&\quad y_{k}=\frac{u_{k}}{\max(u_{k})}=\frac{A^{k}u_{0}}{\max(A^{k}u_{0})}
\end{aligned}
$$
结合乘幂法基本算法中已经证明的结论
$$
A^{k}u_{0}=u_{k}≈\lambda_{1}^{k}\alpha_{1}x_{1}
$$
则有：
$$
\lim_{k\to\infty}y_{k}=\frac{x_{1}}{\max(x_{1})},
\lim_{k\to\infty}\max(u_{k})=\lambda_{1}
$$
这表明此时y_ k可近似看作特征值λ_ 1对应的特征向量（较原特征向量x_ 1而言仅进行了归一化处理，其结果仍为该特征值对应的特征向量），而对应的按模最大特征值λ_ 1的估计值可近似看作未归一化前的迭代向量最大分量值max(u_k)。

基于以上的算法原理，可以编写如下MATLAB程序，实现归一化处理下矩阵按模最大特征值及其对应特征向量的计算：

function calculatePowerMethod(matrixEdit, iterEdit)
    % 获取用户输入的矩阵matrixEdit和最大迭代次数iterEdit
    A = str2num(matrixEdit.Value); % 将字符串转换为矩阵
    n = length(A);
    V = rand(n,1);  % 以随机方式初始化迭代向量u0
    max_iter = iterEdit.Value;    % 获取迭代次数

    Eps = 1E-4; % 迭代精度
    k = 0; % 初始迭代次数
    lambda0 = 0;

    % 乘幂法迭代，求得矩阵按模最大特征值lambda及其对应特征向量v
    while k <= max_iter - 1
        v = A * V;
        [vmax, i] = max(abs(v));
        lambda = v(i);
        V = v / lambda;
        if abs(lambda - lambda0)<Eps
            k = k + 1;
            break;
        end
        lambda0 = lambda;
        k = k + 1;
    end
end

归一化乘幂法的进一步改进

针对问题3，从其基本情形出发（其他条件与先前保持一致），即有：
$$
\mid\lambda_1\mid=\mid\lambda_2\mid>\mid\lambda_3\mid\geqslant\cdotp\cdotp\cdotp\geqslant\mid\lambda_n\mid，\lambda_1=-\lambda_2
$$
在矩阵按模最大特征值为一对相反数的情况下，原迭代向量序列中各相邻向量的比值极限发散：
$$
\lim_{k\to\infty}\frac{\boldsymbol{u}_{k+1}}{\boldsymbol{u}k}=\lambda_1\lim{k\to\infty}\frac{a_1\boldsymbol{x}_1+(-1)^{k+1}a_2\boldsymbol{x}2+\boldsymbol{\varepsilon}{k+2}}{a_1\boldsymbol{x}_1+(-1)^ka_2\boldsymbol{x}_2+\boldsymbol{\varepsilon}_k}
$$
说明在该种情况下，采用原有的乘幂法基本算法，得到的迭代向量序列并不收敛，无法得到稳定的矩阵按模最大特征值数值求解结果。

我们在原有算法的基础上进行了一定的改进，构造了新的迭代向量序列，使其对于矩阵按模最大特征值为一对相反数的情况也能保证收敛，从而可以对该情况下的矩阵按模最大特征值进行有效求解。