応用数理実験

固有値問題 (1) 固有値問題に対する冪乗法
桂田 祐史

2005 年 6 月 28 日
今回は固有値問題を取り扱うための冪乗法とその周辺の種々の手法を実験してもらう。
実験に用いる行列は小さなもので良い。自分で選ぶのが面倒なら、次の行列を使ってよい1 。


1 1 1
A= 1 2 2 .
1 2 3• 選んだ行列のサイズが小さい場合は、固有多項式を計算して固有値に関する情報を調べよ2 。サイズが小さくない場合は、Mathematica か (後で例をあげる)、MATLAB (Octave)の関数 eig() を用いて固有値を調べよ。
• 冪乗法により絶対値最大の固有値とそれに属する固有ベクトルを求めよ。収束の速さを調べよ。固有値を反復ベクトルの成分の比として計算した場合と、 Rayleigh 商で計算した場合とで、収束の速さに違いが出るか？
• 逆反復法により (連立 1 次方程式を解くのには LU 分解を使え)、絶対値最小の固有値とそれに属する固有ベクトルを求めよ。
• シフト法の実験をせよ。すなわち絶対値最大、最小の固有値以外の固有値を見い出せ。
(例えば、上に掲げた行列 A に対しては 0.6 に近い固有値を探してみよ。)• (余裕がある人向け) 絶対値最大の固有値に属する固有ベクトルを求めてから、その成分を含まない初期ベクトルを用いて冪乗法の反復を実行することにより、絶対値が二番目に大きな固有値の固有ベクトルを求めるアルゴリズムがある。適当な文献で調べるか、自分で考えて実験してみよ。
• (余裕がある人向け) 絶対値最大の固有値が二つ以上ある場合、冪乗法で何が起こるか観察せよ。
1

こ の 行 列 自 体 が 一 松 信「 数 値 解 析 」朝 倉 書 店 (1982) か ら 採った も の で あ る 。余 談 だ が 、この 本 は 基 本 的 な 数 値 計 算 ア ル ゴ リ ズ ム に 関 し て 概 観 す る に は 便 利 な 本 で あ る 。固 有 値 は λ1
=5.048917339522305313522214407023370089866 · · ·, λ2 = 0.6431041321077905561056004899786994619201 · · ·,λ3 = 0.3079785283699041303721851029979308815479 · · ·.
2
必要があれば Mathematica 等を利用せよ。その場合のヒント: 行列式を計算する Det[], 指定された次数の単位行列を作る IdentityMatrix[], 固有値を計算する Eigenvalues[], 方程式を解く Solve[], 方程式の解の近似値を求める FindRoot[], 関数のグラフを描く Plot[] 等の手続きがある。

1

¶Mathmatica の使用例

³oyabun% math
Mathematica 4.0 for SolarisCopyright 1988-1999 Wolfram Research, Inc.
-- Motif graphics initialized -In[1]:= A={{1,1,1},{1,2,2},{1,2,3}}
Out[1]= {{1, 1, 1}, {1, 2, 2}, {1, 2, 3}}
In[2]:= Eigenvalues[A]Out[2]= (省略 --- 興味があればやってみて下さい)
In[3]:= N[%,40]-46
Out[3]= {5.048917339522305313522214407023369723596 + 0. 10
I,-42>

0.3079785283699041303721851029979308598027 + 0. 10>

0.6431041321077905561056004899786994166009 + 0. 10

I,-42I}

In[4]:= na=N[A,50]
Out[4]= {{1.0000000000000000000000000000000000000000000000000,>
1.0000000000000000000000000000000000000000000000000,>
1.0000000000000000000000000000000000000000000000000},>{1.0000000000000000000000000000000000000000000000000,>
2.0000000000000000000000000000000000000000000000000,>
2.0000000000000000000000000000000000000000000000000},>{1.0000000000000000000000000000000000000000000000000,>
2.0000000000000000000000000000000000000000000000000,>
3.0000000000000000000000000000000000000000000000000}}
In[5]:= Eigenvalues[na]Out[5]= {5.0489173395223053135222144070233697235963877860565,>
0.64310413210779055610560048997869941660087281326538,>
0.30797852836990413037218510299793085980273940067811}
In[6]:= Eigenvectors[na]Out[6]= {{-0.32798527760568176779603202500454990640870532112133,>-0.59100904850610352545794571799937980844713881955014,>-0.73697622909957824233808630700517009796156650157119},>{0.73697622909957824233808630700517009796156650157119,>
0.32798527760568176779603202500454990640870532112133,>-0.59100904850610352545794571799937980844713881955014},>{-0.59100904850610352545794571799937980844713881955014,>
0.73697622909957824233808630700517009796156650157119,>-0.32798527760568176779603202500454990640870532112133}}
In[7]:= NSolve[Det[x IdentityMatrix[3]-A]==0, x, WorkingPrecision->30]}
Out[7]= {{x -> 0.30797852836990413037218510300},{x -> 0.64310413210779055610560048998},{x -> 5.0489173395223053135222144070}}
µ2

´

A

Octave のサンプル・プログラム

グラフ化の説明もかねて Octave のサンプル・プログラムを以下に示す。なお、MATLABでも同じプログラムが走ることを確認済み。
¶
冪乗法

³

format long
% 表示桁数を多く取る
a=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)% eig() で固有値を求める (カンニング)
x=ones(3,1);% 初期ベクトルを (1,1,1) とする
maxiter=10;
r=zeros(maxiter,1);% 近似固有値を記録するベクトルの用意
for k=1:maxiter
y=a*x;
x=y/norm(y);r(k)=x’*(a*x);% Rayleigh 商で近似固有値を計算し、r() に記録
enderror=r-max(eigen)% 近似固有値の誤差を求める
plot(1:maxiter,abs(error))% 反復回数と誤差をプロット
loglog(1:maxiter,abs(error))% 両側対数目盛りでプロット
semilogy(1:maxiter,abs(error)) % 片側対数目盛りでプロットµB´

MATLAB メモ• 6701 号室の Windows XP マシンでは、M-file は Z:Y か、あるいはその下、例えば
Z:Ymatlab に作ることを勧める。[File] メニューの [Set Path] サブメニューで [AddFolder] しておく。
• 2005 年度の情報処理教室のパソコンには、MATLAB や Octave がインストールされている。前者はライセンスの問題でつねに利用できるとは限らないが…

3

C

自分でやってみた

MATLAB を使って実験してみた。非常に気楽にプログラムを書いて出来てしまう。そろそろ課題を新しくして、この古い課題のプログラムをサンプルとして提示する形にするのが良いかもしれない。

C.1

冪乗法

¶³% eigen1.m --- power method
format longa=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)x=ones(3,1);
maxiter=10;
r=zeros(maxiter,1);
for k=1:maxiter
y=a*x;
x=y/norm(y);r(k)=x’*(a*x);endrerror=r-max(eigen)clf
semilogy(1:maxiter,abs(error))
µ´

10

10

10

10

10

10

10

10

-2-4-6-8-10-12-14-16

12

3

4

5

6

4

78

910

C.2

Rayleigh 商の効果¶³% naive_vs_rayleigh.m
format longa=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)x=ones(3,1);
maxiter=10;
r=zeros(maxiter,1);n=zeros(maxiter,1);
for k=1:maxiter
y=a*x;
x=y/norm(y);ax=a*x;r(k)=x’*ax;n(k)=ax(1)/x(1);endrnerror_r=r-max(eigen)
error_n=n-max(eigen)%subplot(1,2,1)%semilogy(1:maxiter,abs(error_r))%subplot(1,2,2)%semilogy(1:maxiter,abs(error_n))
clf
semilogy(1:maxiter,abs(error_r),1:maxiter,abs(error_n))µ´

1010101010101010100-2-4-6-8-10-12-14-16

123456578910

C.3

逆反復法

¶³% eigen2.m --- inverse iteration method
format longa=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)x=ones(3,1);
maxiter=30;
r=zeros(maxiter,1);[L u p]=lu(a);shift=1;
for k=1:maxiter
y=u\(L\(p*x));x=y/norm(y);r(k)=x’*(a*x);
endr
error=r-min(eigen)clf
semilogy(1:maxiter,abs(error))µ´

1010101010101010

100-2-4-6-8-10-12-14-160510156202530C.4

シフト法¶³% eigen3.m --- shift method
format longa=[1 1 1;1 2 2;1 2 3]
eigen=eig(a)x=ones(3,1);
maxiter=10;
r=zeros(maxiter,1);
shift=0.6;% we cannot use the factorization of a, can we?
[L u p]=lu(a-shift*eye(3,3));
for k=1:maxiter
y=u\(L\(p*x));x=y/norm(y);r(k)=x’*(a*x);
endr
eigen=sort(eigen);error=r-eigen(2)clfsemilogy(1:maxiter,abs(error))µ´

101010101010

1010-2-4-6-8-10-12-14-16

123456778910