5 線形計算の常識: 行列が疎でも逆行列が疎とは限らない

微分方程式の数値シミュレーションに 現われる連立1次方程式の係数行列は、 ほとんどの場合、 成分に 0 が多いです。 そのような行列をそ疎 (sparse) であると言います。 疎行列はその性質をうまく利用することで、効率的な計算が可能になり、 大規模な問題が解けるようになる可能性があります (今だと未知数の個数が億のオーダーの連立1次方程式が解けたりします)。

行列が疎であっても、 その逆行列は疎であるとは限らない (大抵の場合、疎ではない) ことに注意する必要があります。

MATLAB のコマンド・メモ (1)

1:n	[1 2 3 $\cdots$ n]
a'	a の Hermite 共役 (実行列、実ベクトルの場合転置)
a.'	a の転置
eye(m,n)	単位行列
zeros(m,n)	零行列
ones(m,n)	成分がすべて $1$ の行列
rand(m,n)	成分が乱数の行列
diag()	対角行列 (を少しずらした行列)
命令1 ; 命令2	一行に複数の命令を書ける。 ``;'' があると表示が抑制される。

以下の例では、

$$A=\left( \begin{array}{cccccc} 2 & -1 & \\ -1 & 2 & -1\\ & \ddots &\ddots & \ddots \\ & & -1 & 2& -1\\ & & & -1 & 2 \end{array} \right)$$

について (実はある微分方程式を解く際に実際に現れる行列)、 連立1次方程式を解いたり、 逆行列を計算したりしています。

>> n=4
n = 4
>> 1:n
ans =

  1  2  3  4

>> (1:n)'*(1:n)
ans =

   1   2   3   4
   2   4   6   8
   3   6   9  12
   4   8  12  16

>> eye(n,n)
ans =

  1  0  0  0
  0  1  0  0
  0  0  1  0
  0  0  0  1

>> zeros(n,n)
ans =

  0  0  0  0
  0  0  0  0
  0  0  0  0
  0  0  0  0

>> ones(n-1,1)
ans =

  1
  1
  1

>> diag(ones(n-1,1),1)
ans =

  0  1  0  0
  0  0  1  0
  0  0  0  1
  0  0  0  0

>> diag(ones(n-1,1),-1)
ans =

  0  0  0  0
  1  0  0  0
  0  1  0  0
  0  0  1  0

>> a=2*eye(n,n)-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1)
a =

   2  -1   0   0
  -1   2  -1   0
   0  -1   2  -1
   0   0  -1   2

>> x=(1:n)'
x =

  1
  2
  3
  4

>> b=a*x
b =

  0
  0
  0
  5

>> ia=inv(a)
ia =

  0.80000  0.60000  0.40000  0.20000
  0.60000  1.20000  0.80000  0.40000
  0.40000  0.80000  1.20000  0.60000
  0.20000  0.40000  0.60000  0.80000

>> ia*a
ans =

   1.00000  -0.00000   0.00000   0.00000
   0.00000   1.00000   0.00000   0.00000
   0.00000   0.00000   1.00000   0.00000
   0.00000   0.00000   0.00000   1.00000

>> n=8;a=2*eye(n,n)-diag(ones(n-1,1),1)-diag(ones(n-1,1),-1)
a =

   2  -1   0   0   0   0   0   0
  -1   2  -1   0   0   0   0   0
   0  -1   2  -1   0   0   0   0
   0   0  -1   2  -1   0   0   0
   0   0   0  -1   2  -1   0   0
   0   0   0   0  -1   2  -1   0
   0   0   0   0   0  -1   2  -1
   0   0   0   0   0   0  -1   2

>> inv(a)
ans =

  0.88889  0.77778  0.66667  0.55556  0.44444  0.33333  0.22222  0.11111
  0.77778  1.55556  1.33333  1.11111  0.88889  0.66667  0.44444  0.22222
  0.66667  1.33333  2.00000  1.66667  1.33333  1.00000  0.66667  0.33333
  0.55556  1.11111  1.66667  2.22222  1.77778  1.33333  0.88889  0.44444
  0.44444  0.88889  1.33333  1.77778  2.22222  1.66667  1.11111  0.55556
  0.33333  0.66667  1.00000  1.33333  1.66667  2.00000  1.33333  0.66667
  0.22222  0.44444  0.66667  0.88889  1.11111  1.33333  1.55556  0.77778
  0.11111  0.22222  0.33333  0.44444  0.55556  0.66667  0.77778  0.88889

>>

$A$ が疎であっても、$A^{-1}$ の成分には $1$ つも 0 がないことを確認して下さい。

…あまりピンと来ないと思いますが、 これは実際上解決しないといけない壁なのです。

解決策はいくつかありますが、良く知られているのは、 逆行列の代わりに「LU分解」を使う、というものです。

任意の正則行列 $A$ に対して、

$$P A =L U,$$ (1)

$$P = ,\quad L= ,\quad U=$$ (2)

を満たす $P$ , $L$ , $U$ が存在する (付録 B.5を見よ)。

$A x=b$ であれば、 $P b=P A x=L U x$ であるから、

$$x=U^{-1}(L^{-1}(P b))$$

となり²、 $L^{-1}$ , $U^{-1}$ の掛け算は非常に簡単なので (付録 B.5.2 参照)、 $P$ , $L$ , $U$ が求まっていれば連立1次方程式は簡単に解ける。 (1), (2) を満たす $P$ , $L$ , $U$ を 求めることを $A$ を LU 分解するという。

MATLAB のコマンド・メモ (2)

[L u p]=lu(a)	a を LU 分解する。 $\texttt{p}\texttt{a}=\texttt{L} \texttt{u}$ を満たす p, L, u を 求める。
u\(L\(p*b))	$\texttt{a}\texttt{x}=\texttt{b}$ を解く。

>> [L u p]=lu(a)
L =

   1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
  -0.50000   1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000  -0.66667   1.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000  -0.75000   1.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000  -0.80000   1.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000  -0.83333   1.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000  -0.85714   1.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  -0.87500   1.00000

u =

   2.00000  -1.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   1.50000  -1.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   1.33333  -1.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   1.25000  -1.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   1.20000  -1.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   1.16667  -1.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   1.14286  -1.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   1.12500

p =
  1  0  0  0  0  0  0  0
  0  1  0  0  0  0  0  0
  0  0  1  0  0  0  0  0
  0  0  0  1  0  0  0  0
  0  0  0  0  1  0  0  0
  0  0  0  0  0  1  0  0
  0  0  0  0  0  0  1  0
  0  0  0  0  0  0  0  1

>> x=(1:n)'
x =

  1
  2
  3
  4
  5
  6
  7
  8

>> b=a*x
b =

  0
  0
  0
  0
  0
  0
  0
  9

>> u\(L\(p*b))
ans =

  1.0000
  2.0000
  3.0000
  4.0000
  5.0000
  6.0000
  7.0000
  8.0000

>>

$A$ を LU 分解したときの因子 $P$ , $L$ , $U$ が確かに 置換行列、下三角行列、上三角行列であることを確かめよう (このケースでは、$P$ は単位行列である)。 $L$ , $U$ が疎行列であることも確かめよう。

今回は使わなかったが、よく使われるコマンドを紹介しておく。

MATLAB のコマンド・メモ (3)

x(i)	x の第 i 成分
a(i,j)	a の第 $(\texttt{i},\texttt{j})$ 成分
det(a)	a の行列式
y'*x	縦ベクトル x, y の内積
x	x のノルム (成分の絶対値の二乗の和の平方根)
tril(a)	a の下三角部分
triu(a)	a の上三角部分
tic; コマンド; toc	コマンドの実行時間を計測
a(i,:)	a の第 i 行ベクトル
a(:,j)	a の第 j 列ベクトル
a(i1:i2,j1:j2)	a の第 i1〜i2 行、 第 j1〜j2 列の部分のブロック

余裕があれば試してみよう

1. 行列 a が大きいとき、 inv(a), det(a), a の実行時間を測って比較す る。
2. a の実行時間は、未知数の個数 $n$ につれどう変化するか。

    n=10
    for i=1:8  % この 8 をむやみに大きくしないこと
      a=rand(n,n); b=rand(n,1);
      tic; a\b; toc
      n=n*2
    end