ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

線形方程式

計算上の考慮事項

技術計算の中で最も重要な問題の 1 つは、線形方程式を解くことです。

行列表記では一般的な問題は次のような形式になります。[2 つの行列 A と b があるとき、Ax = b または xA = b の条件のいずれかを満たす一意の行列 x は存在しますか。]

1 行 1 列の簡単な例を考えます。たとえば次の方程式について考えます。

7x = 21

上記に一意の解は存在しますか。

もちろん存在します。この方程式には一意の解 x = 3 が存在します。解は除算から簡単に得られます。

x = 21/7 = 3.

解を求める場合、7 の逆数すなわち 7-1 = 0.142857... を計算し、この 7-1 に 21 を乗算する方法は通常 "使いません"。この方法では余計な計算ステップが必要となり、7-1 を有限小数値で表すと精度が低くなります。複数の未知数がある線形方程式にも同様の考え方で、MATLAB® は逆行列を使わずにこのような方程式を解きます。

標準の数学的な表記ではありませんが、MATLAB ではスカラーの場合と同じ除算記号を使用して一般的な連立方程式の解を記述します。2 つの除算記号、"スラッシュ" (/) と、"バックスラッシュ" (\) はそれぞれ MATLAB の 2 つの関数 mrdividemldivide に対応します。mrdividemldivide は、それぞれ係数行列の左または右に未知の行列が現れる場合に使用します。

x = b/A

行列方程式 xA = b の解です。

x = A\b

行列方程式 Ax = b の解です。

Ax = b または xA = b の方程式の両辺を A で割ると考えます。係数行列 A は常に「分母」 になります。

x = A\b に対する次元の整合性の条件は、2 つの行列 Ab の行数が同じであることです。その結果、解 xb と同じ列数になり、行数は A の列数と同じになります。x = b/A に対しては、行と列の処理が反転します。

実際には Ax = b の形式の線形方程式の方が xA = b の形式より頻繁に使用されます。そのためバックスラッシュはスラッシュよりも頻繁に使用されます。この節の以降ではバックスラッシュ演算子を中心に説明します。スラッシュ演算子の特性は次の等式から推定できます。

(b/A)' = (A'\b').

係数行列 A は正方である必要はありません。A が m 行 n 列の行列の場合、次の 3 つのケースが考えられます。

m = n

正方システム。厳密解が得られます。

m > n

過決定システム。最小二乗解を求めます。

m < n

劣決定システム。最大で m 個のゼロでない構成要素をもつ基本解を求めます。

mldivide アルゴリズム

mldivide 演算子は、さまざまなソルバーを使って各種の係数行列を処理します。実行するアルゴリズムは係数行列を調べて自動的に決められます。詳細は、関数 mldivide のリファレンス ページの「アルゴリズム」の節を参照してください。

一般解

線形方程式 Ax = b の一般解は、すべての解を示します。次のようにして一般解を求めることができます。

  1. 対応する同次方程式 Ax  =  0 を解きます。null コマンドを使用して「null(A)」と入力します。Ax = 0 の解空間の基底が返されます。解はすべて基底ベクトルの線形結合になります。

  2. 非同次方程式 Ax  =  b. の特殊解を導出します。

Ax = b の解はすべて、手順 2 の Ax = b の特殊解に手順 1 の基底ベクトルの線形結合を加えたものとして表現できます。

この節の以降では、MATLAB を使用して手順 2 の Ax = b の特殊解を導出する方法を説明します。

正方システム

一般的な例は、正方係数行列 A と右辺に列ベクトル b がある場合です。

特異でない係数行列

行列 A が特異でない場合、x = A\bb と同じサイズになります。たとえば、

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

A*xu と等価であることを確認してください。

Ab が正方で同じサイズの場合、x= A\b も同じサイズになります。

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

A*xb と等価であることを確認してください。

上記 2 つの例は、厳密に整数の解になります。これは係数行列がフル ランク行列 (特異でない) の pascal(3) であるためです。

特異な係数行列

正方行列 A が線形に独立した列をもっていない場合を「特異」といいます。A が特異な場合、Ax = b の解は存在しないかまたは一意ではありません。バックスラッシュ演算 A\b は、A が特異に近い場合は警告を出力し、特異である場合はエラーを出力します。

A が特異で Ax = b が解をもつ場合は、次のように入力することで一意ではない特殊解を導出できます。

P = pinv(A)*b

P は A の疑似逆行列です。Ax = b が厳密解をもたない場合、pinv(A) は最小二乗解を返します。

たとえば、

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

が特異であることは以下のように入力することでわかります。

rank(A)

ans =

     2

A はフル ランクではないため、いくつかの特異値がゼロに等しくなります。

    メモ:   関数 pinv を使用して係数行列が四角形になる方程式を解く方法の詳細は、疑似逆行列を参照してください。

厳密解-  b =[5;2;12] に対して、方程式 Ax = b は次のような厳密解をもちます。

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

以下のように入力すると、pinv(A)*b が厳密解であることを確かめることができます。

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

最小二乗解-  ただし、b = [3;6;0] の場合、Ax = b は厳密解をもちません。この例の場合、pinv(A)*b は最小二乗解を返します。以下のように入力すると、

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

元のベクトル b には戻りません。

Ax = b が厳密解であるかどうかは、拡大係数行列 [A b] を階段型行列に変換することによって判断できます。例では以下のようになります。

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

最終行は最後の要素以外はすべてゼロであるため、方程式は解をもちません。この例の場合、pinv(A) は最小二乗解を返します。

過決定システム

以下の例では、実験データでのさまざまな曲線近似において過決定システムがどのように使われているかを示します。

y は時刻 t のさまざまな値で測定され、観測値は次のとおりです。次のステートメントを使用して、データを入力し、テーブルで表示することができます。

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)
B = 

     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

このデータを指数的に減衰させる次の関数でモデル化します。

$$ y(t) =  c_1 + c_2 e^{-t} $$.

上記の方程式では、ベクトル y が、他の 2 つのベクトルの線形結合によって近似されなければならないことを表します。1 つはすべてが 1 になる定数ベクトルで、もう 1 つは成分 exp(-t) をもつベクトルです。未知の係数 $c_1$ $c_2$ は、モデル データの偏差の二乗和を最小にする最小二乗近似を行うことで計算します。2 つの未知数をもつ 6 つの方程式が、6 行 2 列の行列で表されます。

E = [ones(size(t)) exp(-t)]
E =

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

最小二乗解はバックスラッシュ演算子で求めます。

c = E\y
c =

    0.4760
    0.3413

別の表記では、データの最小二乗近似は次のように表せます。

$$ y(t) = 0.4760 + 0.3413 e^{-t} .$$

次のステートメントは、t を等間隔にインクリメントしてモデルを評価し、元のデータと共に結果をプロットします。

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

E*cy と厳密に一致しませんが、その差は元データに含まれる測定誤差よりもかなり小さくなります。

方形行列 A は線形に独立した列がない場合はランク落ちになります。A がランク落ちの場合、AX = B の最小二乗解は一意になりません。バックスラッシュ演算子 A\B は、A がランク落ちの場合に警告を出力し、システムが解をもたない場合は最小二乗解を作成し、システムが解を無限にもつ場合は基本解を作成します。

劣決定システム

この例では、劣決定システムの解が一意にならない理由について説明します。劣決定の線形システムは方程式より未知数の数が多くなります。MATLAB では行列の左除算演算により基本解を求めます。つまり、mn 列の係数行列には、最大 m 個の非ゼロの成分があります。

ここでは小さな乱数の例をあげます。

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)
R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8      

線形システム Rp = b には、4 つの未知数の中に 2 つの式があります。この行列係数に含まれる整数は小さいため、format コマンドを使用して有理数形式の解を表示するのが適切です。特殊解は次のステートメントで得られます。

format rat
p = R\b
p =

       0       
      17/7     
       0       
     -29/7    

R(:,2)R の最大ノルムをもつ列であるため、ゼロ以外の要素の 1 つは p(2) になります。R(:,2) を取り除くと R(:,4) が最大ノルムをもつため、他のゼロ以外の要素は p(4) になります。

劣決定システムの完全な一般解は、ヌル空間ベクトルの任意の線形結合に p を付加することで特性化され、有理数で表示するオプションにより関数 null を使用して求められます。

Z = null(R,'r')
Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1       

R*Z がゼロになります。残差 R*x - b は、以下に示す任意のベクトル x の場合は小さくなります。

x = p + Z*q.

Z の列はヌル空間ベクトルであるため、積 Z*q はこれらのベクトルの線形結合です。

Z*q=(x1x2)(uw)=ux1+wx2.

説明のため、任意の q を選択し、x を作成します。

q = [-2; 1];
x = p + Z*q;

残差のノルムを計算します。

format short
norm(R*x - b)
ans =

   2.6645e-15

連立線形方程式のマルチスレッド計算

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

  1. 関数は同時実行可能な部分へ簡単に分割できる処理を実行します。この分割部分は各処理間で通信をほとんど行わずに実行されるものでなければなりません。すなわちシーケンシャルな処理がほとんどないものを扱います。

  2. データを分割し、複数の実行スレッドを管理する時間を含めても同時実行する価値があるように、データ サイズは十分に大きいものを扱います。たとえば配列が数千個以上の要素を含む場合は、ほとんどの関数に対して高速化の効果があります。

  3. 処理がメモリに拘束されないものを扱います。すなわち処理時間の大部分がメモリのアクセス時間にならないものを扱います。一般的には、シンプルな処理より複雑な処理をする関数の方が、高速化の効果があります。

関数 invlscovlinsolvemldivide はマルチスレッド計算を有効にすると、大規模な倍精度配列 (10,000 個以上の要素) に対する計算速度が大幅に上昇します。

線形方程式を解く反復法

係数行列 A が大規模でスパース性がある場合、分解する方法は一般に効果がありません。反復法は近似解を生成できます。MATLAB は大規模でスパース性のある行列を扱う、いくつかの反復法を用意しています。

pcg

前処理を使用した共役勾配法。この方法はエルミート正定係数行列 A に適しています。

bicg

双共役傾斜法

bicgstab

双共役傾斜安定化法

bicgstabl

BiCGStab(l) 法

cgs

共役勾配二乗法

gmres

一般化最小残差法

lsqr

LSQR 法

minres

最小残差法。この方法はエルミート係数行列 A に適しています。

qmr

準最小残差法

symmlq

対称 LQ 法

tfqmr

転置なし準最小残差法

この情報は役に立ちましたか?