[解決済み] MATLABによるパワーメソッド

質問

を実装したいのですが パワーメソッド MATLABで行列の支配的な固有値と固有ベクトルを決定する。

ここまで書いたものがこちらです。

%function to implement power method to compute dominant
%eigenvalue/eigenevctor
function [m,y_final]=power_method(A,x);
m=0;
n=length(x);
y_final=zeros(n,1);
y_final=x;
tol=1e-3;
while(1)
    mold=m;
 y_final=A*y_final;
 m=max(y_final);
 y_final=y_final/m;
 if (m-mold)<tol
     break;
 end
end
end

上記のコードで、数値の例を示します。

 A=[1 1 -2;-1 2 1; 0 1 -1]

A =

     1     1    -2
    -1     2     1
     0     1    -1

>> x=[1 1 1];
>> x=x';
>> [m,y_final]=power_method(A,x);
>> A*x

ans =

     0
     2
     0

MATLABで上の行列の固有値、固有ベクトルと比較すると、そうなりました。

[V,D]=eig(A)

V =

    0.3015   -0.8018    0.7071
    0.9045   -0.5345    0.0000
    0.3015   -0.2673    0.7071


D =

    2.0000         0         0
         0    1.0000         0
         0         0   -1.0000

固有値は一致するが、固有ベクトルは接近するはずである [1/3 1 1/3] . ここで、私は得る。

 y_final

y_final =

    0.5000
    1.0000
    0.5000

この不正確さを見てもいいのでしょうか、それとも私が何か勘違いをしているのでしょうか？

解決方法は？

正しい実装をしているが、チェックが甘い 両方 固有ベクトルと固有値が収束しているかどうか。 固有値の収束だけをチェックしているのです。 べき乗法では、顕著な固有ベクトルと固有値の両方を推定するので、おそらく 両方とも は収束した。 そうすると、なんとか [1/3 1 1/3] . これを容易にするために、あなたのコードをどのように修正したかを説明します。

function [m,y_final]=power_method(A,x)
m=0;
n=length(x);
y_final=x;
tol=1e-10; %// Change - make tolerance more small to ensure convergence
while(1)
     mold = m;
     y_old=y_final; %// Change - Save old eigenvector
     y_final=A*y_final;
     m=max(y_final);
     y_final=y_final/m;
     if abs(m-mold) < tol && norm(y_final-y_old,2) < tol %// Change - Check for both
         break;
     end
end
end

上記のコードを例の入力で実行すると、こうなります。

>> [m,y_final]=power_method(A,x)

m =

     2


y_final =

    0.3333
    1.0000
    0.3333

について余談ですが eig MATLABはその固有ベクトルを別のノルムを用いてスケーリングしている可能性が高いです。 固有ベクトルは は一意ではない であり、スケールアップまで正確です。 もし確認したいのであれば、単に V で割ると、Power 法と同じように、1成分が1の値で正規化されます。

>> [V,D] = eig(A);
>> V(:,1) / max(abs(V(:,1)))

ans =

    0.3333
    1.0000
    0.3333

これは、あなたが観察したものと一致します。