Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

分散配列を使用した直接法による線形方程式系の求解

この例では、分散配列を使用して、直接法によって Ax=b の形式の線形方程式系を並列で求解する方法を説明します。

分散配列は大規模な数学的計算、たとえば大規模な線形代数の問題などによく適しています。クライアントのメモリに保存された配列と同様に、分散配列を使用して定義された線形方程式系はmldivideを使用して解くことができるため、コードを変更する必要はありません。

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

直接法による線形方程式の求解では通常、係数行列 A を分解して解を計算します。mldivide は、A の構造と、A が非スパースまたはスパースのいずれであるかに応じて、一連の直接ソルバー メソッドから 1 つを選択します。

この例では、既知の厳密解 x を使用して、Ax=b の形式の単純な線形方程式系を解く方法を示します。この方程式系は、行列 A および列ベクトル b によって定義されます。解 x も列ベクトルです。この例の方程式系は、非スパース行列とスパース行列を使用して定義されています。方程式系の定義に分散配列またはクライアントのメモリ上の配列のいずれを使用しても、必要なコードは同じです。

関連する、反復ソルバーと分散配列の使用方法を示した例は、分散配列を使用した反復法による線形方程式系の求解を参照してください。

非スパース行列系の求解

まず、係数行列 A をクライアントのメモリ内の変数 A として定義してから、この行列を関数作成に渡して、同じ行列の分散バージョン ADist を作成します。関数 distributed を使用すると、MATLAB は既定のクラスター設定を使用して、並列プールを自動的に起動します。

n = 1e3;
A = randi(100,n,n);
ADist = distributed(A);

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

b = sum(A,2);
bDist = sum(ADist,2);

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

xEx = ones(n,1);
xDistEx = ones(n,1,'distributed');

これで線形方程式系が定義されたため、mldivide を使用して直接、方程式系を求解できます。MATLAB では、特殊演算子 \ を使用して mldivide を呼び出すことができます。mldivide は分散配列を自動的にサポートするため、分散された方程式系を解くためにコードを変更する必要はありません。

解を計算した後、得られた結果 x の各要素と xexact の各予測値との誤差をチェックできます。

x = A\b;
err = abs(xEx-x);

xDist = ADist\bDist;
errDist = abs(xDistEx-xDist);

figure
subplot(2,1,1)
semilogy(err,'o');
title('System of Linear Equations with Full Matrices');
ylabel('Absolute Error');
xlabel('Element in x');
ylim([10e-17,10e-13])
subplot(2,1,2)
semilogy(errDist,'o');
title('System of Linear Equations with Distributed Full Matrices');
ylabel('Absolute Error');
xlabel('Element in x');
ylim([10e-17,10e-13])

分散配列とクライアントに保存された配列のいずれでも、x の計算結果と厳密な結果 xexact との絶対誤差はわずかです。いずれの配列タイプでも、解の精度はほぼ同じです。

mean(err)
ans = 1.6031e-13
mean(errDist)
ans =

   1.2426e-13

スパース行列系の求解

分散配列はスパース データも格納できます。係数行列 A を作成するために、sprandspeye を使用して、スパース乱数行列とスパース単位行列を直接生成します。単位行列を加算することで、A が特異行列またはほぼ特異行列 (いずれも分解が困難) として作成されることを回避できます。

n = 1e3;
density = 0.2;
A = sprand(n,n,density) + speye(n);
ADist = distributed(A);

右辺のベクトル bA の行の和として選択した場合、非スパース行列系の解と同じ形式の厳密解が得られます。

b = sum(A,2);
bDist = sum(ADist,2);
xEx = ones(n,1);
xDistEx = ones(n,1,'distributed');

非スパース行列の場合と同様に、これで mldivide を使用して直接この線形方程式系を解き、得られた結果と予測値との誤差をチェックできます。

x = A\b;
err = abs(xEx-x);

xDist = ADist\bDist;
errDist = abs(xDistEx-xDist);

figure
subplot(2,1,1)
semilogy(err,'o');
title('System of Linear Equations with In-Client Sparse Matrices');
ylabel('Absolute Error');
xlabel('Element in x');
ylim([10e-17,10e-13])
subplot(2,1,2)
semilogy(errDist,'o');
title('System of Linear Equations with Distributed Sparse Matrices');
ylabel('Absolute Error');
xlabel('Element in x');
ylim([10e-17,10e-13])

非スパース行列系の場合と同様に、クライアント上の配列と分散配列のいずれを使用して線形方程式系を解いても、ほぼ同じ精度で解が出力されます。

mean(err)
ans = 1.6031e-13
mean(errDist)
ans =

   1.2426e-13

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

delete(gcp('nocreate'));

解の精度の改善

A が大規模でスパースな特定タイプの係数行列である場合、直接の分解よりも効率的に方程式系を解く方法があります。このような場合、反復法で線形方程式系を解く方がより効率的である可能性があります。反復法では一連の近似解が生成されて最終結果に収束します。反復法を使用し、大規模でスパースな入力行列を指定して線形方程式を解く方法の例は、分散配列を使用した反復法による線形方程式系の求解を参照してください。

参考

| |

関連するトピック