Main Content

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

CORDIC を使用した QR 分解の実行

この例で使用するアルゴリズムは、CORDIC (Coordinate Rotation Digital Computer) を使用して実装された QR 分解です。この例では、浮動小数点および固定小数点の両方のデータ型で機能する MATLAB® コードの記述方法を説明します。

固定小数点ターゲットを対象とするアルゴリズムを記述するための優れた方法の 1 つは、アルゴリズムが機能することを検証できるよう MATLAB で組み込み浮動小数点型を使用して記述することです。アルゴリズムが固定小数点型に対して機能するよう改良する場合、最善の方法は同じコードが浮動小数点でも引き続き機能するように記述することです。こうすることにより、デバッグ時に入力を浮動小数点型と固定小数点型との間で切り替えて、動作の違いの原因がオーバーフローや量子化などの固定小数点型の影響なのか、それともアルゴリズムの違いなのかを判断することができます。アルゴリズムが浮動小数点ターゲットに適していない場合 (以下の例で示すような CORDIC を使用するケースなど) でも、デバッグ目的で MATLAB コードを浮動小数点型で機能させることにはメリットがあります。

これとは対照的に、ターゲットが浮動小数点型の場合、戦略はまったく異なります。たとえば、QR アルゴリズムは、ハウスホルダー変換および行ピボットまたは列ピボット付きの浮動小数点型で実行されることがよくあります。しかし、固定小数点型では、CORDIC を使用してピボットなしのギブンス回転を適用する方が効率的であることがよくあります。

この例では、最初のケースについて説明します。ターゲットは固定小数点型ですが、作成とデバッグを簡単に行えるようにデータ型に依存しないアルゴリズムを記述します。

この例では、複数のシステムに適用できるさまざまなコーディング方法を説明します。この例で使用する重要な設計パターンは、以下のとおりです。

  • データ型からの独立性: MATLAB コードがデータ型に依存せずに、fidouble、およびsingleのいずれに対しても同じように機能するようなアルゴリズムを記述します。

  • オーバーフロー防止: オーバーフローを確実に防止する方法。固定小数点型でオーバーフローを回避する方法を説明します。

  • 方程式の求解: 計算の効率性を活用する方法。定義対象を切り分けて、コード範囲を絞り込んでください。

この例の主要な部分は、ギブンス回転に CORDIC を使用する固定小数点演算にqr分解を実装することです。MATLAB コードがデータ型に依存せずに、固定小数点型、倍精度浮動小数点型、単精度浮動小数点型のいずれにも同様に機能するようなアルゴリズムを記述します。

M 行 N 列の行列 A の QR 分解によって、M 行 N 列の上三角行列 R および A = Q*R であるような M 行 M 列の直交行列 Q が生成されます。対角より下がすべてゼロであれば、行列は上三角行列になります。Q'*Q = eye(M)、つまり単位行列である場合、M 行 M 列の行列 Q は直交行列です。

QR 分解は最小二乗問題で広く使用されます。たとえば、適応フィルターでは再帰的最小二乗 (RLS) アルゴリズムが使用されます。

CORDIC アルゴリズムは、固定小数点型で QR アルゴリズムを計算する場合に便利です。シフト演算と加算演算だけを使用して、CORDIC で直交ギブンス回転を適用できるためです。

CORDIC QR アルゴリズムの定義

以下の関数 cordicqr は CORDIC QR アルゴリズムを実装します。ここで、A は M 行 N 列の実数行列、niter は CORDIC の反復回数です。出力 Q は M 行 M 列の直交行列であり、RQ*R = A である M 行 N 列の上三角行列です。

function [Q,R] = cordicqr(A,niter)
  Kn = inverse_cordic_growth_constant(niter);
  [m,n] = size(A);
  R = A;
  Q = coder.nullcopy(repmat(A(:,1),1,m)); % Declare type and size of Q
  Q(:) = eye(m);                          % Initialize Q
  for j=1:n
    for i=j+1:m
      [R(j,j:end),R(i,j:end),Q(:,j),Q(:,i)] = ...
          cordicgivens(R(j,j:end),R(i,j:end),Q(:,j),Q(:,i),niter,Kn);
    end
  end
end

この関数は、データ型から独立するよう記述されました。これは、double および single 浮動小数点型でも固定小数点型 fi オブジェクトでも、同様に機能します。

データ型から独立しているコードを記述することについて、最も面倒な点の 1 つは、新しい変数のデータ型とサイズを指定することです。データ型を明示的に指定することなく保持するため、出力 R は入力 A と同じになるように設定されています。

    R = A;

この関数はデータ型に依存しないようになっていますが、MATLAB Coder™ で効率的な C コードを生成できるようにも記述されています。MATLAB では、ほとんどの場合、変数の宣言および開始値の設定を以下のように 1 ステップで行います。

    Q = eye(m)

ただし、Q = eye(m)Q を常に倍精度浮動小数点変数として生成します。A が固定小数点の場合は Q は固定小数点でなければならず、A が single の場合は Q は single でなければならないといったようになります。

したがって、Q の型とサイズを 1 ステップで宣言してから、次のステップで開始値を設定する必要があります。これにより、MATLAB Coder には正しい型とサイズで効率的な C プログラムを作成するために必要な情報が提供されます。完成したコードで、M 行 M 列の単位行列であり、データ型が A と同じになるように出力 Q の開始値を設定します。

    Q = coder.nullcopy(repmat(A(:,1),1,m)); % Declare type and size of Q
    Q(:) = eye(m);                          % Initialize Q

関数coder.nullcopyは、開始値を設定することなく Q のサイズと型を宣言します。A の最初の列をrepmatで拡張したものは、MATLAB によって生成されたコードには現れません。サイズの指定に使用されるだけです。A(:,1:m) の代わりに関数 repmat が使用されました。これは、A では列数より行数の方が多いからです。これは、最小二乗問題ではよくあることです。coder.nullcopy で宣言するときは、値を配列の各要素に必ず代入してください。そうしないと、メモリが初期化されないからです。

この代入パターンは何度も繰り返し出てきます。これは、コードをデータ型から独立させるためのもう 1 つのキー ポイントです。

この関数の核心は、所定の位置の直交ギブンス回転を R の行に適用して下対角要素をゼロに設定し、こうして上三角行列を形成することです。この同じ回転は、単位行列の列にも適用され、直交行列 Q が形成されます。ギブンス回転は関数 cordicgivens を使用して適用されます。これについては、次のセクションで定義されます。R の行数と Q の列数は、関数 cordicgivens の入力および出力として使用されます。その結果、計算は正しい位置で実行され、RQ が上書きされます。

[R(j,j:end),R(i,j:end),Q(:,j),Q(:,i)] = ...
    cordicgivens(R(j,j:end),R(i,j:end),Q(:,j),Q(:,i),niter,Kn);

CORDIC ギブンス回転の定義

関数 cordicgivens は、x(1) = R(j,j)y(1) = R(i,j) によって定義される角度について CORDIC 反復を実行して、ギブンス回転を x = R(j,j:end) 行と y = R(i,j:end) 行に適用します。ただし、i > j であり、こうして R(i,j) がゼロに設定されます。同じ回転が、u = Q(:,j) 列と v = Q(:,i) 列にも適用されます。こうして、直交行列 Q が形成されます。

function [x,y,u,v] = cordicgivens(x,y,u,v,niter,Kn)
  if x(1)<0
    % Compensation for 3rd and 4th quadrants
    x(:) = -x;  u(:) = -u;
    y(:) = -y;  v(:) = -v;
  end
  for i=0:niter-1
    x0 = x;
    u0 = u;
    if y(1)<0
      % Counter-clockwise rotation
      % x and y form R,         u and v form Q
      x(:) = x - bitsra(y, i);  u(:) = u - bitsra(v, i);
      y(:) = y + bitsra(x0,i);  v(:) = v + bitsra(u0,i);
    else
      % Clockwise rotation
      % x and y form R,         u and v form Q
      x(:) = x + bitsra(y, i);  u(:) = u + bitsra(v, i);
      y(:) = y - bitsra(x0,i);  v(:) = v - bitsra(u0,i);
    end
  end
  % Set y(1) to exactly zero so R will be upper triangular without round off
  % showing up in the lower triangle.
  y(1) = 0;
  % Normalize the CORDIC gain
  x(:) = Kn * x;  u(:) = Kn * u;
  y(:) = Kn * y;  v(:) = Kn * v;
end

固定小数点型で CORDIC を使用することには、標準ギブンス回転を使用することに比べて、CORDIC では平方根演算も除算演算も使用されないというメリットがあります。主要ループで必要とされるのは、ビット シフト、加算、および減算だけであり、スカラー量とベクトル量を最後に一度乗算するだけで、CORDIC ゲインを正規化することができます。また、CORDIC 反復はパイプライン方式のアーキテクチャでも機能します。

各反復においてビット シフトは、0.5 で乗算するか 2 で除算するbitshiftではなく、ビット シフト右算術関数 (bitsra) で実行されます。これは、bitsra が以下の条件を満たすからです。

  • より効率的な組み込みコードを生成します。

  • 正の数値でも負の数値でも、同じように機能します。

  • 浮動小数点数型、固定小数点数型、および整数型で同じように機能します。

  • このコードをデータ型から独立したままに保ちます。

添字を使用して変数 a(:) = b に代入すること (subsasgn) と変数 a = b を上書きすることには、違いがあることに注意してください。変数に対して以下のように添字を使用して代入するとします。

x(:) = x + bitsra(y, i);

これにより、左辺の引数 x のデータ型が常に保持されます。固定小数点型では、これこそまさに推奨されるプログラミング スタイルです。たとえば、固定小数点型では和の語長が大きくなることがよくあります。これは、fimathオブジェクトの SumMode プロパティで決まりますが、この方法によって右辺 x + bitsra(y,i) のデータ型は x とは別のデータ型にすることができます。

今度は、左辺を以下のように上書きしてみましょう。

x = x + bitsra(y, i);

すると、左辺 x のデータ型は右辺の合計の型と同じになります。このプログラミング スタイルだと、固定小数点型コードでの x のデータ型が変更されてしまうので、お勧めできません。

逆 CORDIC 成長定数の定義

関数 inverse_cordic_growth_constant は、反復を niter 回実行した後、CORDIC 成長係数の逆数を返します。これが必要とされるのは、1 回の CORDIC 回転によって値が約 1.6468 倍になり、反復回数に応じて値が成長していくからです。したがって、cordicgivens の最後のステップで逆数 Kn = 1/1.6468 = 0.60725 を乗算することによって、ゲインが正規化されます。

function Kn = inverse_cordic_growth_constant(niter)
  Kn = 1/prod(sqrt(1+2.^(-2*(0:double(niter)-1))));
end

反復回数の関数としての CORDIC 成長率の検討

CORDIC 成長率の関数は、以下のように定義されます。

growth = prod(sqrt(1+2.^(-2*(0:double(niter)-1))))

そして、逆は以下のように定義されます。

inverse_growth = 1 ./ growth

成長率は反復回数 niter の関数であり、約 1.6468 にすぐに収束します。逆は約 0.60725 に収束します。27 回反復した後に 1 回の反復から次の停止までの差が変化することが次の表からわかります。これは、27 回目の反復で計算が double 浮動小数点型の精度の限界に達したためです。

niter       growth            diff(growth)         1./growth        diff(1./growth)
  0    1.000000000000000                   0   1.000000000000000                   0
  1    1.414213562373095   0.414213562373095   0.707106781186547  -0.292893218813453
  2    1.581138830084190   0.166925267711095   0.632455532033676  -0.074651249152872
  3    1.629800601300662   0.048661771216473   0.613571991077896  -0.018883540955780
  4    1.642484065752237   0.012683464451575   0.608833912517752  -0.004738078560144
  5    1.645688915757255   0.003204850005018   0.607648256256168  -0.001185656261584
  6    1.646492278712479   0.000803362955224   0.607351770141296  -0.000296486114872
  7    1.646693254273644   0.000200975561165   0.607277644093526  -0.000074126047770
  8    1.646743506596901   0.000050252323257   0.607259112298893  -0.000018531794633
  9    1.646756070204878   0.000012563607978   0.607254479332562  -0.000004632966330
 10    1.646759211139822   0.000003140934944   0.607253321089875  -0.000001158242687
 11    1.646759996375617   0.000000785235795   0.607253031529134  -0.000000289560741
 12    1.646760192684695   0.000000196309077   0.607252959138945  -0.000000072390190
 13    1.646760241761972   0.000000049077277   0.607252941041397  -0.000000018097548
 14    1.646760254031292   0.000000012269320   0.607252936517010  -0.000000004524387
 15    1.646760257098622   0.000000003067330   0.607252935385914  -0.000000001131097
 16    1.646760257865455   0.000000000766833   0.607252935103139  -0.000000000282774
 17    1.646760258057163   0.000000000191708   0.607252935032446  -0.000000000070694
 18    1.646760258105090   0.000000000047927   0.607252935014772  -0.000000000017673
 19    1.646760258117072   0.000000000011982   0.607252935010354  -0.000000000004418
 20    1.646760258120067   0.000000000002995   0.607252935009249  -0.000000000001105
 21    1.646760258120816   0.000000000000749   0.607252935008973  -0.000000000000276
 22    1.646760258121003   0.000000000000187   0.607252935008904  -0.000000000000069
 23    1.646760258121050   0.000000000000047   0.607252935008887  -0.000000000000017
 24    1.646760258121062   0.000000000000012   0.607252935008883  -0.000000000000004
 25    1.646760258121065   0.000000000000003   0.607252935008882  -0.000000000000001
 26    1.646760258121065   0.000000000000001   0.607252935008881  -0.000000000000000
 27    1.646760258121065                   0   0.607252935008881                   0
 28    1.646760258121065                   0   0.607252935008881                   0
 29    1.646760258121065                   0   0.607252935008881                   0
 30    1.646760258121065                   0   0.607252935008881                   0
 31    1.646760258121065                   0   0.607252935008881                   0
 32    1.646760258121065                   0   0.607252935008881                   0

CORDIC と標準ギブンス回転との比較

関数 cordicgivens は、数値的には Golub と Van Loan の "Matrix Computations" [3] の以下の標準ギブンス回転アルゴリズムと等価です。関数 cordicqr で、cordicgivens の呼び出しを givensrotation の呼び出しで置換すると、標準ギブンス QR アルゴリズムが得られます。

function [x,y,u,v] = givensrotation(x,y,u,v)
  a = x(1); b = y(1);
  if b==0
    % No rotation necessary.  c = 1; s = 0;
    return;
  else
    if abs(b) > abs(a)
      t = -a/b; s = 1/sqrt(1+t^2); c = s*t;
    else
      t = -b/a; c = 1/sqrt(1+t^2); s = c*t;
    end
  end
  x0 = x;             u0 = u;
  % x and y form R,   u and v form Q
  x(:) = c*x0 - s*y;  u(:) = c*u0 - s*v;
  y(:) = s*x0 + c*y;  v(:) = s*u0 + c*v;
end

関数 givensrotation は、除算と平方根演算を使用します。これらは固定小数点型では高コストですが、浮動小数点型アルゴリズムには適しています。

CORDIC 回転の例

アルゴリズムの各ステップで CORDIC 回転に従う 3 行 3 列の行列の例を以下に示します。このアルゴリズムでは、直交回転を使用して R の下対角要素をゼロに設定します。その際、対角要素をピボットとして使用します。この同じ回転は、単位行列にも適用され、Q*R = A であるような直交行列 Q が形成されます。

A をランダムな 3 行 3 列の行列とします。R = A および Q = eye(3) として開始値を設定します。

  R = A = [-0.8201    0.3573   -0.0100
           -0.7766   -0.0096   -0.7048
           -0.7274   -0.6206   -0.8901]
      Q = [ 1         0         0
            0         1         0
            0         0         1]

最初の回転は、R の 1 行目と 2 行目、および Q の 1 列目と 2 列目について実行されます。要素 R(1,1) がピボットであり、R(2,1) は回転して 0 になります。

    R before the first rotation            R after the first rotation
x [-0.8201    0.3573   -0.0100]   ->    x [1.1294   -0.2528    0.4918]
y [-0.7766   -0.0096   -0.7048]   ->    y [     0    0.2527    0.5049]
   -0.7274   -0.6206   -0.8901            -0.7274   -0.6206   -0.8901
    Q before the first rotation            Q after the first rotation
    u         v                            u         v
   [1]       [0]        0                [-0.7261] [ 0.6876]        0
   [0]       [1]        0         ->     [-0.6876] [-0.7261]        0
   [0]       [0]        1                [      0] [      0]        1

以下のプロットでは、CORDIC 反復を 1 回実行するたびに x が成長することがわかります。この成長率は、最後のステップで Kn = 0.60725 を乗算することにより抽出されます。y(1) が反復されてゼロになることがわかります。最初は、点 [x(1), y(1)] は第 3 象限にあり、CORDIC 反復の開始前に第 1 象限に鏡映されます。

2 回目の回転は、R の 1 行目と 3 行目、および Q の 1 列目と 3 列目について実行されます。要素 R(1,1) がピボットであり、R(3,1) は回転して 0 になります。

    R before the second rotation           R after the second rotation
x  [1.1294   -0.2528    0.4918]   ->    x [1.3434    0.1235    0.8954]
         0    0.2527    0.5049                  0    0.2527    0.5049
y [-0.7274]  -0.6206   -0.8901    ->    y [     0   -0.6586   -0.4820]
    Q before the second rotation           Q after the second rotation
    u                       v              u                   v
  [-0.7261]   0.6876       [0]           [-0.6105]   0.6876  [-0.3932]
  [-0.6876]  -0.7261       [0]    ->     [-0.5781]  -0.7261  [-0.3723]
  [      0]        0       [1]           [-0.5415]        0  [ 0.8407]

3 回目の回転は、R の 2 行目と 3 行目、および Q の 2 列目と 3 列目について実行されます。要素 R(2,2) がピボットであり、R(3,2) は回転して 0 になります。

    R before the third rotation            R after the third rotation
    1.3434    0.1235    0.8954             1.3434    0.1235    0.8954
 x       0  [ 0.2527    0.5049]   ->    x       0   [0.7054    0.6308]
 y       0  [-0.6586   -0.4820]   ->    y       0   [     0    0.2987]
    Q before the third rotation            Q after the third rotation
              u         v                            u         v
   -0.6105  [ 0.6876] [-0.3932]           -0.6105  [ 0.6134] [ 0.5011]
   -0.5781  [-0.7261] [-0.3723]   ->      -0.5781  [ 0.0875] [-0.8113]
   -0.5415  [      0] [ 0.8407]           -0.5415  [-0.7849] [ 0.3011]

これで QR 分解が完了しました。R は上三角行列、Q は直交行列です。

R =
   1.3434    0.1235    0.8954
        0    0.7054    0.6308
        0         0    0.2987
Q =
  -0.6105    0.6134    0.5011
  -0.5781    0.0875   -0.8113
  -0.5415   -0.7849    0.3011

乗算によって直交行列になるため、Q が丸め誤差の範囲内であり、しかも単位行列に近いことを確認して、検証できます。

Q*Q' =  1.0000    0.0000    0.0000
        0.0000    1.0000         0
        0.0000         0    1.0000
Q'*Q =  1.0000    0.0000   -0.0000
        0.0000    1.0000   -0.0000
       -0.0000   -0.0000    1.0000

誤差は単位行列を減算することにより確認することができます。

Q*Q' - eye(size(Q)) =           0   2.7756e-16   3.0531e-16
                       2.7756e-16   4.4409e-16            0
                       3.0531e-16            0   6.6613e-16

Q*RA に近いことは、減算して誤差を確認することにより検証することができます。

Q*R - A =  -3.7802e-11  -7.2325e-13  -2.7756e-17
           -3.0512e-10   1.1708e-12  -4.4409e-16
            3.6836e-10  -4.3487e-13  -7.7716e-16

固定語長向けの Q の最適な出力型の判断

Q は直交行列なので、値はすべて -1 と +1 の間になります。浮動小数点数では、Q の型について判断することは何もありません。A と同じ浮動小数点型になるからです。ただし、固定小数点型では、QA と同じ固定小数点型をもたせるよりもさらに優れた方法があります。たとえば、A の語長が 16、小数部の長さが 8 である場合、Q についても語長 16、小数部の長さ 8 とすると、Q の精度を本来より強制的に下げることになり、固定小数点範囲の上半分が無駄になってしまいます。

Q にとって最善の型は、可能な出力の全範囲をもたせることに加えて、中間計算に CORDIC 成長係数 1.6468 を組み入れることです。したがって、Q の語長が入力 A の語長と同じであると仮定すると、Q にとって最善の小数部の長さは語長から 2 ビット引いた値 (1 ビットは 1.6468 の分、もう 1 ビットは符号の分) となります。

したがって、関数 cordicqr における Q の開始値の設定を以下のように改良することができます。

if isfi(A) && (isfixed(A) || isscaleddouble(A))
      Q = fi(one*eye(m), get(A,'NumericType'), ...
             'FractionLength',get(A,'WordLength')-2);
else
  Q = coder.nullcopy(repmat(A(:,1),1,m));
  Q(:) = eye(m);
end

わずかに不便な点があるとすれば、コードの該当のセクションがデータ型に依存してしまうことです。ただし、Q にとって最適な型を選択することにより大きなメリットが得られ、アルゴリズムの主要部は依然としてデータ型から独立しています。この種の入力解析を関数の開始時に実行して、アルゴリズムの主要部をデータ型から独立したままにすることができます。

固定小数点型 R におけるオーバーフローの回避

このセクションでは、オーバーフローを回避できるように R の固定小数点出力型を決定する方法について説明します。出力型を選択するには、R の値の大きさがどれほど大きくなるのかを知る必要があります。

実数行列 A とその QR 分解が、ピボットなしのギブンス回転により計算されるとします。R の要素の大きさの上限は、A の行数の平方根に A の最大要素の大きさをかけた値になります。さらに、1 回の中間計算中にこの成長率がこれより大きくなることはありません。つまり、[m,n] = size(A) かつ [Q,R] = givensqr(A) だとします。以下のようになります。

max(abs(R(:))) <= sqrt(m) * max(abs(A(:))).

これが成り立つのは、R の各要素が A の対応する列からの直交回転により形成されるからです。したがって、任意の要素 R(i,j) が得ることができる最大値は、対応する列 A(:,j) の要素すべてが回転して 1 つの値になるときの値です。つまり、取り得る最大値は、A(:,j) の 2 ノルムによって上限が設定されます。A(:,j) の 2 ノルムは、m 個の要素の二乗和の平方根に等しく、かつ各要素は A の最大要素以下なので、以下のようになります。

norm(A(:,j)) <= sqrt(m) * max(abs(A(:))).

つまり、すべての j について、以下のようになります。

norm(A(:,j))  = sqrt(A(1,j)^2 + A(2,j)^2 + ... + A(m,j)^2)
             <= sqrt( m * max(abs(A(:)))^2)
              = sqrt(m) * max(abs(A(:))).

したがって、すべての i,j について、以下のようになります。

abs(R(i,j)) <= norm(A(:,j)) <= sqrt(m) * max(abs(A(:))).

したがって、これは R の最大要素についても成り立ちます。

max(abs(R(:))) <= sqrt(m) * max(abs(A(:))).

これは、データ型が到達できる最大値に A の要素が非常に近くなることがよくある固定小数点型で役立ちます。その結果、A の値を知らなくとも、差がほとんどない上限を設定することができます。これが重要なのは、A のデータ型の上限しか知らない状況で、最小のビット数で R の出力型を設定する必要があるからです。この値は、fi のメソッドupperboundを使用して取得できます。

そのため、すべての i,j について、以下のようになります。

abs(R(i,j)) <= sqrt(m) * upperbound(A)

sqrt(m)*upperbound(A)A の要素の上限であることに注意してください。

abs(A(i,j)) <= upperbound(A) <= sqrt(m)*upperbound(A)

したがって、固定小数点データ型を選択すると、sqrt(m)*upperbound(A)AR の両方で機能する上限になります。

最大値に達することは簡単で、よくあることです。最大値に達するのは、直交列をもつ以下の行列のように、要素がすべて回転されて 1 つの要素にまとまるときです。

A = [7    -7     7     7
     7     7    -7     7
     7    -7    -7    -7
     7     7     7    -7];

最大値は 7 であり、行数は m = 4 です。したがって、R の最大値は max(abs(A(:)))*sqrt(m) = 7*sqrt(4) = 14 によって上限が設定されます。この例の A は直交行列なので、各列は回転されて対角上の最大値になります。

niter = 52;
[Q,R] = cordicqr(A,niter)
Q = 4×4

    0.5000   -0.5000    0.5000    0.5000
    0.5000    0.5000   -0.5000    0.5000
    0.5000   -0.5000   -0.5000   -0.5000
    0.5000    0.5000    0.5000   -0.5000

R = 4×4

   14.0000    0.0000   -0.0000   -0.0000
         0   14.0000   -0.0000    0.0000
         0         0   14.0000    0.0000
         0         0         0   14.0000

最大成長率に達する別の簡単な例は、すべての要素が等しい行列です。たとえば、要素がすべて 1 である行列です。1 の行列は、回転されて最初の行は 1*sqrt(m) になり、それ以外はゼロになります。たとえば、以下の 9 行 5 列の行列では、R の最初の行はすべて 1*sqrt(9)=3 になります。

m = 9; n = 5;
A = ones(m,n)
A = 9×5

     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1
     1     1     1     1     1

niter = 52;
[Q,R] = cordicqr(A,niter)
Q = 9×9

    0.3333    0.5567   -0.6784    0.3035   -0.1237    0.0503    0.0158    0.0056   -0.0921
    0.3333    0.0296    0.2498   -0.1702   -0.6336    0.1229   -0.3012   -0.5069   -0.1799
    0.3333    0.2401    0.0562   -0.3918    0.4927    0.2048   -0.5395    0.0359    0.3122
    0.3333    0.0003    0.0952   -0.1857    0.2148    0.4923    0.7080   -0.2351   -0.0175
    0.3333    0.1138    0.0664   -0.2263    0.1293   -0.8348    0.2510   -0.2001    0.0610
    0.3333   -0.3973   -0.0143    0.3271    0.4132   -0.0354   -0.2165   -0.0939   -0.6294
    0.3333    0.1808    0.3538   -0.1012   -0.2195         0    0.0824    0.7646   -0.2849
    0.3333   -0.6500   -0.4688   -0.2380   -0.2400         0         0    0.2300    0.2820
    0.3333   -0.0740    0.3400    0.6825   -0.0331         0         0         0    0.5485

R = 9×5

    3.0000    3.0000    3.0000    3.0000    3.0000
         0    0.0000    0.0000    0.0000    0.0000
         0         0    0.0000    0.0000    0.0000
         0         0         0    0.0000    0.0000
         0         0         0         0    0.0000
         0         0         0         0         0
         0         0         0         0         0
         0         0         0         0         0
         0         0         0         0         0

関数 cordicqr の場合と同様に、ギブンス QR アルゴリズムは A の位置に R を代入して上書きして記述されることがよくあります。そのため、アルゴリズムの開始時に AR のデータ型にキャストできると便利です。

さらに、CORDIC でギブンス回転を計算する場合、約 1.6468 にすぐに収束する成長係数があります。この成長係数は、ギブンス回転が 1 回実行されるたびに正規化されますが、これを中間計算に組み入れる必要があります。したがって、ギブンス回転と CORDIC 成長率を含め、必要とされる追加ビット数は、log2(1.6468*sqrt(m)) になります。先頭空間の追加ビットは、語長を増やすか小数部の長さを短くすることで追加することができます。

語長を増やすことのメリットは、所定の語長に対して取り得る最大の精度が考慮されることです。デメリットは、最適なワード サイズがプロセッサ上でのネイティブ型に対応していない可能性があることです (たとえば、16 ビットから 18 ビットに増やすなど)。または、一段上の大きいネイティブ ワード サイズ) に増やさなければならず、これが非常に大きい値になる可能性があることです (たとえば、18 ビットで十分なのに、16 ビットから 32 ビットに増やすなど)。

小数部の長さを短くすることのメリットは、計算を A のネイティブ ワード サイズで実行できることです。デメリットは、精度が低下することです。

別の選択肢は、右シフトを実行して入力をあらかじめスケール変更しておくことです。これは、小数部の長さを短くすることと同じなので、問題のスケーリングを変更する点が不便です。ただし、分数演算または整数演算でのみ処理する場合には、この選択肢の方が便利かもしれません。

R における固定小数点成長率の例

固定小数点入力行列 A がある場合、前のセクションで定義された成長率を使用して固定小数点出力 R を定義することができます。

ランダム行列 X で始めます。

X = [0.0513   -0.2097    0.9492    0.2614
     0.8261    0.6252    0.3071   -0.9415
     1.5270    0.1832    0.1352   -0.1623
     0.4669   -1.0298    0.5152   -0.1461];

X から固定小数点 A を作成します。

A = fi(X,1)
A = 
    0.0513   -0.2097    0.9492    0.2614
    0.8261    0.6252    0.3071   -0.9415
    1.5270    0.1832    0.1352   -0.1623
    0.4669   -1.0298    0.5152   -0.1461

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 14
m = size(A,1)
m = 4

成長係数は、1.6468 に A の行数の平方根をかけた値です。ビット成長率は、2 を底とする成長率の対数より大きい次の整数です。

bit_growth = ceil(log2(cordic_growth_constant*sqrt(m)))
bit_growth = 2

R の開始値を A と同じ値に設定します。すると、語長がビット成長率によって増大します。

R = fi(A,1,get(A,'WordLength')+bit_growth,get(A,'FractionLength'))
R = 
    0.0513   -0.2097    0.9492    0.2614
    0.8261    0.6252    0.3071   -0.9415
    1.5270    0.1832    0.1352   -0.1623
    0.4669   -1.0298    0.5152   -0.1461

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

R を入力として使用して、上書きします。

niter = get(R,'WordLength') - 1
niter = 17
[Q,R] = cordicqr(R, niter)
Q = 
    0.0284   -0.1753    0.9110    0.3723
    0.4594    0.4470    0.3507   -0.6828
    0.8490    0.0320   -0.2169    0.4808
    0.2596   -0.8766   -0.0112   -0.4050

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 18
        FractionLength: 16
R = 
    1.7989    0.1694    0.4166   -0.6008
         0    1.2251   -0.4764   -0.3438
         0         0    0.9375   -0.0555
         0         0         0    0.7214

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

Q*Q' が単位行列に近いことを確認します。

double(Q)*double(Q')
ans = 4×4

    1.0000   -0.0001    0.0000    0.0000
   -0.0001    1.0001    0.0000   -0.0000
    0.0000    0.0000    1.0000   -0.0000
    0.0000   -0.0000   -0.0000    1.0000

Q*R - AA の精度に対して相対的に小さいことを確認します。

err = double(Q)*double(R) - double(A)
err = 4×4
10-3 ×

   -0.1048   -0.2355    0.1829   -0.2146
    0.3472    0.2949    0.0260   -0.2570
    0.2776   -0.1740   -0.1007    0.0966
    0.0138   -0.1558    0.0417   -0.0362

R における精度の向上

前のセクションで、R でのオーバーフローを回避しながら A の精度を維持する方法を説明しました。R の小数部の長さを A と同じにしておくと、R の精度が A の精度を超えることはできないのですが、R の精度がより高くなければならないような精度要件があることがあります。

これの極端な例は、整数固定小数点型 (つまり、小数部の長さがゼロ) で行列を定義することです。行列 X の要素を符号付き 8 ビット整数の全範囲、つまり -128 と +127 の間の値であるとします。

 X = [-128  -128  -128   127
      -128   127   127  -128
       127   127   127   127
       127   127  -128  -128];

固定小数点行列 A を 8 ビット整数と等価になるように定義します。

A = fi(X,1,8,0)
A = 
  -128  -128  -128   127
  -128   127   127  -128
   127   127   127   127
   127   127  -128  -128

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 8
        FractionLength: 0
m = size(A,1)
m = 4

必要な成長率は、1.6468 に A の行数の平方根をかけた値です。

bit_growth = ceil(log2(cordic_growth_constant*sqrt(m)))
bit_growth = 2

R の開始値を A と同じ値に設定し、前のセクションと同様にビット成長率を考慮に入れます。

R = fi(A,1,get(A,'WordLength')+bit_growth,get(A,'FractionLength'))
R = 
  -128  -128  -128   127
  -128   127   127  -128
   127   127   127   127
   127   127  -128  -128

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 10
        FractionLength: 0

QR 分解を計算し、R を上書きします。

niter = get(R,'WordLength') - 1;
[Q,R] = cordicqr(R, niter)
Q = 
   -0.5039   -0.2930   -0.4062   -0.6914
   -0.5039    0.8750    0.0039    0.0078
    0.5000    0.2930    0.3984   -0.7148
    0.4922    0.2930   -0.8203    0.0039

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 10
        FractionLength: 8
R = 
   257   126    -1    -1
     0   225   151  -148
     0     0   211   104
     0     0     0  -180

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 10
        FractionLength: 0

R が整数値で返されたことに注意してください。これは、R の小数部の長さを A の小数部の長さと同じ 0 のままにしたためです。

A の最下位ビット (LSB) のスケーリングは 1 であり、誤差は LSB に比例することがわかります。

err = double(Q)*double(R)-double(A)
err = 4×4

   -1.5039   -1.4102   -1.4531   -0.9336
   -1.5039    6.3828    6.4531   -1.9961
    1.5000    1.9180    0.8086   -0.7500
   -0.5078    0.9336   -1.3398   -1.8672

小数部の長さを伸ばすことによって、QR 分解の精度を上げることができます。前の例では、整数部に 10 ビット (開始に 8 ビット、成長率に 2 ビット) が必要でした。したがって、小数部の長さを伸ばす場合にも、整数部に 10 ビットを維持する必要があります。たとえば、語長を 32 に増やし、小数部の長さを 22 に設定することができます。これにより、整数部に 10 ビットが確保されます。

R = fi(A,1,32,22)
R = 
  -128  -128  -128   127
  -128   127   127  -128
   127   127   127   127
   127   127  -128  -128

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 32
        FractionLength: 22
niter = get(R,'WordLength') - 1;
[Q,R] = cordicqr(R, niter)
Q = 
   -0.5020   -0.2913   -0.4088   -0.7043
   -0.5020    0.8649    0.0000    0.0000
    0.4980    0.2890    0.4056   -0.7099
    0.4980    0.2890   -0.8176    0.0000

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 32
        FractionLength: 30
R = 
  255.0020  127.0029    0.0039    0.0039
         0  220.5476  146.8413 -147.9930
         0         0  208.4793  104.2429
         0         0         0 -179.6037

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 32
        FractionLength: 22

これで、R の小数部を確認でき、Q*R-A が小さいことがわかります。

err = double(Q)*double(R)-double(A)
err = 4×4
10-5 ×

   -0.1234   -0.0014   -0.0845    0.0267
   -0.1234    0.2574    0.1260   -0.1094
    0.0720    0.0289   -0.0400   -0.0684
    0.0957    0.0818   -0.1034    0.0095

小数部の長さに指定するビット数は、使用するアルゴリズムの精度要件で決まります。

既定の反復回数の選択

反復回数は、目的の精度で決まりますが、A の語長によって制限されます。1 回反復するたびに、値は 1 ビットだけ右にシフトします。最終ビットがシフトされ、値がゼロになると、それ以上回転させるための追加の値はなくなります。したがって、最大限の精度を得るには、niter を語長より 1 ビットだけ小さい値にします。

浮動小数点型の場合、反復回数は仮数のサイズにより限界が設定されます。double では、同じ指数をもつ何かに続けて追加するためにできる最多の反復回数は 52 回です。single では 23 回です。浮動小数点の精度の詳細については、epsを参照してください。

このように、反復回数を入力しなくても済むようにし、niter にこの既定の設定を使用するよう関数 cordicqr を変更して最大限の精度を得ることによって、コードをさらに使いやすくすることができます。

function [Q,R] = cordicqr(A,varargin)
  if nargin>=2 && ~isempty(varargin{1})
     niter = varargin{1};
  elseif isa(A,'double') || isfi(A) && isdouble(A)
    niter = 52;
  elseif isa(A,'single') || isfi(A) && issingle(A)
    niter = single(23);
  elseif isfi(A)
    niter = int32(get(A,'WordLength') - 1);
  else
    assert(0,'First input must be double, single, or fi.');
  end

この方法の不便な点は、コードの一部のセクションがデータ型に依存するようになってしまうことです。ただし、必ずしも niter を指定する必要がなく、しかもアルゴリズムの主要部は依然としてデータ型から独立しているので、関数の利便性が高まるというメリットがあります。Q の最適な出力型を選択する場合と同様に、この種の入力解析を関数の開始時に実行して、アルゴリズムの主要部をデータ型から独立したままにすることができます。

前の節に示した例を再度示します。最適な niter を指定する必要がないことに注意してください。

A = [7    -7     7     7
     7     7    -7     7
     7    -7    -7    -7
     7     7     7    -7];
[Q,R] = cordicqr(A)
Q = 4×4

    0.5000   -0.5000    0.5000    0.5000
    0.5000    0.5000   -0.5000    0.5000
    0.5000   -0.5000   -0.5000   -0.5000
    0.5000    0.5000    0.5000   -0.5000

R = 4×4

   14.0000    0.0000   -0.0000   -0.0000
         0   14.0000   -0.0000    0.0000
         0         0   14.0000    0.0000
         0         0         0   14.0000

例: 一意でない QR 分解

MATLAB の関数 QR の結果と関数 cordicqr の結果を比較すると、QR 分解が一意ではないことに気づきます。重要なのは、Q が直交行列、R が上三角行列、Q*R - A が小さいということだけです。

その差を示す簡単な例を示します。

m = 3;
A = ones(m)
A = 3×3

     1     1     1
     1     1     1
     1     1     1

MATLAB の組み込み QR 関数では、別のアルゴリズムが使用されます。

[Q0,R0] = qr(A)
Q0 = 3×3

   -0.5774    0.8165   -0.0000
   -0.5774   -0.4082   -0.7071
   -0.5774   -0.4082    0.7071

R0 = 3×3

   -1.7321   -1.7321   -1.7321
         0   -0.0000   -0.0000
         0         0    0.0000

関数 cordicqr と比較します。

[Q,R] = cordicqr(A)
Q = 3×3

    0.5774    0.7495    0.3240
    0.5774   -0.6553    0.4871
    0.5774   -0.0942   -0.8110

R = 3×3

    1.7321    1.7321    1.7321
         0    0.0000    0.0000
         0         0   -0.0000

関数 cordicqr による Q の要素が、組み込み QR による Q0 と異なることに注意してください。ただし、どちらの結果も Q は直交行列である、という要件を満たしています。

Q0*Q0'
ans = 3×3

    1.0000    0.0000   -0.0000
    0.0000    1.0000   -0.0000
   -0.0000   -0.0000    1.0000

Q*Q'
ans = 3×3

    1.0000    0.0000    0.0000
    0.0000    1.0000   -0.0000
    0.0000   -0.0000    1.0000

しかも、どちらも Q*R - A が小さい、という要件を満たしています。

Q0*R0 - A
ans = 3×3
10-15 ×

   -0.1110   -0.1110   -0.1110
   -0.1110   -0.1110   -0.1110
   -0.1110   -0.1110   -0.1110

Q*R - A
ans = 3×3
10-15 ×

   -0.2220    0.2220    0.2220
    0.4441         0         0
    0.2220    0.2220    0.2220

Q を形成しない方程式系の求解

行列 AB が与えられた場合、QR 分解を使用して以下の式の X を求めることができます。

A*X = B.

A において列数より行数の方が多い場合、X は最小二乗解になります。XB の列数が 1 を超える場合、一度に複数の解を求めることができます。A = Q*RA の QR 分解である場合、以下をバックソルブすることによって解を求めることができます。

R*X = C

ただし、C = Q'*B です。Q を形成し乗算して C = Q'*B を得るよりも、C を直接計算する方が効率的です。C を直接計算するには、単位行列の列ではなく B の行に回転を適用します。C = B の開始値を設定し、Q の列ではなく C の行を演算するというように少し変更するだけで、新しいアルゴリズムを作成することができます。

  function [R,C] = cordicrc(A,B,niter)
    Kn = inverse_cordic_growth_constant(niter);
    [m,n] = size(A);
    R = A;
    C = B;
    for j=1:n
      for i=j+1:m
        [R(j,j:end),R(i,j:end),C(j,:),C(i,:)] = ...
            cordicgivens(R(j,j:end),R(i,j:end),C(j,:),C(i,:),niter,Kn);
      end
    end
  end

アルゴリズムをこの例を使用して検証することができます。A をランダムな 3 行 3 列の行列、B をランダムな 3 行 2 列の行列であるとします。

A = [-0.8201    0.3573   -0.0100
     -0.7766   -0.0096   -0.7048
     -0.7274   -0.6206   -0.8901];

B = [-0.9286    0.3575
      0.6983    0.5155
      0.8680    0.4863];

A の QR 分解を計算します。

[Q,R] = cordicqr(A)
Q = 3×3

   -0.6105    0.6133    0.5012
   -0.5781    0.0876   -0.8113
   -0.5415   -0.7850    0.3011

R = 3×3

    1.3434    0.1235    0.8955
         0    0.7054    0.6309
         0         0    0.2988

C = Q'*B を直接計算します。

[R,C] = cordicrc(A,B)
R = 3×3

    1.3434    0.1235    0.8955
         0    0.7054    0.6309
         0         0    0.2988

C = 3×2

   -0.3068   -0.7795
   -1.1897   -0.1173
   -0.7706   -0.0926

減算します。すると、誤差が丸め誤差にほぼ等しいことがわかります。

Q'*B - C
ans = 3×2
10-15 ×

   -0.0555    0.3331
         0         0
    0.1110    0.2914

それでは、この同じ例を固定小数点型で試してみましょう。AB を固定小数点型として宣言します。

A = fi(A,1)
A = 
   -0.8201    0.3573   -0.0100
   -0.7766   -0.0096   -0.7048
   -0.7274   -0.6206   -0.8901

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 15
B = fi(B,1)
B = 
   -0.9286    0.3575
    0.6983    0.5155
    0.8680    0.4863

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 16
        FractionLength: 15

必要な成長率は、1.6468 に A の行数の平方根をかけた値です。

bit_growth = ceil(log2(cordic_growth_constant*sqrt(m)))
bit_growth = 2

R の開始値を A と同じ値に設定し、ビット成長率を考慮に入れます。

R = fi(A,1,get(A,'WordLength')+bit_growth,get(A,'FractionLength'))
R = 
   -0.8201    0.3573   -0.0100
   -0.7766   -0.0096   -0.7048
   -0.7274   -0.6206   -0.8901

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

C の成長率は R と同じなので、同じように C の開始値を設定しビット成長率を考慮に入れます。

C = fi(B,1,get(B,'WordLength')+bit_growth,get(B,'FractionLength'))
C = 
   -0.9286    0.3575
    0.6983    0.5155
    0.8680    0.4863

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

C = Q'*B を直接計算し、RC を上書きします。

[R,C] = cordicrc(R,C)
R = 
    1.3435    0.1233    0.8954
         0    0.7055    0.6308
         0         0    0.2988

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 18
        FractionLength: 15
C = 
   -0.3068   -0.7796
   -1.1898   -0.1175
   -0.7706   -0.0926

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

このアルゴリズムの興味深い使い方は、B に開始値を設定して単位行列にすると、出力引数 CQ' になることです。この機能を使用して、Q のデータ型をさらに自由に操れるようにしても良いでしょう。次に例を示します。

A = [-0.8201    0.3573   -0.0100
     -0.7766   -0.0096   -0.7048
     -0.7274   -0.6206   -0.8901];
B = eye(size(A,1))
B = 3×3

     1     0     0
     0     1     0
     0     0     1

[R,C] = cordicrc(A,B)
R = 3×3

    1.3434    0.1235    0.8955
         0    0.7054    0.6309
         0         0    0.2988

C = 3×3

   -0.6105   -0.5781   -0.5415
    0.6133    0.0876   -0.7850
    0.5012   -0.8113    0.3011

すると、C は直交行列になります。

C'*C
ans = 3×3

    1.0000    0.0000    0.0000
    0.0000    1.0000   -0.0000
    0.0000   -0.0000    1.0000

そして、R = C*A は以下のようになります。

R - C*A
ans = 3×3
10-15 ×

    0.6661   -0.0139   -0.1110
    0.5551   -0.2220    0.6661
   -0.2220   -0.1110    0.2776

この例で使用された関数

この例で使用された MATLAB 関数は以下のとおりです。

関数 CORDICQR は CORDIC を使用して QR 分解を計算します。

  • [Q,R] = cordicqr(A) は、A の型に基づいて CORDIC の反復回数を選択します。

  • [Q,R] = cordicqr(A,niter) は、CORDIC 反復回数 niter を使用します。

関数 cordicrc は、A の QR 分解から R を計算し、Q を計算することなく C = Q'*B を返します。

  • [R,C] = cordicrc(A,B) は、A の型に基づいて CORDIC の反復回数を選択します。

  • [R,C] = cordicrc(A,B,niter) は、CORDIC 反復回数 niter を使用します。

関数 cordic_growth_constant は、CORDIC 成長定数を返します。

  • cordic_growth = cordic_growth_constant(niter) は、CORDIC 成長定数を CORDIC 反復回数 niter の関数として返します。

関数 GIVENSQR は 標準ギブンス回転を使用して QR 分解を計算します。

  • A が M 行 N 列の行列である場合、[Q,R] = givensqr(A) は、M 行 N 列の上三角行列 R および A = Q*R であるような M 行 M 列の直交行列 Q を生成します。

CORDICQR_MAKEPLOTS は MATLAB コマンド ラインから以下のコマンドを実行することによって、この例のプロットを作成します。

load A_3_by_3_for_cordicqr_demo.mat
niter = 32;
[Q,R] = cordicqr_makeplots(A,niter)

参考文献

  1. Ray Andraka, "A survey of CORDIC algorithms for FPGA based computers," 1998, ACM 0-89791-978-5/98/01.

  2. Anthony J Cox and Nicholas J Higham, "Stability of Householder QR factorization for weighted least squares problems," in Numerical Analysis, 1997, Proceedings of the 17th Dundee Conference, Griffiths DF, Higham DJ, Watson GA (eds). Addison-Wesley, Longman: Harlow, Essex, U.K., 1998; 57-73.

  3. Gene H. Golub and Charles F. Van Loan, Matrix Computations, 3rd ed, Johns Hopkins University Press, 1996, section 5.2.3 Givens QR Methods.

  4. Daniel V. Rabinkin, William Song, M. Michael Vai, and Huy T. Nguyen, "Adaptive array beamforming with fixed-point arithmetic matrix inversion using Givens rotations," Proceedings of Society of Photo-Optical Instrumentation Engineers (SPIE) -- Volume 4474 Advanced Signal Processing Algorithms, Architectures, and Implementations XI, Franklin T. Luk, Editor, November 2001, pp. 294--305.

  5. Jack E. Volder, "The CORDIC Trigonometric Computing Technique," Institute of Radio Engineers (IRE) Transactions on Electronic Computers, September, 1959, pp. 330-334.

  6. Musheng Wei and Qiaohua Liu, "On growth factors of the modified Gram-Schmidt algorithm," Numerical Linear Algebra with Applications, Vol. 15, issue 7, September 2008, pp. 621-636.

%#ok<*MNEFF,*NASGU,*NOPTS,*ASGLU>

参考

|

関連するトピック