『線形方程式解法の計算速度比較』

線形連立方程式を解くための方法にはさまざまな手法があります。Maple では、それら線形方程式系を解くための関数として通常の solve 関数以外に LinearSolve 関数が線形代数用パッケージ LinearAlgebra に用意されています。
ここでは、LinearSolve のオプションとして指定できるいくつかの算法で、計算速度がどのように変わるかを見ていきます。なお、当然のことながら計算速度は使用するCPUやメモリ量などによっても若干変わってくることに注意してください。

>    restart:

 LinearSolve 関数とは?

LinearSolve 関数は、LinearAlgebra (線形代数)用パッケージに収められています。

>    with(LinearAlgebra):

LinearSolve 関数は、次のような線形連立方程式の問題を解く場合に利用する関数です。

>    A := Matrix([[2, -5, 1],[-3, 4, 3],[0, -2, 1]]);

A := Matrix(%id = 12653704)

>    b := Vector([-3, -8, 1]);

b := Vector(%id = 12607048)

元の方程式は以下のようになっています。

>    A.Vector([x1,x2,x3]) = b;

Vector(%id = 12854912) = Vector(%id = 12607048)

LinearSolve 関数に係数行列とベクトルを渡すことで、この方程式の解を計算することができます。

>    sol := LinearSolve(A, b);

sol := Vector(%id = 12817420)

この解は、確かに方程式を満たしていることがわかります。

>    A.sol = b;

Vector(%id = 12901624) = Vector(%id = 12607048)

 大規模疎行列に対する計算時間の比較


LinearSolve 関数は、係数行列が実数(浮動小数点)の場合でも利用できるよう、さまざまなオプションが用意されています。この節では各解法に要する計算時間と結果の比較を行ってみましょう。
特に、係数行列は疎行列(三重行列、ただし非対称)を用います。

行列のサイズを表す変数を定義しておきます。

>    size := 1000:

乱数により行列を生成するために、乱数用ツールパッケージを読み込んでおきます。

>    with(RandomTools):

対角成分を生成しておきましょう。

>    diag := Generate(list(float(range=0.1..10.0), size)):

対角成分の上および下成分を生成しておきます。

>    semidiag1 := Generate(list(float(range=0.5..0.7), size-1)):
semidiag2 := Generate(list(float(range=0.2..0.4), size-1)):

係数行列を形成します。オプションに、sotrage=sparse を与え、疎行列用として宣言していることに注意してください。

>    A := Matrix([semidiag1, diag, semidiag2],
            storage=sparse, datatype=float[8], scan=band[1,2]);

A := Matrix(%id = 12626220)

定数項のベクトルを与えます。

>    b := Vector(Generate(list(float(range=0.0001..0.5), size)),
            datatype=float[8]);

b := Vector(%id = 13044216)

 オプション 『solve』 の場合

>    st := time():
sol_solve := LinearSolve(A, b, method='solve'):
time() - st;

73.325

>   

 オプション 『QR』 の場合

>    st := time():
sol_QR := LinearSolve(A, b, method='QR'):
time() - st;

3.585

>   

 オプション 『SparseLU』 の場合

>    st := time():
sol_SparseLU := LinearSolve(A, b, method='SparseLU'):
time() - st;

.20e-1

>   

 オプション 『SparseDirect』 の場合

>    st := time():
sol_SparseDirect := LinearSolve(A, b, method='SparseDirect'):
time() - st;

.20e-1

>   

 オプション 『hybrid』 の場合

※ オプション hybrid の効果を出すために、一時的に計算精度をあげてることに注意してください。

>    oldDigits := Digits: Digits := 30:
st := time():
sol_hybrid := LinearSolve(A, b, method='hybrid'):
time() - st;
Digits := oldDigits:

11.236

>   

 検算

各オプションで計算された結果を確かめてみましょう。

>    verifyError := proc(v::Vector)
  global A, b;
  local ret;

  Digits := 30:
  ret := max(op(map(abs, convert(A.v-b, list))));
  Digits := 10:
  return(ret);
end proc:

以下が各解法ごとの誤差の絶対値です。

>    errors :=
map(verifyError,
   [sol_solve, sol_QR, sol_SparseLU, sol_SparseDirect, sol_hybrid]
):

>    solvers :=
["solve", "QR", "SparseLU", "SparseDirect", "hybrid"]:
for i from 1 to nops(errors) do
  printf("\t%s = %e\n\n", solvers[i], errors[i]);
end do:

	solve = 4.580931e-08

	QR = 9.738641e-15

	SparseLU = 2.673344e-15

	SparseDirect = 2.673344e-15

	hybrid = 1.000000e-28

>   

計算速度も考慮すると、明らかに Sparse 型の処理算法が有利であることがわかります。

ここでの例題のように、特に疎な行列に対する線形方程式解法については、LinearSolve のオプションをうまく利用することで、高速にそしてより正確な解を得ることができるようになります。

閉じる