Next: Functions and Variables for lapack, Previous: lapack, Up: lapack [Contents][Index]
lapackは
SLATECプロジェクトから得られるようなFortranライブラリLAPACKの
(プログラム f2cを介した) Common Lisp翻訳です。
(訳者注意書き: lapackを使用するには、load("lapack"); load("eigensys");を実行してください。load("lapack")には、初回だけコンパイルで時間がかかるかもしれません。)
Previous: Introduction to lapack, Up: lapack [Contents][Index]
行列Aの固有値と、オプションで固有ベクトルを計算します。 Aの要素はすべて整数か浮動小数点数でなければいけません。 Aは平方(行と列が同じ数)でなければいけません。 Aは対称であってもなくてもいいです。
dgeev(A)はAの固有値だけを計算します。
dgeev(A, right_p, left_p)は
Aの固有値と、
right_p = trueの時、右固有ベクトル、
left_p = trueの時、左固有ベクトルを計算します。
3項目のリストが返されます。
最初の項目は固有値のリストです。
二番目の項目はfalseか右固有ベクトルの行列です。
三番目の項目はfalseか左固有ベクトルの行列です。
右固有ベクトルv(j) (右固有ベクトル行列のj番目の列)は
A . v(j) = lambda(j) . v(j)
を満たします。ここでlambda(j)は対応する固有値です。
左固有ベクトルv(j) (左固有ベクトル行列のj番目の列)は
u(j)**H . A = lambda(j) . u(j)**H
を満たします。ここでu(j)**Hはu(j)の共役転置を意味します。
Maxima関数ctransposeが共役転置を計算します。
計算された固有ベクトルは、 Euclideanノルムが1に等しく、 最大成分の虚部が0になるように規格化されます。
例:
(%i1) load ("lapack")$
(%i2) fpprintprec : 6;
(%o2) 6
(%i3) M : matrix ([9.5, 1.75], [3.25, 10.45]);
[ 9.5 1.75 ]
(%o3) [ ]
[ 3.25 10.45 ]
(%i4) dgeev (M);
(%o4) [[7.54331, 12.4067], false, false]
(%i5) [L, v, u] : dgeev (M, true, true);
[ - .666642 - .515792 ]
(%o5) [[7.54331, 12.4067], [ ],
[ .745378 - .856714 ]
[ - .856714 - .745378 ]
[ ]]
[ .515792 - .666642 ]
(%i6) D : apply (diag_matrix, L);
[ 7.54331 0 ]
(%o6) [ ]
[ 0 12.4067 ]
(%i7) M . v - v . D;
[ 0.0 - 8.88178E-16 ]
(%o7) [ ]
[ - 8.88178E-16 0.0 ]
(%i8) transpose (u) . M - D . transpose (u);
[ 0.0 - 4.44089E-16 ]
(%o8) [ ]
[ 0.0 0.0 ]
行列 AのQR分解します。 All elements of Aのすべての要素は整数か浮動小数点数でなければいけません。 Aは行と列の数は同じかもしれませんし違うかもしれません。
2つの項目のリストを返します。
一番目の項目は行列 Qで、それはAと同じ行数を持つ平方正規直交行列です。
二番目の項目は行列 Rで、それはAtお同じサイズで、
対角以下のすべての要素がが零に等しいものです。
積 Q . Rは(浮動小数点の丸め誤差を除いて)Aに等しい。
ここで "."は非可換乗算演算子です。
(%i1) load ("lapack") $
(%i2) fpprintprec : 6 $
(%i3) M : matrix ([1, -3.2, 8], [-11, 2.7, 5.9]) $
(%i4) [q, r] : dgeqrf (M);
[ - .0905357 .995893 ]
(%o4) [[ ],
[ .995893 .0905357 ]
[ - 11.0454 2.97863 5.15148 ]
[ ]]
[ 0 - 2.94241 8.50131 ]
(%i5) q . r - M;
[ - 7.77156E-16 1.77636E-15 - 8.88178E-16 ]
(%o5) [ ]
[ 0.0 - 1.33227E-15 8.88178E-16 ]
(%i6) mat_norm (%, 1);
(%o6) 3.10862E-15
線形方程式 A x = bの 解 xを計算します。 ここで、 Aは平方行列、 bはAと同じ数の行と任意の長さの列を持つ行列です。 戻り値 xは bと同じサイズです。
Aと bの要素は
floatを介して実の浮動小数点数に評価されなければいけません;
従って、要素は任意の数値型か、数値定数のシンボルか、
浮動小数点に評価される式であり得ます。
xの要素はいつも浮動小数点数です。
すべての算術は浮動小数演算として実行されます。
dgesvは
AのLU分解を介して
解を計算します。
例:
dgesvは
線形方程式 A x = bの解を計算します。
(%i1) A : matrix ([1, -2.5], [0.375, 5]);
[ 1 - 2.5 ]
(%o1) [ ]
[ 0.375 5 ]
(%i2) b : matrix ([1.75], [-0.625]);
[ 1.75 ]
(%o2) [ ]
[ - 0.625 ]
(%i3) x : dgesv (A, b);
[ 1.210526315789474 ]
(%o3) [ ]
[ - 0.215789473684211 ]
(%i4) dlange (inf_norm, b - A.x);
(%o4) 0.0
bはAと同じ数の行と任意の長さの列を持つ行列です。 xは bと同じサイズです。
(%i1) A : matrix ([1, -0.15], [1.82, 2]);
[ 1 - 0.15 ]
(%o1) [ ]
[ 1.82 2 ]
(%i2) b : matrix ([3.7, 1, 8], [-2.3, 5, -3.9]);
[ 3.7 1 8 ]
(%o2) [ ]
[ - 2.3 5 - 3.9 ]
(%i3) x : dgesv (A, b);
[ 3.103827540695117 1.20985481742191 6.781786185657722 ]
(%o3) [ ]
[ - 3.974483062032557 1.399032116146062 - 8.121425428948527 ]
(%i4) dlange (inf_norm, b - A . x);
(%o4) 1.1102230246251565E-15
Aと bの要素は 実の浮動小数点数に評価されなければいけません;
(%i1) A : matrix ([5, -%pi], [1b0, 11/17]);
[ 5 - %pi ]
[ ]
(%o1) [ 11 ]
[ 1.0b0 -- ]
[ 17 ]
(%i2) b : matrix ([%e], [sin(1)]);
[ %e ]
(%o2) [ ]
[ sin(1) ]
(%i3) x : dgesv (A, b);
[ 0.690375643155986 ]
(%o3) [ ]
[ 0.233510982552952 ]
(%i4) dlange (inf_norm, b - A . x);
(%o4) 2.220446049250313E-16
特異値から成る行列 Aの特異値分解(SVD)を計算します。 オプションで左および右特異ベクトルを取ります。
Aの要素はすべて整数か浮動小数点数でなければいけません。 Aは(行と列が同じ数の)平方かもしれませんし、そうでないかもしれません。
mを Aの行数、nを列数とします。 Aの特異値分解は A = U . Sigma . V^T のような3つの行列 U, Sigma, V^Tから構成されます。 ここで、 Uは m掛けmのユニタリ行列、 Sigmaは m掛けnの対角行列、 V^Tは n掛けnのユニタリ行列です。
sigma[i]を
Sigmaの対角要素、すなわち、 Sigma[i, i] = sigma[i]と
します。
要素 sigma[i]は Aのいわゆる特異値です;
これらは実数で、非負で、降順で返されます。
Uと Vの最初の min(m, n)列は
Aの左と右特異ベクトルです。
dgesvdは、V自身ではなくVの転置を返すことに注意してください。
dgesvd(A)は Aの特異値だけを計算します。
dgesvd(A, left_p, right_p)は
Aの特異値と、
left_p = trueの時、左特異ベクトル、
right_p = trueの時、右特異ベクトルを計算します。
3つの項目のリストが返されます。
一つ目の項目は特異値のリストです。
二つ目の項目は falseか、左特異ベクトルの行列です。
三つ目の項目は falseか、右特異ベクトルの行列です。
例:
(%i1) load ("lapack")$
(%i2) fpprintprec : 6;
(%o2) 6
(%i3) M: matrix([1, 2, 3], [3.5, 0.5, 8], [-1, 2, -3], [4, 9, 7]);
[ 1 2 3 ]
[ ]
[ 3.5 0.5 8 ]
(%o3) [ ]
[ - 1 2 - 3 ]
[ ]
[ 4 9 7 ]
(%i4) dgesvd (M);
(%o4) [[14.4744, 6.38637, .452547], false, false]
(%i5) [sigma, U, VT] : dgesvd (M, true, true);
(%o5) [[14.4744, 6.38637, .452547],
[ - .256731 .00816168 .959029 - .119523 ]
[ ]
[ - .526456 .672116 - .206236 - .478091 ]
[ ],
[ .107997 - .532278 - .0708315 - 0.83666 ]
[ ]
[ - .803287 - .514659 - .180867 .239046 ]
[ - .374486 - .538209 - .755044 ]
[ ]
[ .130623 - .836799 0.5317 ]]
[ ]
[ - .917986 .100488 .383672 ]
(%i6) m : length (U);
(%o6) 4
(%i7) n : length (VT);
(%o7) 3
(%i8) Sigma:
genmatrix(lambda ([i, j], if i=j then sigma[i] else 0),
m, n);
[ 14.4744 0 0 ]
[ ]
[ 0 6.38637 0 ]
(%o8) [ ]
[ 0 0 .452547 ]
[ ]
[ 0 0 0 ]
(%i9) U . Sigma . VT - M;
[ 1.11022E-15 0.0 1.77636E-15 ]
[ ]
[ 1.33227E-15 1.66533E-15 0.0 ]
(%o9) [ ]
[ - 4.44089E-16 - 8.88178E-16 4.44089E-16 ]
[ ]
[ 8.88178E-16 1.77636E-15 8.88178E-16 ]
(%i10) transpose (U) . U;
[ 1.0 5.55112E-17 2.498E-16 2.77556E-17 ]
[ ]
[ 5.55112E-17 1.0 5.55112E-17 4.16334E-17 ]
(%o10) [ ]
[ 2.498E-16 5.55112E-17 1.0 - 2.08167E-16 ]
[ ]
[ 2.77556E-17 4.16334E-17 - 2.08167E-16 1.0 ]
(%i11) VT . transpose (VT);
[ 1.0 0.0 - 5.55112E-17 ]
[ ]
(%o11) [ 0.0 1.0 5.55112E-17 ]
[ ]
[ - 5.55112E-17 5.55112E-17 1.0 ]
行列 Aのノルムもしくはノルムのような関数を計算します。
maxmax(abs(A(i, j)))を計算します。 ここで iと jはそれぞれ行と列を行き渡ります。 この関数は適切な行列ノルムではないことに注意してください。
one_normAの L[1]ノルム、 すなわち、それぞれの列の要素の絶対値の和の最大値 を計算します。
inf_normAの L[inf]ノルム、 すなわち、それぞれの行の要素の絶対値の和の最大値 を計算します。
frobeniusAのFrobeniusノルム、すなわち、 行列要素の平方の和の平方根 を計算します。
2つの行列の積を計算します。オプションで積を三つ目の行列に足し算します。
最も簡単な形式では、
dgemm(A, B)は
2つの実行列 Aと Bの積を計算します。
二番目の形式では、
dgemmは
alpha * A * B + beta * C
を計算します。
ここで A, B, Cは
適当なサイズの実行列であり、
alphaと betaは実数です。
オプションで、 Aと/もしくは Bは
積を計算する前に転置を取ることができます。
追加のパラメータはオプションのキーワード引数で指定できます:
キーワード引数はオプションで、
どんな順番でも指定できます。
それらはすべて、形式 key=valを取ります。
キーワード引数は以下の通りです:
C足すべき行列 C。
デフォルトは falseであり、行列を足さないことを意味します。
alphaAと Bの積がこの値に掛けられます。 デフォルトは1です。
betaもし行列 Cが与えられたら、 この値は、足される前にCに掛けられます。 デフォルト値は0で、これは, たとえCが与えられても Cが足されないことを意味します。 故に、必ずbetaに零でない値を指定してください。
transpose_aもし trueなら、
Aの代わりにAの転置が積に使われます。
デフォルトは falseです。
transpose_bもし trueなら
Bの代わりにBの転置が積に使われます。
デフォルトは falseです。
(%i1) load ("lapack")$
(%i2) A : matrix([1,2,3],[4,5,6],[7,8,9]);
[ 1 2 3 ]
[ ]
(%o2) [ 4 5 6 ]
[ ]
[ 7 8 9 ]
(%i3) B : matrix([-1,-2,-3],[-4,-5,-6],[-7,-8,-9]);
[ - 1 - 2 - 3 ]
[ ]
(%o3) [ - 4 - 5 - 6 ]
[ ]
[ - 7 - 8 - 9 ]
(%i4) C : matrix([3,2,1],[6,5,4],[9,8,7]);
[ 3 2 1 ]
[ ]
(%o4) [ 6 5 4 ]
[ ]
[ 9 8 7 ]
(%i5) dgemm(A,B);
[ - 30.0 - 36.0 - 42.0 ]
[ ]
(%o5) [ - 66.0 - 81.0 - 96.0 ]
[ ]
[ - 102.0 - 126.0 - 150.0 ]
(%i6) A . B;
[ - 30 - 36 - 42 ]
[ ]
(%o6) [ - 66 - 81 - 96 ]
[ ]
[ - 102 - 126 - 150 ]
(%i7) dgemm(A,B,transpose_a=true);
[ - 66.0 - 78.0 - 90.0 ]
[ ]
(%o7) [ - 78.0 - 93.0 - 108.0 ]
[ ]
[ - 90.0 - 108.0 - 126.0 ]
(%i8) transpose(A) . B;
[ - 66 - 78 - 90 ]
[ ]
(%o8) [ - 78 - 93 - 108 ]
[ ]
[ - 90 - 108 - 126 ]
(%i9) dgemm(A,B,c=C,beta=1);
[ - 27.0 - 34.0 - 41.0 ]
[ ]
(%o9) [ - 60.0 - 76.0 - 92.0 ]
[ ]
[ - 93.0 - 118.0 - 143.0 ]
(%i10) A . B + C;
[ - 27 - 34 - 41 ]
[ ]
(%o10) [ - 60 - 76 - 92 ]
[ ]
[ - 93 - 118 - 143 ]
(%i11) dgemm(A,B,c=C,beta=1, alpha=-1);
[ 33.0 38.0 43.0 ]
[ ]
(%o11) [ 72.0 86.0 100.0 ]
[ ]
[ 111.0 134.0 157.0 ]
(%i12) -A . B + C;
[ 33 38 43 ]
[ ]
(%o12) [ 72 86 100 ]
[ ]
[ 111 134 157 ]