本文已参与「新人创作礼」活动,一起开启掘金创作之路。
Power Method Matlab Code
知道什么是power method,为什么需要,懂原理,并且想直接用的同学直接看下面的代码即可。如果不懂原理得同学,可以先跳过这一节,看下一节的内容。
另外,我这个例子是我上《矩阵分析》课程的一个作业,请随意使用。
clear;clc;
% generate A:
m = 10;
A = cell(m,m);
for i=1:m
for j=1:m
if i==j
A{i,j}=i+j;
else
A{i,j}=i*j;
end
end
end
A = cell2mat(A);
eigen_pairs = cell(m,2);
convergence_error = 0.0001;
max_iterate_num = 1000;
% initialize random positive vector x:
x = ones(m,1);
lambda_dominant = 0;
A_pre = A;
v_dominant = x;
% calculate all eigenvalues and eigenvectors:
for i=1:m
A_cur = A_pre - lambda_dominant * (v_dominant * v_dominant');
disp(A_cur);
iterate_count = 0;
x_pre = x;
% iterate until convergence or reach the iterate max limit
while iterate_count < max_iterate_num
x_cur = A_cur * x_pre;
error = norm(abs(x_cur/x_cur(end)) - abs(x_pre/x_pre(end)), 1) ;
% convergence condition:
if error < convergence_error
break
end
x_pre = x_cur/norm(x_cur);
iterate_count = iterate_count+ 1;
end
v = x_pre;
% using Rayleigh quotient to calculate lambda
lambda = (v' * A_cur' * v) / (v'*v);
% store results:
eigen_pairs{i, 1} = v;
eigen_pairs{i, 2} = lambda;
lambda_dominant = lambda;
v_dominant = v;
A_pre = A_cur;
end
% print out results:
for i=1:m
fprintf('the %d th eigenvalue: %f \n', i, eigen_pairs{i, 2});
fprintf(' and eigenvector: \n');
disp(eigen_pairs{i, 1});
end
什么是power method
首先,对一矩阵,我们要求它的特征值和特征向量,怎么求?一般当然是直接调包了。。。不过你知道是包里面怎么实现的吗?另外,对于有些问题,掉包可能反而效率低,所以我们可以自己实现一个方法来求特征值和特征向量。一个经典的方法就是power method。不过呢,这个方法一般就是用来找最大的一组特征值和特征向量。
在某些问题中,我们只需要找到最大的主特征值及其对应的特征向量。在这种情况下,我们可以使用幂法——一种将收敛到最大特征值的迭代方法。下面让我们看看 power 方法是如何工作的。
考虑一个𝑛×𝑛n×n矩阵𝐴A具有𝑛n线性独立实特征值𝜆1,𝜆2,…,𝜆𝑛λ1,λ2,…,λn和相应的特征向量𝑣1,𝑣2,…,𝑣𝑛v1,v2,…,vn. 由于特征值是标量,我们可以对它们进行排序,以便|𝜆1|>|𝜆2|>⋯>|𝜆𝑛||λ1|>|λ2|>⋯>|λn|(实际上,我们只需要|𝜆1|>|𝜆2||λ1|>|λ2|,其他特征值可能彼此相等)。
可以直接参考