分散配列を使用した反復法による連立線形方程式の求解

ライブ スクリプトを開く

大規模な数学計算の場合、反復法は直接法より効率的なことがあります。この例では、分散配列を使用して、反復法によって Ax=b の形式の連立線形方程式を並列的に求解する方法を説明します。

この例は、分散配列を使用した直接法による連立線形方程式の求解で説明されているトピックの続きです。mldivideで実装される直接ソルバー メソッドは、連立線形方程式を分散して並列で求解するために使用できますが、大規模でスパースな特定の連立方程式では効率的でない可能性があります。反復法では、初期推定から一連の解を生成し、いくつかのステップを経て最終結果に収束させます。これらのステップの方が、解を直接計算する場合よりも計算量が少なくなることがあります。

分散配列は、クライアント ワークスペースのデータをローカル マシンまたはクラスターの並列プールに分散します。各ワーカーは配列の一部をそのメモリに保存しますが、他のワーカーと通信して配列のすべてのセグメントにアクセスすることもできます。分散配列は、完全な行列、スパース行列などのさまざまなデータ型を格納できます。

この例では、関数pcgを使用して、共役勾配法と前処理付き共役勾配法によって、大規模な連立線形方程式を解く方法を説明します。反復法は、密行列とスパース行列のいずれとも使用できますが、スパース行列方程式の場合により効率的です。

スパース行列を使用した連立線形方程式の定義

distributed関数を使用すると、ソフトウェアは既定のクラスター設定を使用して、並列プールを自動的に起動します。この例では、関数galleryからの Wathen 行列を使用します。この行列は、スパースで対称の乱数行列であり、その全体の次元は N=3n2+4n+1 です。

n = 400;
A = distributed(gallery("wathen",n,n));
N = 3*n^2+4*n+1
N = 
481601

ここで、右辺のベクトル b を定義できます。この例で、b は A の行の和として定義されています。これにより、Ax=b の厳密解は xexact=[1,...,1]T の形式になります。

b = sum(A,2);

sum は分散配列を処理するため、b も分散され、そのデータは並列プールのワーカーのメモリに保存されます。最後に、反復法を使用して求めた解と比較するために、厳密解を定義できます。

xExact = ones(N,1,"distributed");

共役勾配法による連立線形方程式の求解

MATLAB 関数 pcg は、共役勾配 (CG) 法を提供します。これは、ステップごとに解を改善させながら、x の一連の近似解を反復的に生成します。 

[xCG_1,flagCG_1,relres_CG1,iterCG_1,resvecCG_1] = pcg(A,b);

方程式を解いたら、得られた結果 xCG_1 の各要素と予測値 xExact との誤差をチェックできます。計算結果の誤差は比較的大きなものです。

errCG_1 = abs(xExact-xCG_1);

figure(1)
hold off
semilogy(errCG_1,"o");
title("System of Linear Equations with Sparse Matrix");
ylabel("Absolute Error");
xlabel("Element in x");

反復的な計算は、一連の近似解が特定の許容誤差に収束するか、反復ステップの最大回数に達すると終了します。pcg は、分散配列とクライアント上の配列のいずれにも同じ既定の設定を使用します。

  • 既定の最大許容誤差は、10-6 です。

  • 反復ステップの既定の最大回数は 20、または係数行列 A の次数 (20 未満の場合) です。

関数 pcg は 2 番目の出力引数として収束フラグも返します。収束フラグは、計算された解が目標の許容誤差に収束したかどうかなど、得られた結果についてより多くの情報を提供します。たとえば、値 0 は、解が適切に収束したことを示します。

flagCG_1
flagCG_1 = 
1

この例では、解が既定の最大反復回数内で収束せず、大きな誤差が出力されます。

収束の確率を高めるために、許容誤差と反復ステップの最大回数の設定をカスタマイズできます。

tolerance = 1e-12;
maxit = N;

tCG = tic;
[xCG_2,flagCG_2,relresCG_2,iterCG_2,resvecCG_2] = pcg(A,b,tolerance,maxit);
tCG = toc(tCG);

flagCG_2
flagCG_2 = 
0

カスタム設定を使用すると、解は収束します。この解の絶対誤差は、前述の解より改善されています。

errCG_2 = abs(xExact-xCG_2);
figure(2)
hold off
semilogy(errCG_1,"o");
hold on
semilogy(errCG_2,"d");
title("Comparison of Absolute Error");
ylabel("Absolute Error");
xlabel("Element in x");
legend("Default tolerance and iterations", ...
    "Improved tolerance and iterations");
hold off

pcg メソッドは、各反復ステップで残差ノルム norm(b-A*x)/norm(b) のベクトルも返します。相対残差ノルムは、連続する反復ステップ間の精度の比を示します。反復処理中の残差の変化は、カスタム設定を使用しない場合に解が収束しなかった理由を理解する上で役立ちます。 

figure(3)
f=semilogy(resvecCG_2./resvecCG_2(1));
hold on
semilogy(f.Parent.XLim,[1e-6 1e-6],"--")
semilogy([20 20], f.Parent.YLim,"--")
semilogy(f.Parent.XLim,[1e-12 1e-12],"--")
title("Evolution of Relative Residual");
ylabel("Relative Residual");
xlabel("Iteration Step");
legend("Residuals of CG","Default Tolerance", ...
    "Default Number of Steps","Custom Tolerance")
hold off

明らかに、この方程式の良好な解を得るには既定のステップ数が不十分です。

前処理付き共役勾配法による連立線形方程式の求解

前処理付き共役勾配法 (PCG) 法により、方程式の求解の効率を向上できます。まず、前処理行列 M を使用して、連立線形方程式の前処理を行います。次に、CG 法を使用して、前処理された方程式を求解します。PCG 法の反復回数は CG 法よりかなり少なくなります。

PCG 法では、関数 pcg も使用されます。追加の入力として、適切な前処理行列 M を指定できます。

理想的な前処理行列は、逆行列 M-1 が係数行列の逆行列 A-1 に近い近似になる一方、計算がより容易な行列です。この例では、A の対角を使用して、連立線形方程式の前処理を行います。 

M = spdiags(spdiags(A,0),0,N,N);
tPCG = tic;
    [xPCG,flagPCG,relresPCG,iterPCG,resvecPCG]=pcg(A,b,tolerance,maxit,M);
tPCG = toc(tPCG);

figure(4)
hold off;
semilogy(resvecCG_2./resvecCG_2(1))
hold on;
semilogy(resvecPCG./resvecPCG(1))
title("Evolution of Relative Residual");
ylabel("Relative Residual");
xlabel("Iteration Step");
legend("Residuals of CG","Residuals of PCG with M \approx diag(A)")

前の図は、前処理をしない方程式と比較して、PCG 法では収束に要するステップ数が著しく少ないことを示しています。この結果は、実行回数にも反映されます。

fprintf("\nTime to solve system with CG:  %d s" + ...
    " \nTime to solve system with PCG: %d s",tCG,tPCG);
Time to solve system with CG:  8.819875e+00 s 
Time to solve system with PCG: 1.307804e+00 s

PCG 法は、この例の方程式をより少ない反復ステップ数で解くだけでなく、より正確な解を返します。

errPCG = abs(xExact-xPCG);
figure(5)
hold off
semilogy(errCG_1,"o");
hold on
semilogy(errCG_2,"d");
semilogy(errPCG,"x");
title("Comparison of absolute error");
ylabel("Absolute error");
xlabel("Element in x");
legend("CG default","CG custom","PCG");

計算が完了したら、並列プールを削除できます。関数gcpは現在の並列プール オブジェクトを返すため、現在のプールを削除できます。

delete(gcp("nocreate"))
Parallel pool using the 'Processes' profile is shutting down.

この例で使用した Wathen 行列は、適切な前処理行列によって求解の効率を著しく向上できることをよく示しています。Wathen 行列の非対角要素は比較的小さいため、M=diag(A) を選択すると適切な前処理行列が得られます。任意の行列 A の場合、前処理行列を求めることはそう単純ではありません。

線形システムにより微分方程式を近似し、マルチグリッド前処理行列をもつ分散反復ソルバーを使用してこれを解く方法の例は、分散離散化でマルチグリッド前処理行列を使用して微分方程式を解くを参照してください。

参考

| |

トピック