3.2.3 実例: LU分解を用いて $A x=b$ を解く

(ある二日間の試行錯誤の記録。最初に MATLAB のやり方を踏襲しようとして， イマイチな結果となって，後でよりまともそうなやり方を見つけました， というストーリーです。お急ぎの方は最後の例だけ見て下さい。)

MATLABでは $A x=b$ はこうやって解く

n=10000;
a=rand(n,n);
[L,U,P]=lu(a);
x=ones(n,1);
b=a*x;
x2=U\(L\(P*b));
norm(x-x2)

あるいは

MATLABでは $A x=b$ はこうやって解く (置換をベクトルで扱う版)

n=10000;
a=rand(n,n);
[L,U,p]=lu(a,'vector');
x=ones(n,1);
b=a*x;
x2=U\(L\(b(p)));
norm(x-x2)

これを Scipy でやってみる。

>>> from numpy import *
>>> import scipy as sci
>>> import scipy.linalg

Numpy の関数はダイレクトに名前 (mat() とか linalg.solve() とか) で使える。 scipy の関数は sci. を先頭につけて (sci.linalg.lu() とか) 使える。 scipy の多くのサブパッケージは， 一々インポートしないと使えないものが多い。 ここでも scipy.linalg をインポートする。

>>> a=mat([[1,2],[3,4]])
>>> help(scipy.linalg.lu)
  help(sci.linalg.lu) でもOK
>>> P,L,U=sci.linalg.lu(a)
>>> P=mat(P)
>>> L=mat(L)
>>> U=mat(U)
>>> P
matrix([[ 0.,  1.],
        [ 1.,  0.]])
>>> L
matrix([[ 1.        ,  0.        ],
        [ 0.33333333,  1.        ]])
>>> U
matrix([[ 3.        ,  4.        ],
        [ 0.        ,  0.66666667]])
>>> P*L*U
matrix([[ 1.,  2.],
        [ 3.,  4.]])
>>> x=mat(ones((2,1)))
>>> b=a*x
>>> linalg.solve(U,linalg.solve(L,P.T*b))
matrix([[ 1.],
        [ 1.]])

さて大きい行列でやってみよう。

>>> n=10000
>>> a=mat(random.random((n,n)))
>>> P,L,U=sci.linalg.lu(a)
>>> P=mat(P)
>>> L=mat(L)
>>> U=mat(U)
>>> x=mat(ones((n,1)))
>>> b=a*x
>>> x2=linalg.solve(U,linalg.solve(L,P.T*b))
  (遅いね。もう一度解いているのかも。)
>>> linalg.norm(x-x2)
7.1042857390400037e-09

scipy.linalg.lu() は Scipy 用に書き下ろされたものだそうだ (LAPACK とかではなくて)。 これは出来がイマイチだ。 (MATLABでは，三角行列による割算 (L\    , U\    ) は高速に出来るので， この計算手順がまっとうなものとなる。)

(翌日)

はて？sci.linalg.lu_factor(), sci.linalg.lu_solve() を使うのか？

import numpy as np
import scipy as sci
import scipy.linalg

a=sci.mat([[1,2],[3,4]])
lu,piv=sci.linalg.lu_factor(a)
x=sci.mat([1,2]).T
b=a*x
x2=sci.linalg.lu_solve((lu,piv),b)

n=10000
a=sci.mat(np.random.random((n,n)))
lu,piv=sci.linalg.lu_factor(a)
x=np.ones((n,1))
b=a*x
x2=sci.linalg.lu_solve((lu,piv),b)
sci.linalg.norm(x-x2)

どうもそうみたい。こちらはずっと速い。