はじめに

この節で説明する 3 つの行列分解では、対角要素の上部または下部の要素がすべてゼロである "三角" 行列を使います。三角行列を含む線形方程式は、"前進代入" または "後退代入" のいずれかを使うと、簡単かつ高速に解くことができます。

コレスキー分解

コレスキー分解は、対称行列を三角行列と転置行列との積として表現します。

ここで、R は上三角行列です。

対称行列すべてがこの方法で分解できるわけではなく、このような分解をもつ行列は正定値と考えられます。これは、A の対角要素がすべて正で、非対角要素が「大きすぎない」ことを意味しています。パスカル行列を例に行列分解を行います。この章全体を通して行列の例 A は 3 行 3 列のパスカル行列を扱いますが、ここでは一時的に 6 行 6 列の行列を作成します。

A = pascal(6)

A =
       1     1     1     1     1     1
       1     2     3     4     5     6
       1     3     6    10    15    21
       1     4    10    20    35    56
       1     5    15    35    70   126
       1     6    21    56   126   252

A の要素は二項係数になります。各要素は上と左の要素の和になります。コレスキー分解は次のようになります。

R = chol(A)

R =
     1     1     1     1     1     1
     0     1     2     3     4     5
     0     0     1     3     6    10
     0     0     0     1     4    10
     0     0     0     0     1     5
     0     0     0     0     0     1

結果の要素は再び二項係数になります。R'*R が A に等しいということは、これが二項係数の積の和を含んでいることを示しています。

コレスキー分解を使用すると、次の線形方程式

を次の式で置き換えることができます。

バックスラッシュ演算子は三角行列を認識するため、MATLAB^® 環境で次のようにすばやく解くことができます。

x = R\(R'\b)

A が n 行 n 列の行列の場合、chol(A) の計算量は O(n³) になりますが、バックスラッシュの計算量は O(n²) です。

LU 分解

LU 分解またはガウスの消去法は任意の正方行列 A を、下三角行列と上三角行列の置換行列の積として表します。

ここで L は対角要素に 1 をもつ下三角行列を並べ替えたもので、U は上三角行列を並べ替えたものです。

並べ替えは理論上および計算上の面で必要になります。行列

は行を並べ替えないと、三角行列の積として表現することはできません。ただし、行列は次のようになります。

は三角行列の積として表すことができます。ε が小さい場合、因子の要素は大きくなり誤差も大きくなります。よってこの並べ替えは厳密には必要ありませんが、推奨されます。部分ピボットは L の要素の大きさを 1 以下に制限し、U の要素が A の要素より大きくならないようにします。

以下に例を示します。

[L,U] = lu(B)

L =
    1.0000         0         0
    0.3750    0.5441    1.0000
    0.5000    1.0000         0

U =
    8.0000    1.0000    6.0000
         0    8.5000   -1.0000
         0         0    5.2941

A を LU 分解すると、線形方程式

A*x = b

は次の式で高速に解くことができます。

x = U\(L\b)

行列式と逆行列は、次の関係を使って LU 分解から計算されます。

det(A) = det(L)*det(U)

および

inv(A) = inv(U)*inv(L)

det(A) = prod(diag(U)) を使用して行列式を計算することもできますが、行列式の符号が逆になる場合があります。

QR 分解

"直交" 行列または直交性の列がある行列とは、各列が単位長さをもち相互に直交関係になる実数行列です。Q が直交であれば、

ここで、I は単位行列です。

最も簡単な直交行列は 2 次元の座標回転変換です。

複素数行列に対応する概念は "ユニタリ" です。直交ユニタリ行列は長さと角度が保存され誤差が大きくならないため、数値計算に効果的です。

直交分解または QR 分解は任意の方形行列を、直交行列またはユニタリ行列と上三角行列の積で表現します。列の置換も含まれます。

または

ここで、Q は直交行列またはユニタリ行列、R は上三角行列、P は置換行列です。

QR 分解には、フル サイズとエコノミー サイズ、列置換をもつものともたないものの 4 種類あります。

過決定線形システムは、列よりも多くの行をもつ方形行列を含んでいます。すなわち、m 行 n 列で m > n です。フル サイズの QR 分解は、m 行 m 列の正方直交行列 Q と m 行 n 列の上三角行列 R を生成します。

C=gallery('uniformdata',[5 4], 0);
[Q,R] = qr(C)

Q =

    0.6191    0.1406   -0.1899   -0.5058    0.5522
    0.1506    0.4084    0.5034    0.5974    0.4475
    0.3954   -0.5564    0.6869   -0.1478   -0.2008
    0.3167    0.6676    0.1351   -0.1729   -0.6370
    0.5808   -0.2410   -0.4695    0.5792   -0.2207


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648
         0         0         0         0

多くの場合、Q の最後の m – n 列は R の下部にゼロを乗算することになるため必要ありません。したがってエコノミー サイズの QR 分解は、直交要素をもつ m 行 n 列の方形行列 Q と、正方の n 行 n 列の上三角行列 R を生成します。この 5 行 4 列の例ではあまり大きな影響はありませんが、非常に大きい方形行列の場合は、時間とメモリが節約されるため数値計算ではこの手法が非常に効果的になります。

[Q,R] = qr(C,0)	
Q =

    0.6191    0.1406   -0.1899   -0.5058
    0.1506    0.4084    0.5034    0.5974
    0.3954   -0.5564    0.6869   -0.1478
    0.3167    0.6676    0.1351   -0.1729
    0.5808   -0.2410   -0.4695    0.5792


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648

LU 分解とは異なり、QR 分解はピボットや置換は必要ありません。ただし 3 番目の出力引数を設定すると、オプションの列置換が出力され、特異性やランク落ちを調べる場合に便利です。分解の各ステップでは、分解していない残りの行列の列の中で最も大きなノルムをもつ列がそのステップでの基底として使われます。この計算方法により R の対角要素を小さい順に並べることができ、要素間の関係性を調べることによって各列の間に存在する線形依存性を明らかにすることができます。この簡単な例では、C の 2 列目は 1 列目のノルムよりも大きいノルムをもちます。そこで 2つの列を交換します。

[Q,R,P] = qr(C)

Q =
   -0.3522    0.8398   -0.4131
   -0.7044   -0.5285   -0.4739
   -0.6163    0.1241    0.7777

R =
  -11.3578   -8.2762
         0    7.2460
         0         0

P =
     0     1
     1     0

エコノミー サイズと列の置換を組み合わせると、3 番目の出力引数は置換行列ではなく置換ベクトルになります。

[Q,R,p] = qr(C,0)

Q =
   -0.3522    0.8398
   -0.7044   -0.5285
   -0.6163    0.1241

R =
  -11.3578   -8.2762
         0    7.2460


p =
     2     1

QR 分解では過決定の線形方程式を等価な三角行列に変換します。次の式をみてみましょう。

norm(A*x - b)

は以下と等価です。

norm(Q*R*x - b)

直交行列の乗算はユークリッド ノルムを保存するため、この式は以下と同じです。

norm(R*x - y)

ここで y = Q'*b です。R の最後の m-n 行はゼロであるため、この式は 2 つに分けられます。

norm(R(1:n,1:n)*x - y(1:n))

および

norm(y(n+1:m))

A がフル ランクの場合は x に対して解くことができ、最初の式はゼロになります。次に 2 番目の式が残差のノルムを指定します。A がフル ランクでない場合は、R の三角構造から最小二乗問題に対する基本解を計算します。

行列分解のマルチスレッド計算

MATLAB は多数の線形関数と要素ごとに計算する数値関数のマルチスレッド計算をサポートしています。これらの関数は自動的にマルチスレッドで実行されます。関数や式に対し、複数の CPU を使用して計算を高速化するにはさまざまな条件があります。

関数 lu と関数 qr を使うと、大規模な倍精度配列 (10,000 個以上の要素) に対する計算を大幅に高速化します。