Sylvester行列のQR分解による最大公約多項式の計算過程についての覚書

（Greatest Common Divisor of Polynomial:
GCD）を計算する手法の一つとして二つの多項式$f,g$の$i$次シルベスター行列$S_i(f,g)$の唯一の$0$特異値に対応する右特異ベクトルからGCDの係数ベクトルを計算するものがある．この手法では$S_i(f,g)$の特異値分解を複数回行う必要があり計算量が莫大になってしまう．これを避ける手段として$S_i(f,g)$の$QR$分解から最小特異値を求めるものが存在する．本稿では$QR$から最小特異値を計算する方法と$S_i(f,g)$の$QR$分解から$S_{i-1}(f,g)$の$QR$分解を計算する方法について述べる．

最小特異値の計算

$S_i(f,g)$の二番目に小さい特異値$\sigma_2(S_i(f,g))$が$0$でなく，$QR=S_i(f,g)$が$S_i(f,g)$の$QR$分解であるとする．この時正しいサイズのランダムな初期値${\bf z}_0$に対して以下の手順

1. 前進代入によって$R^{H}{\bf y}_j = {\bf z}_{j-1}$をとく

2. 後退代入によって$R{\bf z}_j = {\bf y}_j$をとく

3. ${\bf z}_j$の正規化を行う

の繰り返しによってベクトルの系列${\bf z}_j, j=1,2,\ldots$は${\bf z}_\ast$に収束し，
$$\|S_k(f,g){\bf z}_\ast \| = \|R{\bf z}_\ast \| = \sigma_1(S_i(f,g)) \tag{1}$$
（ただし$\sigma_2(S_i(f,g))$は$(S_i(f,g)$の最小特異値）を満たす．またその収束率は
$$\| {\bf z}_j - {\bf z}_\ast \| \leq \left\{ \sigma _1(S_i(f,g)) / \sigma _2(S_i(f,g)) \right\} ^{2j} \| {\bf z}_0 - {\bf z}_\ast \|$$となるこの補題の証明は

S.V.Huffel ”Iterative algorithms for computing the singular subspace of a matrix associated with its smallest singular values” Linear Algebra and its Applications, Volumes 154-156, 1991, pp 675-709.

を参照されたい．
前進代入及び後退代入は$n\times m$行列に対して$O(nm)$で実行可能である．また実用上は３回程度の繰り返しによって${\bf z}_3 \approx {\bf z}_\ast $とできるため，一度$QR$分解が与えられれば$O(nm)$で最小特異値を計算することが可能である．これは$O(n^3)$の特異値分解と比べると大きな改善である．

ユニタリ変換による特異値の不変性

$A \in \mathbb{C}^{n \times m}$にユニタリ行列$P$と$Q^{H}$をそれぞれ左右からかけてもその特異値$\sigma_i(A), i=1,\ldots {\rm min}(n,m)$は変化しない

$A$の特異値分解$A = U\Sigma V^{H}$がある時，ユニタリ行列は積について閉じているため$(PU),\ (QV)^{H}$はどちらもユニタリ行列である．そのため$PAQ^{H}=(PU)\Sigma (QV)^{H}$は$PAQ^H$の特異値分解となっており，$A$と同じ特異値を持つ

列の入れ替え操作による特異値の不変性

右からかけることで列の入れ替えを行う行列を$P_i$と書く．するとこれはユニタリ行列であるため，この性質によってシルベスター行列$S_k(f,g)$の最小特異値を求める代わりに$S_k(f,g)P_k$の最小特異値を求めても同じであることがわかる．

シルベスター行列の更新

第一段階

$S_kP_k$の最下行にゼロベクトルを追加した時，以下の行列分解はQR分解となっている
$$\begin{bmatrix} S_kP_k\\ {\bf 0}^{\rm T}\\ \end{bmatrix} = \begin{bmatrix} Q_k & \\ & 1\\ \end{bmatrix} \begin{bmatrix} R_k\\ {\bf 0}^{\rm T}\\ \end{bmatrix}$$

第二段階

$\begin{bmatrix} S_kP_k\\ {\bf 0}^T\\ \end{bmatrix} =:[ {\bf a}_1, \ldots ,{\bf a}_{n+m-2k+2} ]$の第$n-k+2$列と行列の最も右に列ベクトル${\bf b}_1, {\bf b}_2$を追加すると$k-1$次シルベスター行列$S_{k-1}$を得る．
$$S_{k-1} = \left[ {\bf a}_1, \ldots , {\bf a}_{n-k+1},{\bf b}_1, {\bf a}_{n-k+2},\ldots,{\bf a}_{n+m-2k+2},{\bf b}_2 \right] \nonumber$$この行列に対して${\bf b}_1$が${\bf b}_2$の右隣に配置されるように置換行列$P_{k-1}$を右側からかけ，$$\begin{eqnarray} S_{k-1}P_{k-1} &=& [ {\bf a}_1, \ldots ,{\bf a}_{n+m-2k+2} ,{\bf b}_1,{\bf b}_2]\\ &=& \left[ \begin{pmatrix}S_{k-1}\\{\bf 0}^T\end{pmatrix}, {\bf b}_1, {\bf b}_2 \right]\end{eqnarray}$$を得る．上式の両辺に左側から
$$\begin{bmatrix} Q_k & \\ & 1\\ \end{bmatrix}^{\rm H} \left(= \begin{bmatrix} Q_k^{\rm H} & \\ & 1\\ \end{bmatrix} \right) \nonumber$$ をかけ ${\bf u}_1 = \begin{bmatrix} Q_k^{\rm H} & \\ & 1\\ \end{bmatrix}{\bf b}_1,\ {\bf u}_2 = \begin{bmatrix} Q_k^{\rm H} & \\ & 1\\ \end{bmatrix}{\bf b}_2$とおくと $$\begin{eqnarray} \begin{bmatrix} Q_k^{\rm H} & \\ & 1\\ \end{bmatrix}S_{k-1}P_{k-1} &=& \left[ \begin{pmatrix} Q_k^{\rm H} & \\ & 1\\ \end{pmatrix}\begin{pmatrix}S_{k-1}\\{\bf 0}^{\rm T}\end{pmatrix}, {\bf u}_1, {\bf u}_2 \right]\\ &=& \left[ \begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix}, {\bf u}_1, {\bf u}_2 \right] \tag{2}\end{eqnarray}$$である．$S_{k-1}P_{k-1}$の$QR$分解としては，$\left[ \begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix}, {\bf u}_1, {\bf u}_2 \right]$の対角成分より下が$0$になるように${\bf u}_1,{\bf u}_2$をユニタリ行列で変形すれば良い．

ギブンズ回転

行列またはベクトルのある成分の零化はギブンズ回転によって行われる．

ギブンズ回転とは，正方行列 $$G(i,j,\theta) = \begin{pmatrix} 1 \\ & \ddots \\ & & \cos \theta & \dots & \sin \theta \\ & & \vdots & & \vdots \\ & & -\sin \theta & \dots & \cos \theta \\ & & & & & \ddots \\ & & & & & & 1 \end{pmatrix} \nonumber$$をベクトルの右側からかけることで鏡像回転を行うハウスホルダー変換の一種であり，明らかに
$$G^{T}(i,j,\theta) = G^{-1}(i,j,\theta) = G(i,j,-\theta)$$である．
$$\cos \theta = \frac{u_i}{\sqrt{u_i^2 + u_j^2}},\ \sin \theta = \frac{u_j}{\sqrt{u_i^2 + u_j^2}} \tag{3}$$とすると
$${\bf v} = G(i,j,\theta){\bf u}$$は $$v_l = \begin{cases} \sqrt{u_i^2 + u_j^2} & l = i \\ 0 & l = j\\ u_l & l \neq i,j \end{cases}$$となり，${\bf u}$の第$j$要素を利用して第$i$要素を零化できる．

QR分解を完了する

ギブンズ回転を利用して（2）式の三角化を引き続き行うと$\begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix}$は既に上三角行列となっているため，$S_{k-1}P_{k-1}$の$QR$分解は${\bf u}_1$,${\bf u}_2 $の対角以下の成分が$0$になるように右側からギブンズ行列をかければ良い．ベクトル${\bf v }$の第$j$要素$u_j$を用いて第$j+1$要素を消去する行列を（3）の$\theta$に対して$G_{{\bf u},j}^{T}=G^{T}(j,j+1,\theta)$とかくことにする．このギブンズ回転を複数回行うことで${\bf u}_1$,${\bf u}_2 $は
$$\begin{eqnarray} {\bf u}'_1 &=&G_{{\bf u}_1,l}^{T}G_{{\bf u}_1,l+1}^{T}\cdots G_{{\bf u}_1,n}^{T}{\bf u}_1,\nonumber\\ {\bf u}'_2 &=& G_{{\bf u}_1,l}^{T}G_{{\bf u}_1,l+1}^{T}\cdots G_{{\bf u}_1,n}^{T}{\bf u}_2,\nonumber\\ {\bf u}'_1 &=&G_{{\bf u}_2,l+1}^{T}G_{{\bf u}_2,l+2}^{T}\cdots G_{{\bf u}_2,n}^{T}{\bf u}'_1,\nonumber\\ {\bf u}''_2 &=& G_{{\bf u}_2,l+1}^{T}G_{{\bf u}_2,l+2}^{T}\cdots G_{{\bf u}_2,n}^{T}{\bf u}'_2\nonumber\end{eqnarray}$$と余分な成分の零化を実行できる．なお$\begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix}$の第$l$成分以下は既に全て$0$であるためこれらのギブンズ回転を行っても変化しない．つまり
$$\begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix} = G_{{\bf u}_2,l+1}^{T}G_{{\bf u}_2,l+2}^{T}\cdots G_{{\bf u}_2,n}^{T}\begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix} \nonumber$$である．$G_{{\bf u}_2,l+1}^{T}G_{{\bf u}_2,l+2}^{T}\cdots G_{{\bf u}_2,n}^{T}G_{{\bf u}_1,l}^{T}G_{{\bf u}_1,l+1}^{T}\cdots G_{{\bf u}_1,n}^{T}$を（2）の両辺に右側からかけることで
$$G_{{\bf u}_2,l+1}^{T}G_{{\bf u}_2,l+2}^{T}\cdots G_{{\bf u}_2,n}^{T}G_{{\bf u}_1,l}^{T}G_{{\bf u}_1,l+1}^{T}\cdots G_{{\bf u}_1,n}^{T} \begin{bmatrix} Q_k^{\rm H} & \\ & 1\\ \end{bmatrix}S_{k-1}P_{k-1} = \left[ \begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix}, {\bf u}'_1, {\bf u}''_2 \right] \nonumber$$を得る．式変形を行うと
$$S_{k-1}P_{k-1} = \underbrace{\begin{bmatrix} Q_k & \\ & 1\\ \end{bmatrix} G_{{\bf u}_1,n} \cdots G_{{\bf u}_1,l}G_{{\bf u}_2,n} \cdots G_{{\bf u}_2,l+1}} _{Q_{k-1}}\ \underbrace{\left[ \begin{pmatrix}R_{k}\\{\bf 0}^{\rm T}\end{pmatrix}, {\bf u}'_1, {\bf u}''_2 \right]} _{R_{k-1}}$$が$S_{k-1}P_{k-1}$の$QR$分解となっている．