GNU Octave 2.1.x 日本語マニュアル/線形代数
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 ]