このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
シストリック アレイでの CORDIC を使用したハードウェア効率に優れた QR 分解の実装
この例では、Simulink® でハードウェア効率に優れた MATLAB® コードを使用して行列の QR 分解を計算する方法を示します。
方程式系を解く、つまり QR 分解を使用して行列方程式 AX = B
の最小二乗解を計算するには、R
と Q'B
を計算します。ここで QR = A
および RX = Q'B
です。R
は上三角行列、Q
は直交行列です。Q
と R
だけが必要な場合は、B
を単位行列に設定します。
この例では、CORDIC アルゴリズムを使用してギブンス回転を適用して、行列 A
から R
が計算されています。C = Q'B
は、同じギブンス回転を適用して行列 B
から計算されています。このアルゴリズムでは反復的なシフトと加算のみを使用してこれらの計算を実行します。
モデルを開く
fxpdemo_real_4x4_systolic_array_QR_model
モデルを開きます。この例では、A
の行数と列数が 4
でなければなりません。行列 B
には 4
つの行がなければなりませんが、列数は任意です。シストリック アレイは、同じパターンに従ってどのサイズにも拡張できます。このモデルでは、固定小数点、倍精度、および単精度のデータ型がサポートされています。この実装は、実数入力行列でのみ動作します。
open_system('fxpdemo_real_4x4_systolic_array_QR_model')
独自の入力行列 A
と B
を入力するには、対応する Enabled Subsystem ブロックのブロック パラメーターをモデルの左側に開きます。シミュレーション後、モデルは計算された出力行列 C = Q'B
および R
をワークスペースに返します。"transpose(Q)*B, R 4x4 Real CORDIC Systolic Array" サブシステムのブロック パラメーターで CORDIC 反復数を指定できます。入力が固定小数点の場合、CORDIC 反復数は語長より小さくなければなりません。計算の精度は反復ごとに 1 ビットずつ、入力の語長まで向上します。
アルゴリズムでの因数分解の実行方法を確認するには、"transpose(Q)*B, R 4x4 Real CORDIC Systolic Array" サブシステムのマスク内を確認します。注釈には、下対角要素をゼロ設定するために演算されている行列の行が示されます。
CORDIC を使用してギブンス回転を実行する MATLAB Function ブロックの MATLAB コードを確認するには、引き続きブロック マスク内を確認します。
QR を使用して行列方程式 Ax = B
を解く
行列方程式 AX = B
を解く最初のステップは RX = Q'B
を計算することです。ここで、R は上三角行列、Q
は直交行列で Q*R = A
です。
以下の入力は倍精度浮動小数点型です。反復数は double の仮数のビット数である 52
に設定します。
NumberOfCORDICIterations = 52; A = 2*rand(4,4)-1; B = 2*rand(4,4)-1;
R
と C = Q'B
を計算するには、モデルをシミュレートします。
sim fxpdemo_real_4x4_systolic_array_QR_model
R
R = 4×4
1.5149 -0.0519 1.7292 -0.3224
0 0.9593 -0.0259 -0.0879
0 0 0.2565 1.0888
0 0 0 -0.6429
C
C = 4×4
0.5942 -0.2382 0.0676 -0.9370
-0.8887 0.6146 -0.5758 0.3051
0.1725 0.7339 0.5409 0.5374
0.8540 1.1078 -0.2183 -0.5620
R
および C = Q'B
で後退代入すると MATLAB と同じ結果になることを確認します。
X = R\C
X = 4×4
-7.1245 -12.1131 -0.6637 1.4236
-0.8779 0.7572 -0.5511 0.3545
6.3113 10.1759 0.6673 -1.6155
-1.3283 -1.7231 0.3396 0.8741
X_should_be = A\B
X_should_be = 4×4
-7.1245 -12.1131 -0.6637 1.4236
-0.8779 0.7572 -0.5511 0.3545
6.3113 10.1759 0.6673 -1.6155
-1.3283 -1.7231 0.3396 0.8741
組み込みの MATLAB の解と CORDIC QR の解の差のノルムは小さくなります。
norm(X - X_should_be)
ans = 3.2578e-14
Q
と R
の計算
Q
および R
を計算するには、B
を単位行列と等価に設定します。
NumberOfCORDICIterations = 52; A = 2*rand(4,4)-1; B = eye(size(A,1),'like',A); sim fxpdemo_real_4x4_systolic_array_QR_model
B が単位行列と等しい場合、Q = C' になります。
Q = C';
理論的な QR 分解は QR==A
であるため、計算された QR と A
の差は小さくなります。
norm(Q*R - A)
ans = 2.2861e-15
QR は一意ではない
QR 分解が一意であるのは R
の行と Q
の列の符号までです。R
の対角要素をすべて正にすると、一意の QR 分解にすることができます。
D = diag(sign(diag(R))); Qunique = Q*D
Qunique = 4×4
-0.3086 0.1224 -0.1033 -0.9376
-0.6277 -0.7636 -0.0952 0.1174
-0.5573 0.3930 0.7146 0.1559
0.4474 -0.4975 0.6852 -0.2877
Runique = D*R
Runique = 4×4
1.4459 -0.8090 0.1547 0.3977
0 1.1441 0.0809 -0.2494
0 0 0.8193 0.1894
0 0 0 0.4836
モデルから計算された QR と MATLAB 関数 qr
を比較します。
[Q0,R0] = qr(A); D0 = diag(sign(diag(R0))); Q0 = Q0*D0
Q0 = 4×4
-0.3086 0.1224 -0.1033 -0.9376
-0.6277 -0.7636 -0.0952 0.1174
-0.5573 0.3930 0.7146 0.1559
0.4474 -0.4975 0.6852 -0.2877
R0 = D0*R0
R0 = 4×4
1.4459 -0.8090 0.1547 0.3977
0 1.1441 0.0809 -0.2494
0 0 0.8193 0.1894
0 0 0 0.4836
ハードウェア効率の高い実装のための固定小数点の使用
固定小数点入力データ型を使用して、ASIC や FPGA デバイス向けの効率的な HDL コードを生成します。
A
および B
の 3 次元配列を作成して、モデルから多くのテスト入力を実行できます。
n_test_inputs = 100;
行列 A
と B
の乱数入力を -1
~ 1
でスケーリングします。固定小数点の語長を 18 ビット、小数部の長さを 14 ビットに設定して QR 分解の成長と CORDIC アルゴリズムの中間計算を可能にします。
word_length = 18; fraction_length = 14;
CORDIC 反復の最高精度数は語長から 1 を引いた数です。CORDIC 反復数が word_length - 1
より小さい値に設定されると、次の Ready 信号への遅延と時間刻みが短くなりますが、精度は低くなります。生成されたコードは固定小数点の語長よりも大きいシフトをサポートしていないため、CORDIC 反復数はこれより大きい値に設定できません。
NumberOfCORDICIterations = word_length - 1
NumberOfCORDICIterations = 17
乱数テスト入力は連結されて、時間 k
における入力は A(:,:,k)
および B(:,:,k)
になります。A
および B
の各要素は、-1
~ +1
の一様な確率変数になります。
A = 2*rand(4,4,n_test_inputs)-1;
Q=C'
となるように単位行列として B
を選択します。
B = eye(4); B = repmat(B,1,1,n_test_inputs);
A
を固定小数点にキャストし、B
を A
と同様にキャストします。
A = fi(A,1,word_length,fraction_length);
B = cast(B,'like',A);
モデルのシミュレーション
固定小数点の入力を使用してモデルをシミュレートします。
sim fxpdemo_real_4x4_systolic_array_QR_model
誤差の計算とプロット
誤差を計算します。
norm_error = zeros(1,size(R,3)); for k = 1:size(R,3) Q_times_R_minus_A = double(C(:,:,k))'*double(R(:,:,k)) - double(A(:,:,k)); norm_error(k) = norm(Q_times_R_minus_A); end
誤差をプロットします。誤差は 10^-3
と同程度になります。
clf plot(norm_error,'o-') grid on title('norm(QR - A)')
%#ok<*NASGU,*NOPTS>