20 線形代数 編集

この章では,Octave の線形代数に関する関数について述べています。 これら関数の多くに対してのリファレンスは,以下を参照してください:

Golub and Van Loan, Matrix Computations, 2ndEd., Johns Hopkins, 1989, and in Lapack Users’ Guide, SIAM, 1992.

20.1 基本的な行列関数 編集

aa = balance (a, opt) 編集

                                                          [Loadable Function]

[dd, aa] = balance (a, opt) 編集

                                                          [Loadable Function]

[cc, dd, aa, bb] = balance (a, b, opt) 編集

                                                          [Loadable Function]

[dd, aa] = balance (a)

は以下の計算結果を返します

aa = dd \ a * dd.

aaは、行と列のノルムの大きさにほぼ等しい行列であり、そして

dd = p * d,

pはpermute matrix(置換行列)であり、dは2のべき乗の対角行列です。 丸め誤差なく計算するするため、balanceをとることができます。 固有値計算の結果は、通常、最初にbalance行うことにより改善されます。

[cc, dd, aa, bb] = balance (a, b)

は、

aa = cc*a*dd

を返します。そして、

bb = cc*b*dd

です。aaとbbはほぼ同じ大きさの非ゼロ要素を持ちます、ここで ddと同様に、ccとddは代数的固有値問題のためのPermuted diagonal matrices(順序交換された対角行列)です。 固有値balancingオプションoptは、以下のように選択されます:

"N", "n" No balancing;

"N", "n":balancingを行いません。

引数はコピーされ、変換はそれぞれ個別にセットされます。

"P"、"p":可能であれば、固有値を分離するために引数の順序を変えます。

"S"、"s":固有値の計算精度を向上させるためにスケールします。

"B"、"b":B、bを用いて順序変更とスケールを、この順番で行います。

(aとb)の行/列は、順序変更から分離され、スケーリングされません。 これはデフォルトの動作です。 代数的固有値バランシングは、標準Lapackルーチンを使用しています。 一般化固有値問題のbalancingには、Wardのアルゴリズムを使用しています (SIAM Journal on Scientific and Statistical Computing, 1981).

cond (a) 編集

                                                           [Function File]

行列の2ノルム条件数(condition number)を計算します。 cond(a)はnorm(a) * norm (inv (a))と定義され,特異値分解を通して計算されます。

[d, rcond] = det (a) 編集

                                                           [Loadable Function]

Lapack を用いてa の行列式(determinant)を計算します。 必要であれば,reciprocal condition number の推定値を返します。

dmult (a, b) 編集

                                                           [Function File]

aが長さrows (b)のベクトルであれば,diag (a) * bを返す(この関数の方がずっと効率がよい)。

dot (x, y, dim) 編集

                                                           [Function File]

2つのベクトルの内積(dot product)を計算します。 xとyが行列であれば,first non-singleton dimension に沿って内積を計算します。 オプション引数dim を与えるならば,この次元に沿って内積を計算します。

lambda = eig (a) 編集

                                                           [Loadable Function]

[v, lambda] = eig (a) 編集

                                                           [Loadable Function]

行列の固有値(および固有ベクトル)を計算します。 この計算には,最初にHessenberg分解を行い,次に固有値が求まるまでSchur分解を行う。 もし望むなら,Schur分解をさらに実行することにより,固有ベクトルを計算します。

g = givens (x, y) 編集

                                                           [Loadable Function]

[c, s] = givens (x, y) 編集

                                                           [Loadable Function]

スカラxとyについて

 


となるような,2×2 の直交行列

 

を返します。

以下に例を示します。

givens (1, 1)
)
0.70711 0.70711
-0.70711 0.70711

[x, rcond] = inv (a) 編集

                                                           [Loadable Function]

[x, rcond] = inverse (a) 編集

                                                           [Loadable Function]

正方行列a の逆行列を計算します。 必要であれば,r条件数を推定します。 一方,条件数が小さいならば,悪条件という警告を表示します。

norm (a, p) 編集

                                                           [Function File]

行列a のp-ノルムを計算します。 もし2つめの引数を指定しないならば,p = 2と仮定します。

a が行列ならば:

p = 1 1-ノルム,a の絶対値の最大の列和
p = 2 a の最大の特異値
p = Inf 無限大ノルム,a の絶対値の最大の行和
p = "fro"
Frobenius ノルム,sqrt (sum (diag (a' * a)))

a がベクトルまたはスカラならば:

p = Inf max (abs (a)).
p = -Inf min (abs (a)).

その他aのp-ノルム,

(sum (abs (a) .^ p)) ^ (1/p)

null (a, tol) 編集

                                                           [Function File]

a の零空間(null space)の直交基底を返します。 零空間の次元は,a の特異値をとり,tol よりも大きくない。 tol を指定しないならば,それは以下のように計算します。

max (size (a)) * max (svd (a)) * eps

orth (a, tol) 編集

                                                           [Function File]

a のrange space の直交基底を返します。 range space の次元は,a の特異値をとり,tol よりも大きくない。 tol を指定しないならば,それは以下のように計算します。

max (size (a)) * max (svd (a)) * eps

pinv (x, tol) 編集

                                                           [Loadable Function]

x の疑似逆行列(一般化逆行列)を返します。 tol より小さな特異値は無視されます。 もし2 番めの引数を省略したならば,以下のように仮定します。

tol = max (size (x)) * sigma_max (x) * eps,

ここでsigma_max (x)は,x は最大の特異値です。

rank (a, tol) 編集

                                                           [Function File]

特異値分解を用いて,a の階数(rank)を計算します。 階数は,指定した基準値tol よりも大きなaの特異値とします。 2番めの引数を省略したならば,これは以下のように計算します。

tol = max (size (a)) * sigma(1) * eps;

ここでepsは計算機の精度であり,sigma(1)はa の最大の特異値です。

trace (a) 編集

                                                           [Function File]

a の対角和sum (diag (a))を計算します。

 20.2 行列の分解  編集

chol (a) 編集

                                                           [Loadable Function]

正定値である対称行列a を,以下のようにCholesky 分解します。

RTR = A

h = hess (a) 編集

                                                           [Loadable Function]

[p, h] = hess (a) 編集

                      [Loadable Function]

行列 a のHessenberg分解を行います。

通常,Hessenberg 分解は固有値計算の最初のステップとして使用されます。 しかし,他にも充分に応用できます。

(詳しくはGolub, Nash, and Van Loan, IEEETransactions on Automatic Control, 1979 を参照してください)

Hessenberg 分解は,以下のような分解です。

 

ここでP は正方ユニタリ行列(Unitary matrix)であり,H はUpper Hessenberg

 

です。

[l, u, p] = lu (a) 編集

                                                           [Loadable Function]

Lapack からのサブルーチンを用いてa のLU 分解を計算します。 その結果は,オプション戻り値p による並べ替え形式で返されます。 たとえば,行列a = [1, 2; 3, 4]を与えると,以下のようになります。

[l, u, p] = lu (a)
returns
l =
1.00000 0.00000
0.33333 1.00000
u =
3.00000 4.00000
0.00000 0.66667
p =
0 1
1 0

この行列は正方行列である必要はありません。

[q, r, p] = qr (a) 編集

                                                           [Loadable Function]

Lapack からのサブルーチンを用いてa のQR 分解を計算します。 その結果は,オプション戻り値p による並べ替え形式で返されます。 たとえば,行列a = [1, 2; 3, 4]を与えると,以下のようになります。

[q, r] = qr (a)
returns
q =
-0.31623 -0.94868
-0.94868 0.31623
r =
-3.16228 -4.42719
0.00000 -0.63246

qr分解は,過剰決定方程式系に対する最小二乗問題の解法の応用です。

min
x
kAx ! bk2

(過剰決定方程式系とはすなわち,A がtall,thin 行列です。)

QR 分解は,

QR = A

である。 ここでQ は直行行列であり,R は上三角行列です。

並べ替えQR 分解[q, r, p] = qr (a)は,rの対角要素がその大きさにおいて減少するよ うなQR 分解を形成します。 たとえば,行列a = [1, 2; 3, 4]を与えると,以下のようになります。

[q, r, p] = qr(a)
returns
q =
-0.44721 -0.89443
-0.89443 0.44721
r =
-4.47214 -3.13050
0.00000 0.44721
p =
0 1
1 0

並べ替えqr分解([q, r, p] = qr (a))は,span (a)の直交基底を構築します。

lambda = qz (a, b) 編集

                                                           [Loadable Function]

一般化固有値問題Ax = SBX、QZ分解。 この関数を呼び出すには、3つの方法があります:

1. lambda = qz(A,B)

(A!SB)の一般化固有値を計算します。


2. [AA, BB, Q, Z, V, W, lambda] = qz (A, B)

(A!SB)をQZ分解し、一般化固有ベクトルと一般化固有値を計算します。

AV = BV diag(,)
WTA = diag(,)WTB
AA = QTAZ;BB = QTBZ

 Q、Zあるいは、直交(ユニタリ)=I


3. [AA,BB,Z{, lambda}] = qz(A,B,opt)

フォーム[2]のように、例えば離散時間代数Riccati方程式の解のような一般固有値対の順序付けができます。

フォーム3は、複素数行列で利用可能ではありません。一般化直交行列Qだけでなく一般化固有ベクトルV、Wも計算しません


GEP束(GEP pencil)の固有値を選択するためのオプションopt。

更新された束の先頭のブロックは、それが満足する全ての固有値を含みます。

"N" = 順序付けされない(デフォルト)
"S" = small: 先頭ブロックは、|lambda| <= 1 をすべて含みます。
"B" = big: 先頭ブロックは、|lambda| >= 1 をすべて含みます
" - "=負の実数部:先頭のブロックは、左開半平面内のすべての固有値を持ちます。
"+" = nonnegative real part(非負の実部):先頭ブロックは右閉鎖半平面にあるすべての固有値を持ちます。

注:qzはスケーリング(scaling)(balanceを参照)ではなく、順序バランス(permutation balancing)を実行します。 出力引数の順序はMATLABとの互換性を持つよう選択されました

See also: balance, dare, eig, schur

[aa, bb, q, z] = qzhess (a, b) 編集

                                                           [Function File]

行列の組(a、b)の、Hessenberg三角分解( Hessenberg-triangular decomposition)を,計算し以下を返します:

aa = q * a * z,
bb = q * b * z,

q と z は直交しています. 例えば,

[aa, bb, q, z] = qzhess ([1, 2; 3, 4], [5, 6; 7, 8])
) aa = [ -3.02244, -4.41741; 0.92998, 0.69749 ]
) bb = [ -8.60233, -9.99730; 0.00000, -0.23250 ]
) q = [ -0.58124, -0.81373; -0.81373, 0.58124 ]
) z = [ 1, 0; 0, 1 ]

Hessenberg三角分解は、モウラーとスチュワート(Moler and Stewart)のQZ分解アルゴリズムの最初のステップです。 以下のアルゴリズムを採用しています:

Golub and Van Loan, Matrix Computations, 2nd edition.

s = schur (a) 編集

                                                           [Loadable Function]

[u, s] = schur (a, opt) 編集

                                                           [Loadable Function]

シューア分解(Schur decomposition)は正方行列の固有値を計算するために使用され、アプリケーションとして制御理論における代数的リカッチ方程式(algebraic Riccati equations)の解法があります(AREとDAREを参照)。 関数schurは常にS= UTAUを返します。ここでUはユニタリ行列で、Sは上三角行列です(UTUは恒等です)。 A(とS)の固有値は、もし行列が実数の場合、Sの対角要素であり、その場合、実数シューア分解が計算されます。 行列Uは直交であり、Sはほとんどが対角線に沿った2x2のサイズのブロックで構成されるブロック上三角行列です。 Sの対角要素はAとSの固有値です。(または、適切であれば2×2ブロックの固有値) 固有値は、オプションoptの値に応じて、対角線に沿って並べられます。 OPT ="a"は、負の実部を持つすべての固有値はSの先端ブロックに移動しなければならないことを示し(AREで使用される)、 OPT ="d"は、1より小さいすべての固有値はSの先端ブロックに移動しなければならないことを示し(DAREで使用される)、 そしてopt="u"は、 デフォルトで、固有値の順序付けが発生しないことを示しています。 Uの先頭のk列は、Sの固有値の先頭K個に対応し、常にA-不変部分空間を張ります。

s = svd (a) 編集

                                                           [Loadable Function]

[u, s, v] = svd (a) 編集

[Loadable Function]

行列Aの特異値分解(Singular value decomposition)を計算する。ここで、特異値分解とはAを3つの行列U、Σ、 VHに分解することである。また、VHはVの随伴行列(共役転置行列)である。

A = U Σ VH

結果を1つだけ返すように指定すると関数svdは特異値ベクトルΣを計算する。もし3項を返す様に指示すると、行列UとΣ(=s)、Vを計算する。例えば、

A=hilb (3);
s=svd (A)

returns

s =
1.4083189
0.1223271
0.0026873

そして

[u, s, v] = svd (A)
returns
u =
-0.82704 0.54745 0.12766
-0.45986 -0.52829 -0.71375
-0.32330 -0.64901 0.68867
s =
1.40832 0.00000 0.00000
0.00000 0.12233 0.00000
0.00000 0.00000 0.00269
v =
-0.82704 0.54745 0.12766
-0.45986 -0.52829 -0.71375
-0.32330 -0.64901 0.68867

もし、第二の引数が与えられると、svdはエコノミーサイズの分解を行い、uとvから不要な行と列を削除する。

[housv, beta, zer] = housh (x, j, z) 編集

                                                           [Function File]

j番目の列と同じにxを鏡映するためにHouseholder反射ベクトルhousvを計算する。すなわち、

(I - beta * housv * housv')x = e(j)

入力

x:ベクトル
j:ベクトルのインデックス
Z:ゼロのしきい値(通常は番号0でなければなりません)

出力: (Golub and Van Loanを参照してください)

beta:beta = 0なら、無反射が適用される必要がある(zerは0に設定)
housv:Householderベクトル

[u, h, nu] = krylov (a, v, k, eps1, pflg); 編集

[Function File]

ブロックkrylov部分空間の直交基底uを構築します。

[v a*v a^2*v ... a^(k+1)*v]

使用されているメソッド:直交性が失われるのを防ぐための鏡映(householder reflections)。

eps1:0と置くための閾値(:1e-12デフォルト)

pflg:行ピボットティングを使用するフラグ(数値の動作を向上させます)

0 [デフォルト]:ピボットティングなし。トリビアル・ヌルスペースが破損している場合、警告メッセージを出力します
1:ピボットティングを実行した出力

u:ブロックkrylov部分空間の直交基底

h:hessenberg行列は、vがベクトルであれば

a u = u h

、それ以外は無意味です

nu:スパンkrylov部分空間の次元(eps1に基づく)

もしbがベクトルで、k> m-1の場合には、krylov関数は、

h =hessenberg行列分解

を返します。

Reference: Hodel and Misra, "Partial Pivoting in the Computation of Krylov Subspaces ", to be submitted to Linear Algebra and its Applications

20.3 行列に対する関数 編集

expm (a) 編集

                                                           [Loadable Function]

無限次元のテイラー級数で定義される行列のべき乗を返す。

 

テイラー級数は、行列のべき乗を計算する方法としては使われていません。詳細は以下の文献の参照してください。

参考文献:Moler and Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, SIAM Review, 1978.

このルーチンは、ウォードのダイアゴナルPad'e近似法を3ステップの前処理演算と共に使用しています(SIAM Journal on Numerical Analysis, 1977)。 ディアゴナルPad'e近似は行列の有理多項式です。

 

そのテイラー級数は、上記のテイラー級数の最初の2Q+1項にマッチします。 同じ前処理演算ステップを行うテイラー級数の直接評価は、Dq(a)が悪条件である場合、Pad'e近似の代わりに、望ましいかもしれません。

logm (a) 編集

                                                           [Function File]

正方行列aの、行列の対数を計算します。

ノート:現在のところ、インプリメントに固有値のエクスパンションを使用している、よりロバストにする必要がある。

[result, error_estimate] = sqrtm (a) 編集

                                                           [Loadable Function]

正方行列aの、行列の平方根を計算します。 参考文献:

Nicholas J. Higham. A new sqrtm for MATLAB. Numerical Analysis Report No. 336, Manchester Centre for Computational Mathematics, Manchester, England,January 1999.

kron (a, b) 編集

                                                           [Function File]

2 つの行列のKronecker 積を計算します。 ブロック同士の積は,以下のように定義します。

x = [a(i, j) b]

以下に例を示す。

kron (1:4, ones (3, 1))
)
1 2 3 4
1 2 3 4
1 2 3 4

x = syl (a, b, c) 編集

                                                           [Loadable Function]

Sylvester方程式を解きます。

AX + XB + C = 0

標準Lapackサブルーチンを用います。 以下に、例題を示します。、

syl ([1, 2; 3, 4], [5, 6; 7, 8], [9, 10; 11, 12])
) [ -0.50000, -0.66667; -0.66667, -0.50000 ]