Matlab/Scilabにはexpm()という行列指数関数が準備されている。Matlab/Scilabのexpm()では、スケーリングと正方化(scaling and squaring)されたパデ近似が使われている。
行列を対角化可能な場合は固有値と固有ベクトルに分解するmatrix decomposition methodがもっとも効率的。matrix decomposition methodはどこかで見たことがあると思ったら、中村先生のコンピュテーショナル・ファイナンスの2章の『至高のgenerator作り』でやったやつだな。あのときは1年後の格付け推移行列を求めた。
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
-->//例1
-->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
-->
-->//例2
-->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
-->
-->