3 行列が疎でも逆行列が疎とは限らない

微分方程式の数値シミュレーションに 現われる連立1次方程式の係数行列は、 ほとんどの場合、 成分に $0$ が多い。 そのような行列をそ疎 (sparse) であるという。 疎行列はその性質をうまく利用することで、効率的な計算が可能になり、 大規模な問題が解けるようになる。 行列が疎であっても、 その逆行列は疎であるとは限らない (大抵の場合、疎ではない) ことに注意する。

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

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\... ...dots \\ & & -1 & 2& -1\\ & & & -1 & 2 \end{array} \right)$$

について、連立1次方程式を解いたり、 逆行列を計算したりしている。

samba03% octave
GNU Octave, version 2.0.17 (sparc-sun-solaris2.8).
Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 John W. Eaton.
This is free software with ABSOLUTELY NO WARRANTY.
For details, type `warranty'.

octave:1> n=4
n = 4
octave:2> 1:n
ans =

  1  2  3  4

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

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

octave:4> eye(n,n)
ans =

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

octave:5> zeros(n,n)
ans =

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

octave:6> ones(n-1,1)
ans =

  1
  1
  1

octave:7> diag(ones(n-1,1),1)
ans =

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

octave:8> diag(ones(n-1,1),-1)
ans =

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

octave:9> 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

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

  1
  2
  3
  4

octave:11> b=a*x
b =

  0
  0
  0
  5

octave:12> 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

octave:13> 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

octave:14> 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

octave:15> 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

octave:16>

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

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

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

$$P = \mbox{置換行列},\quad L=\mbox{下三角行列},\quad U=\mbox{上三角行列}$$ (2)

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

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

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

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

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

[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}$ を解く。

octave:16> [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

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

  1
  2
  3
  4
  5
  6
  7
  8

octave:18> b=a*x
b =

  0
  0
  0
  0
  0
  0
  0
  9

octave:19> u\(L\(p*b))
ans =

  1.0000
  2.0000
  3.0000
  4.0000
  5.0000
  6.0000
  7.0000
  8.0000

octave:20> quit
samba03%

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

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

Octave のコマンド・メモ (4)

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