最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

ハードウェア効率に優れた Real Burst Matrix Solve の実装

この例では、Simulink® モデルに組み込まれている、ハードウェア効率に優れた MATLAB® コードを使用して方程式系を解く方法を示します。この例で使用されるモデルは、固定小数点入力の HDL コード生成に適しています。アルゴリズムでは、ソルバー アルゴリズムのさまざまな手順間で計算単位およびメモリ単位を共有する計算アーキテクチャが使用されます。これは、固定小数点アルゴリズムをリソースに制約がある FPGA デバイスや ASIC デバイスに展開する際に便利です。このソルバーのスループットが完全にパイプライン化された実装より小さい一方で、オンチップ フットプリントも小さくなるため、リソースに配慮した設計に適しています。

方程式系の求解

線形方程式系は、次のような行列形式で表現できます。

$$A x = b,$$

ここで $A$ は既知のデータの $m$$n$ 列の行列で、$x$ は未知数の $n$ 行 1 列のベクトルで、$b$$m$ 行 1 列のベクトルです。

一般に、線形システムは、最初にシステムを下三角形式に分解し、次に後退代入により未知数を検索することで最も効率的に求解します。この方法は、数値的に安定していて、かつ効率的です。したがって、これは組み込みアプリケーションに適しています。

任意の行列 $A$ は、次のように分解できます。

$$A = Q R,$$

ここで $Q'Q = I$$R$$m$$n$ 列の上三角行列であり、$I$$m$$m$ 列の単位行列です。これにより、方程式系を上三角形式に変換できるようになります。

$$R = Q'b.$$

この方程式は、$A$ が非特異である場合、後退代入により解くことができます。この方法の利点は、$Q$ を計算する必要がないことです。代わりに、$A$$R$ への変換は $b$ に同時に適用できます。これにより、小さなメモリと計算フットプリントを使用する前の方程式と同じ形式が得られます。

CORDIC QR 分解アルゴリズム

QR 分解は一連のギブンス回転により実行されます。これらは Coordinate Rotation Digital Computer (CORDIC) 手法を使用してハードウェアに実装されます。CORDIC QR アルゴリズムの詳細な説明については、CORDIC を使用した QR 分解の実行の例を参照してください。

ハードウェア効率に優れた固定小数点行列の計算

Real Burst Matrix Solve using QR Decomposition は、2 進小数点でスケーリングされた固定小数点データの HDL コード生成をサポートします。これはこの適用例に留意して設計されていて、ハードウェア固有のセマンティクスと最適化を採用します。これらの最適化の 1 つにリソース共有があります。

複雑なアルゴリズムを FPGA デバイスや ASIC デバイスに展開する場合、計算のリソースの使用率とその合計スループットの間にトレードオフがある場合が多いです。完全にパイプライン化され、並列化されたアルゴリズムのスループットは最大となる一方で、実際のデバイスで展開するには過度にリソース集約的である場合があります。1 つ以上のコア計算回路にスケジューリング ロジックを実装することで、計算全体でリソースの再利用が可能になります。結果として、合計スループットのコストにもかかわらず、実装ではフットプリントがより小さくなります。リソースを共有する設計は全体的なレイテンシ要件を満たすことができるため、これはたいていの場合、許容可能なトレードオフです。

Real Burst Matrix Solve using QR Decomposition の主要な計算単位のすべては、計算ライフ サイクル全体で再利用されます。これには、QR 分解中のギブンス回転の実行に使用される CORDIC 回路だけでなく、後退代入中に使用される加算器と乗算器も含まれます。これにより、FPGA デバイスや ASIC デバイスへの展開時に DSP リソースとファブリック リソースの両方を節約できます。

サポートされるデータ型

シミュレーションでは、単精度、倍精度、2 進小数点でスケーリングされた固定小数点、2 進小数点でスケーリングされた倍精度データがサポートされます。ただし、HDL コード生成の場合、2 進小数点でスケーリングされた固定小数点型のみがサポートされます。

Real Burst Matrix Solve のパラメーター

Real Burst Matrix Solve using QR Decomposition は、行列 A および B をモデルまたはベース ワークスペースのいずれかから使用して、解を求める方程式系を定義します。そのマスクでは、行列 A および B の行数、行列 A の列数、行列 B の列数、および出力データ型の 4 つの入力パラメーターを取ります。ブロックのマスクを以下に示します。

これらのパラメーターは、それぞれ変数 mnp および OutputType によってワークスペースで定義されます。さらに、変数 num_samples をモデルまたはベース ワークスペースで定義する必要があります。

モデル データの設定

最初に、解を求める一連の方程式系が必要です。実際のデータはガウス アンサンブルから得られることが多いため、この例では、標準正規分布から得られる乱数行列が使用されます。この標準正規分布の平均値は 0 であり、分散は 1 です。このプロパティは、自然にみられるランダム性のソースの多くについて、すべてではありませんが説明します。これにより、Real Burst Matrix Solve using QR Decomposition の使用方法を説明するために、この分布から得られるデータを例として使用します。

最初に、解を求める系のサイズを選択します。

m = 4;
n = 4;
p = 1;
num_samples = 100;

次に、倍精度のランダム データを生成します。このデータは、必要なダイナミック レンジの数値ベンチマークとして使用されます。これは、固定小数点への変換時に役立ちます。

A = randn(m,n,num_samples);
B = randn(m,p,num_samples);

ここで、最初の 2 つのインデックスは、各サンプルの行列のインデックスに該当しますが、3 番目のインデックスはサンプル番号に該当します。たとえば、A(1,2,3) は、3 番目のサンプルの 1 行 2 列目の要素に該当します。

浮動小数点演算のためのモデル例の生成

固定小数点に変換する前に、浮動小数点演算系の解を求めます。次に、モデルの固定小数点バージョンの適切な入力データ型と出力データ型を選択するために、このシミュレーションから取得した範囲情報を使用します。

モデルの浮動小数点バージョンを生成するには、次の関数呼び出しを使用します。

mdl = fixed.getMatrixSolveModel(A,B);

この関数では、モデル ワークスペースにデータを設定します。したがって、$A$ または $B$ に格納されているデータを変更する場合、モデル ワークスペース内のデータを変更するか、固定小数点データを使用して新しいモデルを生成します。この例の後の方で、後者の方法を使用します。

ソースから Matrix Solve ブロックへのデータの受け渡し

上記の関数によって生成されたモデルは、修正を行わずにシミュレーションを実行するように既に設定されています。生成されたモデルは、データ ハンドラーと呼ばれるブロックを使用してデータを行列ソルバーに送信します。ソルバーの ready 端子をリッスンし、ready の立ち上がりエッジで validIn を High に引き上げることで行列ソルバーと通信を行います。これにより、端子 A(i,:) と B(i,:) の現在の行が有効であることを Real Burst Matrix Solve using QR Decomposition に通知します。ready および valid が同時に High になると、これらの端子のデータは、行列の次の有効な行として扱われます。

ready および validIn の両方が true であるとソルバーが認識するたびに、ソルバーは、端子 A(i,:) と B(i,:) のデータ入力を消費します。そのようにする場合は、validIn を High のままにしておかなければなりません。次に、Real Burst Matrix Solve using QR Decomposition は、処理の準備ができるまで、この行をバッファー内に保持します。この演算のモードを以下に示します。ここでは、2 つのサンプルが消費されています。

ready が High になると、2 つの行を消費できます。2 つの行が送信されると、Real Burst Matrix Solve using QR Decomposition は 2 番目の行を内部バッファーに格納します。ready がアサートされる場合、2 つの行の送信は必須ではなく、validIn がデアサートされるために 1 つのサンプルのみが送信される場合、ready はデアサートされ、計算のために追加のデータが必要な場合にのみ再アサートされます。

ready からデータ ハンドラーへのパスか、データ ハンドラーと validIn の間のパスのいずれかで遅延がさらに追加されると、予期しないバック プレッシャーが追加されます。このシナリオでは、サンプルは除外されます。ready の立ち上がりエッジのデータのみを送信するか、ブロックとデータ ソース間に追加の遅延を挿入しないことを推奨します。さらに、このブロックを使用した HDL コードの生成を計画している場合、追加の遅延がこれらのパスに追加される可能性があるためパイプライン方式の最適化を無効にすることを推奨します。

モデルのシミュレーション

モデルのシミュレーションを実行するには、[シミュレーション] ボタンをクリックするか、sim コマンドを使用します。

out = sim(mdl);

出力の処理

Real Burst Matrix Solve using QR Decomposition は一度に 1 行ずつデータを出力します。ブロックが結果行を出力するたびに、validOut をアサートします。解は、変数 out のフィールド X_out の出力です。このデータは 1 × p × n*num_samples の行の配列として格納されたままになります。MATLAB® で解の精度を評価する前に、それを n x p x num_samples の配列に変換すると便利です。これは、ユーティリティ関数 matrixSolveOutputToMatrix を使用して行うことができます。

X = matrixSolveOutputToMatrix(out.X_out,n,p,num_samples);

MATLAB® での予期される解の計算

Simulink® モデルで計算される解の精度は、MATLAB® で容易に計算できます。理論では、系の解の相対誤差は、行列の条件数が増えるほど増加します。これを確認するために、

$$||e||_{2} = \frac{||Ax - b||_{2}}{||b||_{2}}$$

で与えられる誤差を $A$ の条件数に対してプロットします。

condNumber = zeros(1,num_samples);
relativeError = zeros(1,num_samples);
for idx = 1:num_samples
   condNumber(idx) = cond(A(:,:,idx));
   relativeError(idx) = norm(B(:,:,idx) -  A(:,:,idx)*X(:,:,idx))./norm(B(:,:,idx));
end

出力のプロット

次に、誤差計算の出力をプロットします。条件数と誤差の両方が数桁にわたるため、対数スケールが両方の座標で使用されることに注意してください。ここでは 2 つの値の間の正の相関が非常に明確です。

figure(1);
clf;
h1 = axes();
hold on;
plot(condNumber,relativeError,'bo');
h1.XScale = 'log';
h1.YScale = 'log';
xlabel('cond(A)');
ylabel('Relative Error');

モデルの固定小数点バージョンの作成

上記でシミュレーションを行った倍精度アルゴリズムは、非常に正確です。設計の固定小数点バージョンの数値精度を高く維持するには、QR 分解手順または後退代入手順のいずれかで発生する大きさの増加が原因で起こる可能性があるオーバーフローの回避に十分な大きさの入力データ型と出力データ型が必要です。

入力データ型の選択

Real Burst QR Decomposition の場合、R(i,:) の型は A(i,:) の型と同じです。同様に、C(i,:) の型は B(i,:) の型と同じです。これらの型は、計算全体で使用されます。したがって、入力データ型は、R(i,:) および C(i,:) それぞれの計算で必要な取り得る最大値に対応する必要があります。

計算中の入力データの増加を引き起こす効果は 2 つあります。それは、QR 分解で使用される回転の性質を保持するユークリッド ノルムに起因する A の要素の増加と、これらの回転の CORDIC 実装の内部的な拡大です。これらの 2 つの効果は、出力の要素の大きさを次のように制限します。

ceil(log2(1.65*sqrt(m)*max(abs(A(:)))))

(R(i,:) の場合)

ceil(log2(1.65*sqrt(m)*max(abs(B(:)))))

(C(i,:) の場合)この例では、A および B は同じ統計分布から得られたため、同じ語長と小数部の長さの両方を設定します。一般に、A および B の語長は同じであるため、それらの小数部の長さは同じである必要はありません。

maximumAbsValueInDataset = max(max(abs(A(:))),max(abs(B(:))))
inputIntegerLength = ceil(log2(1.65*sqrt(m)*maximumAbsValueInDataset))
maximumAbsValueInDataset =

    3.5784


inputIntegerLength =

     4

入力と同じ 18 ビットの語長型が使用されます。これらを変更して、異なる用途に適するようにします。以下で計算される変数 inputType は、NumericType オブジェクトの型情報を符号化します。

inputWordLength = 18
inputFractionLength = 18 - inputIntegerLength - 1
inputType = numerictype(1, inputWordLength, inputFractionLength)
inputWordLength =

    18


inputFractionLength =

    13


inputType =


          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 18
        FractionLength: 13

CORDIC QR 分解アルゴリズムでの増加の上限の詳細については、Prevent Overflow in Fixed Point Rを参照してください。

後退代入に必要なダイナミック レンジの計算

データは、小さい値による除算が原因のアルゴリズムの後退代入の一部でも大きくなる場合があります。関連する入力データを使用したシミュレーションにより、適切な出力データ型を選択することができます。この手順では、X(i,:) の出力値の大きさは以下により制限されています。

1.65*sqrt(m)*max(abs(B(:)))/sigma_min,

ここで sigma_min はサンプル行列のすべてにわたり観測される最も小さい特異値です。これらの値は次のように計算されます。

sigma = zeros(1, num_samples);
for idx = 1:num_samples
    sigma(idx) = min(svd(A(:,:,idx)));
end
sigma_min = min(sigma);
backsubstituteWordlengthGrowth = ceil(log2(1.65*sqrt(m)/sigma_min));

これにより、以下に示すように出力データ型が得られます。

outputWordlength = inputWordLength + backsubstituteWordlengthGrowth
outputFractionLength = inputFractionLength
OutputType = numerictype(1, outputWordlength, outputFractionLength)
outputWordlength =

    30


outputFractionLength =

    13


OutputType =


          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 30
        FractionLength: 13

モデルの固定小数点バージョンの作成

固定小数点のシミュレーションを実行するには、まず入力データを計算される入力型に変換します。

A_Fixed = fi(A, inputType);
B_Fixed = fi(B, inputType);

モデル ワークスペースの変数 ABOutputType を更新することで、既存のモデルを更新できます。ただし、モデル作成関数にも指定した出力型を使用してバージョンを作成するオプションがあります。

mdl_fi = fixed.getMatrixSolveModel(A_Fixed, B_Fixed, OutputType);

固定小数点のモデルのシミュレーション

out_fi = sim(mdl_fi);

固定小数点出力の MATLAB® 配列での保存

X_fi = matrixSolveOutputToMatrix(out_fi.X_out,n,p,num_samples);

MATLAB® での予期される解の計算

入力計算からの入力行列の条件数は既知であるため、以下の値を使用して固定小数点の解の相対誤差を計算できます。

relativeErrorFixedPoint = zeros(1, num_samples);
for idx = 1:num_samples
   relativeErrorFixedPoint(idx) = norm(double(B_Fixed(:,:,idx) -  A_Fixed(:,:,idx)*X_fi(:,:,idx)))./norm(double(B_Fixed(:,:,idx)));
end

出力のプロット

計算の誤差を再度プロットしますが、今回は固定小数点の結果を使用します。量子化のため、誤差はさらに大きいことに注意してください。

figure(2);
clf;
h2 = axes();
hold on;
plot(condNumber,relativeErrorFixedPoint,'bo');
h2.XScale = 'log';
h2.YScale = 'log';
xlabel('cond(A)');
ylabel('Relative Error');