行列指数関数（２）

Matlab/Scilabにはexpm()という行列指数関数が準備されている。Matlab/Scilabのexpm()では、スケーリングと正方化(scaling and squaring)されたパデ近似が使われている。

行列を対角化可能な場合は固有値と固有ベクトルに分解するmatrix decomposition methodがもっとも効率的。matrix decomposition methodはどこかで見たことがあると思ったら、中村先生のコンピュテーショナル・ファイナンスの2章の『至高のgenerator作り』でやったやつだな。あのときは1年後の格付け推移行列を求めた。

フリーのMatlabクローンであるScilabのサンプル・コードは以下のとおり。

*** sample code ******************************

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

Asave=A

expm(A) //Scilabの元々の行列指数関数

A=Asave

//スケーリングと正方化　scaling and squaring

[f,e]=log2(norm(A,'inf'))

s=max(0,e+1)

A=A/2^s

X=A

c=1/2

E=eye(3,3)+c*A

D=eye(3,3)-c*A

q=6

p=1

for k=2:q

c=c*(q-k+1)/(k*(2*q-k+1))

X=A*X

cX=c*X

E=E+cX

if p

D=D+cX

else

D=D-cX

end

p=~p

end

E=D\E

for k=1:s

E=E*E

end

E1=E

//固有値と固有ベクトルによる行列の指数 matrix decomposition methods

//例１

B=[0 1 2; 0.5 0 1; 2 1 0]

[D,V]=bdiag(B) //Matlabのeig()の場合、DとVが逆の可能性があるので後でチェック

E=diag(exp(diag(D)))

V*E*inv(V)

//例２

B=[-49 24;-64 31]

[D,V]=bdiag(B)

E=diag(exp(diag(D)))

V*E*inv(V)

*** sample code *********************************

このScilabでの実行結果は以下のとおり

スケーリングと正方化つきパデ近似の結果と行列分解法の結果が同じであることが分かる。

*** 出力 ****************************************

Scilabによる行列指数関数

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

A =

0. 1. 2.

0.5 0. 1.

2. 1. 0.

-->Asave=A

Asave =

0. 1. 2.

0.5 0. 1.

2. 1. 0.

-->expm(A) //Scilabの元々の行列指数関数

ans =

5.3090813 4.001203 5.5778403

2.8087901 2.8845155 3.1930144

5.173746 4.001203 5.7131756

-->A=Asave

A =

0. 1. 2.

0.5 0. 1.

2. 1. 0.

-->

-->//スケーリングと正方化　scaling and squaring

-->

-->[f,e]=log2(norm(A,'inf'))

e =

2.

f =

0.75

-->s=max(0,e+1)

s =

3.

-->A=A/2^s

A =

0. 0.125 0.25

0.0625 0. 0.125

0.25 0.125 0.

-->X=A

X =

0. 0.125 0.25

0.0625 0. 0.125

0.25 0.125 0.

-->c=1/2

c =

0.5

-->E=eye(3,3)+c*A

E =

1. 0.0625 0.125

0.03125 1. 0.0625

0.125 0.0625 1.

-->D=eye(3,3)-c*A

D =

1. - 0.0625 - 0.125

- 0.03125 1. - 0.0625

- 0.125 - 0.0625 1.

-->q=6

q =

6.

-->p=1

p =

1.

-->for k=2:q

--> c=c*(q-k+1)/(k*(2*q-k+1))

--> X=A*X

--> cX=c*X

--> E=E+cX

--> if p

--> D=D+cX

--> else

--> D=D-cX

--> end

--> p=~p

-->end

c =

0.1136364

X =

0.0703125 0.03125 0.015625

0.03125 0.0234375 0.015625

0.0078125 0.03125 0.078125

cX =

0.0079901 0.0035511 0.0017756

0.0035511 0.0026634 0.0017756

0.0008878 0.0035511 0.0088778

E =

1.0079901 0.0660511 0.1267756

0.0348011 1.0026634 0.0642756

0.1258878 0.0660511 1.0088778

D =

1.0079901 - 0.0589489 - 0.1232244

- 0.0276989 1.0026634 - 0.0607244

- 0.1241122 - 0.0589489 1.0088778

p =

F

c =

0.0151515

X =

0.0058594 0.0107422 0.0214844

0.0053711 0.0058594 0.0107422

0.0214844 0.0107422 0.0058594

cX =

0.0000888 0.0001628 0.0003255

0.0000814 0.0000888 0.0001628

0.0003255 0.0001628 0.0000888

E =

1.0080788 0.0662139 0.1271011

0.0348825 1.0027521 0.0644383

0.1262133 0.0662139 1.0089666

D =

1.0079013 - 0.0591116 - 0.1235500

- 0.0277802 1.0025746 - 0.0608872

- 0.1244377 - 0.0591116 1.0087891

p =

T

c =

0.0012626

X =

0.0060425 0.0034180 0.0028076

0.0030518 0.0020142 0.0020752

0.0021362 0.0034180 0.0067139

cX =

0.0000076 0.0000043 0.0000035

0.0000039 0.0000025 0.0000026

0.0000027 0.0000043 0.0000085

E =

1.0080865 0.0662182 0.1271046

0.0348864 1.0027547 0.0644409

0.126216 0.0662182 1.0089751

D =

1.0079089 - 0.0591073 - 0.1235464

- 0.0277764 1.0025771 - 0.0608846

- 0.1244350 - 0.0591073 1.0087975

p =

F

c =

0.0000631

X =

0.0009155 0.0011063 0.0019379

0.0006447 0.0006409 0.0010147

0.0018921 0.0011063 0.0009613

cX =

5.780D-08 6.984D-08 0.0000001

4.070D-08 4.046D-08 6.406D-08

0.0000001 6.984D-08 6.069D-08

E =

1.0080865 0.0662183 0.1271048

0.0348864 1.0027547 0.0644410

0.1262161 0.0662183 1.0089752

D =

1.0079089 - 0.0591074 - 0.1235465

- 0.0277764 1.0025771 - 0.0608846

- 0.1244352 - 0.0591074 1.0087975

p =

T

c =

0.0000015

X =

0.0005536 0.0003567 0.0003672

0.0002937 0.0002074 0.0002413

0.0003095 0.0003567 0.0006113

cX =

1.0D-09 *

0.8321428 0.5361264 0.5518949

0.4415159 0.3117848 0.3626738

0.4651685 0.5361264 0.9188691

E =

1.0080865 0.0662183 0.1271048

0.0348864 1.0027547 0.0644410

0.1262161 0.0662183 1.0089752

D =

1.0079089 - 0.0591074 - 0.1235465

- 0.0277764 1.0025771 - 0.0608846

- 0.1244352 - 0.0591074 1.0087975

p =

F

-->

-->E=D\E

E =

1.036393 0.1425675 0.2615269

0.0791531 1.0127849 0.1346981

0.2575922 0.1425675 1.0403277

-->

-->for k=1:s

--> E=E*E

-->end

E =

1.1527624 0.3294314 0.5623219

0.1968960 1.0562215 0.2972511

0.5462318 0.3294314 1.1688526

E =

1.7008831 0.9129553 1.4034189

0.5973082 1.2783914 0.7721247

1.3330036 0.9129553 1.7712983

E =

5.3090813 4.001203 5.5778403

2.8087901 2.8845155 3.1930144

5.173746 4.001203 5.7131756

-->E1=E

E1 =

5.3090813 4.001203 5.5778403

2.8087901 2.8845155 3.1930144

5.173746 4.001203 5.7131756

-->

-->//固有値と固有ベクトルによる行列の指数 matrix decomposition methods

-->//例１

-->B=[0 1 2; 0.5 0 1; 2 1 0]

B =

0. 1. 2.

0.5 0. 1.

2. 1. 0.

-->[D,V]=bdiag(B) //Matlabのeig()場合、DとVが逆の可能性があるので後で注意

V =

- 0.6540388 - 0.6360390 0.3508493

- 0.3800874 - 0.2120130 - 0.9055908

- 0.6540388 0.7420455 0.3508493

D =

2.5811388 0. 0.

0. - 2. 0.

0. 0. - 0.5811388

-->E=diag(exp(diag(D))) //要素ごとのexponentialを計算

E =

13.212176 0. 0.

0. 0.1353353 0.

0. 0. 0.5592611

-->V*E*inv(V)

ans =

5.3090813 4.001203 5.5778403

2.8087901 2.8845155 3.1930144

5.173746 4.001203 5.7131756

-->

-->//例２

-->B=[-49 24;-64 31]

B =

- 49. 24.

- 64. 31.

-->[D,V]=bdiag(B)

V =

- 0.4160251 - 1.8027756

- 0.5547002 - 3.6055513

D =

- 17. 0.

0. - 1.

-->E=diag(exp(diag(D)))

E =

4.140D-08 0.

0. 0.3678794

-->V*E*inv(V)

ans =

- 0.7357588 0.5518191

- 1.4715176 1.1036382

-->

-->