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について

:<math>
G 
\begin{pmatrix}
x \\
y \\
\end{pmatrix}
=
\begin{pmatrix}
* \\
0 \\
\end{pmatrix}
</math>


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

:<math>
A = 
\begin{pmatrix}
c & a \\
-s' & c \\
\end{pmatrix}
</math>

を返します。

以下に例を示します。

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

====[x, rcond] = inv (a)  ====
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　[Loadable Function]

====[x, rcond] = inverse (a)  ====
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　[Loadable Function]

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

====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 分解は，以下のような分解です。

:<math>
A = PHP^{T}
</math>
ここでP は正方ユニタリ行列（Unitary matrix）であり，H はUpper Hessenberg

:<math>(H_{i,j}=0,   
\forall{i} \geq j+1)</math>

です。

==== [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)（[[GNU Octave 2.1.x 日本語マニュアル/線形代数#aa = balance (a, opt)|balance]]を参照）ではなく、順序バランス(permutation balancing)を実行します。
出力引数の順序はMATLABとの互換性を持つよう選択されました

See also: [[GNU Octave 2.1.x 日本語マニュアル/線形代数#aa = balance (a, opt)|balance]], dare, [[GNU Octave 2.1.x 日本語マニュアル/線形代数#lambda = eig (a)|eig]], [[GNU Octave 2.1.x 日本語マニュアル/線形代数#s = schur (a)|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= U<sup>T</sup>AUを返します。ここでUはユニタリ行列で、Sは上三角行列です（U<sup>T</sup>Uは恒等です）。
A（とS）の固有値は、もし行列が実数の場合、Sの対角要素であり、その場合、実数シューア分解が計算されます。
行列Uは直交であり、Sはほとんどが対角線に沿った2ｘ2のサイズのブロックで構成されるブロック上三角行列です。
Sの対角要素はAとSの固有値です。（または、適切であれば2×2ブロックの固有値）
固有値は、オプションoptの値に応じて、対角線に沿って並べられます。
OPT ="a"は、負の実部を持つすべての固有値はSの先端ブロックに移動しなければならないことを示し（AREで使用される）、
OPT ="d"は、１より小さいすべての固有値は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を３つの行列U、Σ、 V<sup>H</sup>に分解することである。また、V<sup>H</sup>はVの随伴行列（共役転置行列）である。

:A = U Σ V<sup>H</sup>

結果を１つだけ返すように指定すると関数svdは特異値ベクトルΣを計算する。もし３項を返す様に指示すると、行列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]

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

:<math>
exp(A) = I + A + \frac{A^2}{2!}+ \frac{A^3}{3!}+ ...
</math>

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

参考文献：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近似は行列の有理多項式です。

:<math>
Dq(a)^{-1}Nq(a) 
</math>

そのテイラー級数は、上記のテイラー級数の最初の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 ]

[[Category:GNU Octave|2.1.x にほんこまにゆある せんけいたいすう]]
[[カテゴリ:線形代数学]]