QR分解(Givens回転による、複素数の場合)

[←]　[↑]　[→]

QR分解のアルゴリズム

A=QR ただしQ自身でなくQ†bを求める

Aの要素は複素数とする。

【Givens回転、回転結果が実数】

/* AをRに入れる */

for j=1 to n

  for i=1 to m

    rij = aij

/* RとQ†bの計算 */

for j=1 to n

  for i=m-1 to j

    d = sqrt(|rij|2 + |ri+1j|2)

    c = rij/d , s = ri+1j/d

    rij = d , ri+1j = 0

    for k=j+1 to n

      tmp = c* rik + s* ri+1k

      ri+1k = -s rik + c ri+1k

      rik = tmp

    tmp = c* bi + s* bi+1

    bi+1 = -s bi + c bi+1

    bi = tmp

【Givens回転,回転行列の左上要素が実数,CORDIC】

/* AをRに入れる */

for j=1 to n

  for i=1 to m

    rij = aij

/* RとQ†bの計算 */

for j=1 to n

  for i=m-1 to j

    d = sqrt(|rij|2 + |ri+1j|2)

    θ= atan2(|rij|,|ri+1j|)

    φ1= atan2(Re(rij),Im(rij))

    φ2= atan2(Re(ri+1j),Im(ri+1j))

    φ=φ1 -φ2

    rij = d exp(iφ), ri+1j = 0

    for k=j+1 to n

      tmp = rik cosθ+ ri+1k exp(iφ)sinθ

      ri+1k = - rik sinθ + ri+1k exp(iφ)cosθ

      rik = tmp

    tmp = bi cosθ+ bi+1 exp(iφ)sinθ

    bi+1 = - bi sinθ+ bi+1 exp(iφ)cosθ

    bi = tmp

連立1次方程式Ax=bを解くには

/* Rx= Q†bをxについて解く(後退代入) */

xn = bn / rnn

for i=n-1 to 1

   tmp = 0

   for j=i+1 to n

      tmp = tmp + rij xj

   xi = (bi - tmp)/ rii

前項のGivensの回転行列を実数から複素数に拡張する。

◇

Givensの回転行列は「与えられた列ベクトルの1つの要素をゼロにする」ものであるから、複素数の列ベクトルに対応できるように偏角をつけて次のようにおく。

ただし、Qはユニタリ(正規直交)でなければならないから偏角は次の条件(1)を満たさなければならない。

さて、与えられた列ベクトル

にQをかけて要素z2をゼロにするには

したがってθは次の(2)で決まり

偏角は(1)と合わせて次の(3)を満たせばよいから、4つの偏角のうち2つを自由に選べる。

偏角の選び方として、計算量が減るようになるべく実数にした方がよいので、回転行列の1行1列目が実数になるように選ぶ(α=0 , β=φ1-φ2)、または回転結果が実数になるように選ぶ(α=-φ1 , β= -φ2)の2とおりが考えられる。

残りのγとδは次の2とおりが考えられる。回転行列の列ベクトルが同じ偏角になるように選ぶ(γ=α , δ=β)、または対角要素が同じ偏角になるように選ぶ(γ=-β , δ=-α)。

これら4つの場合を整理すると

 

列ベクトルが同じ偏角(γ=α , δ=β)

対角要素が同じ偏角(γ=-β , δ=-α)

回転行列の

(1,1)要素が実数

(α=0 , β=φ1-φ2)

回転結果が実数

(α=-φ1 , β= -φ2)

 

◇◇

4つの偏角の選び方のうち、式の対称性がよいのは次の場合である。

このとき、手順は次のようになる。

・Aを作業領域Rに入れる。

・j列(ｊ=1 to n)のステップにおいて

i行目(i=m-1 to j)をゼロにする行列G(ij)の

パラメータを求める

j列の回転結果は

j列の右側を回転する

bを回転する

【Givens法、ボトムアップ順、回転結果が実数】

/* AをRに入れる */

for j=1 to n

  for i=1 to m

    rij = aij

/* RとQ†bの計算 */

for j=1 to n

  for i=m-1 to j

    d = sqrt(|rij|2 + |ri+1j|2)

    c = rij/d , s = ri+1j/d

    rij = d , ri+1j = 0

    for k=j+1 to n

      tmp = c* rik + s* ri+1k

      ri+1k = -s rik + c ri+1k

      rik = tmp

    tmp = c* bi + s* bi+1

    bi+1 = -s bi + c bi+1

    bi = tmp

◇◇◇

演算量が少なくなるように選ぶならば、回転行列を求める計算よりも回転行列をかける方が回数が多いので、回転行列の要素が実数になるようにした方がよい。

そこで通常は、複素数のGivens回転行列として次を選ぶようである。

実数のGivens回転行列のときと同様に、固定小数点でよければCORDICが使える。

まず、与えられた列ベクトルzに対して

となる回転行列Qの回転角の計算はVectorモードを使い

求めたQによる回転はRotationモードを使えばよい。

ここで、複素数の偏角の回転もRotationモードで計算できることに注意。

[←]　[↑]　[→]