(ある二日間の試行錯誤の記録。最初に MATLAB のやり方を踏襲しようとして, イマイチな結果となって,後でよりまともそうなやり方を見つけた, というストーリー。お急ぎの方は最後の例だけ見て下さい。)
MATLABでは はこうやって解く |
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では はこうやって解く (置換をベクトルで扱う版) |
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 |
>>> 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 とかではなくて)。 これは出来がイマイチみたい (, , を使って連立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)) |