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

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

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

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

以下の例では、

$$A=\left( \begin{array}{ccccc} 2 & -1 & & & \bigzerou \\ -1 & ... ...s & \ddots \\ & & -1 & 2& -1\\ \bigzerol & & & -1 & 2 \end{array} \right)$$

という三重対角行列 (tridiagonal matrix) について、 連立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.8000    0.6000    0.4000    0.2000
    0.6000    1.2000    0.8000    0.4000
    0.4000    0.8000    1.2000    0.6000
    0.2000    0.4000    0.6000    0.8000

>> ia*a

ans =

    1.0000   -0.0000   -0.0000    0.0000
         0    1.0000   -0.0000    0.0000
         0         0    1.0000    0.0000
         0         0         0    1.0000

>> 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.8889    0.7778    0.6667    0.5556    0.4444    0.3333    0.2222    0.1111
    0.7778    1.5556    1.3333    1.1111    0.8889    0.6667    0.4444    0.2222
    0.6667    1.3333    2.0000    1.6667    1.3333    1.0000    0.6667    0.3333
    0.5556    1.1111    1.6667    2.2222    1.7778    1.3333    0.8889    0.4444
    0.4444    0.8889    1.3333    1.7778    2.2222    1.6667    1.1111    0.5556
    0.3333    0.6667    1.0000    1.3333    1.6667    2.0000    1.3333    0.6667
    0.2222    0.4444    0.6667    0.8889    1.1111    1.3333    1.5556    0.7778
    0.1111    0.2222    0.3333    0.4444    0.5556    0.6667    0.7778    0.8889

>> 

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

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

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

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

を満たす $P$ , $L$ , $U$ が存在する (LU分解 -- 付録 D.5を見よ)。

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

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

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

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

[L u]=lu(a)	a を LU 分解する。 $\texttt{a}=\texttt{L} \texttt{u}$ を満たす L, u を求める。
u\(L\b)	$\texttt{a}\texttt{x}=\texttt{b}$ を解く ( $U^{-1}\left(L^{-1} b\right)$ を計算して)。
[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}$ を解く ( $U^{-1}\left(L^{-1}(P b)\right)$ を計算して)。

>> [L u]=lu(a)

L =

    1.0000         0         0         0         0         0         0         0
   -0.5000    1.0000         0         0         0         0         0         0
         0   -0.6667    1.0000         0         0         0         0         0
         0         0   -0.7500    1.0000         0         0         0         0
         0         0         0   -0.8000    1.0000         0         0         0
         0         0         0         0   -0.8333    1.0000         0         0
         0         0         0         0         0   -0.8571    1.0000         0
         0         0         0         0         0         0   -0.8750    1.0000


u =

    2.0000   -1.0000         0         0         0         0         0         0
         0    1.5000   -1.0000         0         0         0         0         0
         0         0    1.3333   -1.0000         0         0         0         0
         0         0         0    1.2500   -1.0000         0         0         0
         0         0         0         0    1.2000   -1.0000         0         0
         0         0         0         0         0    1.1667   -1.0000         0
         0         0         0         0         0         0    1.1429   -1.0000
         0         0         0         0         0         0         0    1.1250

>> 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\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 のコマンド・メモ (4)

x(i)	x の第 i 成分
a(i,j)	a の第 $(\texttt{i},\texttt{j})$ 成分
det(a)	a の行列式
y'*x	縦ベクトル x, y の内積
norm(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