5.3.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, L, U は array である。)

>>> 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.]])

scipy.linalg.lu() は Scipy 用に書き下ろされたものだそうだ (LAPACK とかではなくて)。 これは出来がイマイチみたい ($P$, $L$, $U$ を使って連立1次方程式を解くとき、 あまり速く解けない)。

sci.linalg.lu_factor(), sci.linalg.lu_solve() というのもある。

>>> import scipy as sci
>>> import scipy.linalg
>>> a=sci.mat([[1,2],[3,4]])
>>> lu,piv=sci.linalg.lu_factor(a)
>>> lu
array([[ 3.        ,  4.        ],
       [ 0.33333333,  0.66666667]])
>>> piv
array([1, 1], dtype=int32)

(この lu, piv を見て「なるほど」という感じがする。)

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

大きい問題を解いてみる。

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

n=10000
print(n,"次の乱数行列生成")
a=sci.mat(np.random.random((n,n)))
print("LU分解")
lu,piv=sci.linalg.lu_factor(a)
x=np.ones((n,1))
print("掛け算")
b=a*x
print("方程式を解く")
x2=sci.linalg.lu_solve((lu,piv),b)
print("誤差=",sci.linalg.norm(x-x2))

こちらはマルチコアをちゃんと使って速い (time で測って 580% という数値が出た)。