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

G{\begin{pmatrix}x\\y\\\end{pmatrix}}={\begin{pmatrix}*\\0\\\end{pmatrix}}

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

A={\begin{pmatrix}c&a\\-s'&c\\\end{pmatrix}}

を返します。

以下に例を示します。

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

A=PHP^{T}

ここでP は正方ユニタリ行列（Unitary matrix）であり，H はUpper Hessenberg

(H_{i,j}=0,\forall {i}\geq j+1)

です。

[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との互換性を持つよう選択されました

[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^TAUを返します。ここでUはユニタリ行列で、Sは上三角行列です（U^TUは恒等です）。 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^Hに分解することである。また、V^HはVの随伴行列（共役転置行列）である。

A = U Σ V^H

結果を１つだけ返すように指定すると関数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]

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

exp(A)=I+A+{\frac {A^{2}}{2!}}+{\frac {A^{3}}{3!}}+...

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

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

Dq(a)^{-1}Nq(a)

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

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 ]

GNU Octave 2.1.x 日本語マニュアル/線形代数

20 線形代数

20.1 基本的な行列関数

aa = balance (a, opt)

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

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

cond (a)

[d, rcond] = det (a)

dmult (a, b)

dot (x, y, dim)

lambda = eig (a)

[v, lambda] = eig (a)

g = givens (x, y)

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

[x, rcond] = inv (a)

[x, rcond] = inverse (a)

norm (a, p)

null (a, tol)

orth (a, tol)

pinv (x, tol)

rank (a, tol)

trace (a)

20.2 行列の分解

chol (a)

h = hess (a)

[p, h] = hess (a)

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

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

lambda = qz (a, b)

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

s = schur (a)

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

s = svd (a)

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

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

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

20.3 行列に対する関数

expm (a)

logm (a)

[result, error_estimate] = sqrtm (a)

kron (a, b)

x = syl (a, b, c)

　20.2 行列の分解　