はじめに

行列の低ランク近似といえば特異値分解、ということになっている。Eckart-Youngの定理からそれがフロベニウスノルムおよびスペクトルノルムの意味で最良の近似であることが保証されているので、よく使われるのは理の当然だ。ただし、いかんせん計算コストが大きく、大規模行列には適用しにくいという側面もある。

そこで、QR分解による低ランク近似も色々研究されているようだ。同じ条件下では特異値分解より若干精度は劣るものの、その分計算コストも低くなっており、大規模行列に適した手法ということらしい。

以下、自分用のメモとしてQR分解から派生する低ランク近似のアルゴリズムについて簡単にまとめる。

QR分解

行列$\mathbf{A} \in \mathbb{R}^{n \times m}$のQR分解とは以下の形の行列分解である（$n > m$とする）。

\mathbf{A} = \mathbf{Q} \mathbf{R}

ここで、$\mathbf{Q} \in \mathbb{R}^{n \times m}$は直交行列、$\mathbf{R} \in \mathbb{R}^{m \times m}$は上三角行列である。特に$\mathbf{Q}$のサイズを$n \times m$とする定式化はthin QR分解とかeconomic QR分解といわれ、数値計算用の各種ライブラリもおおむねこれに対応している。QR分解はGram-Schmidt直交化と数学的に等価であり、ベクトルを順次直交化することで$\mathbf{Q}$と$\mathbf{R}$を求めることができる。ただし、一般的にはより直交性の精度に優れるHouseholder変換が用いられる。

特に$n \gg m$となる縦長行列の場合、前処理としてQR分解を行うことによって各種の行列計算のコストを削減することができる。たとえば、最初にQR分解を行い、その後サイズの小さい$\mathbf{R}$の特異値分解$\mathbf{R} = \mathbf{U}'\mathbf{\Sigma}\mathbf{V}^{\top}$を求めれば、$\mathbf{A}$の特異値分解は$\mathbf{A} = \mathbf{Q}\mathbf{U}'\mathbf{\Sigma}\mathbf{V}^{\top}$で得られる（$\mathbf{U} = \mathbf{Q}\mathbf{U}'$）。QR分解の計算コストは特異値分解よりも小さいため、こうすることで$\mathbf{A}$の特異値分解を直接計算するよりもトータルのコストを小さくできる。

列ピボット選択付きQR分解

特異値分解は特異値と特異ベクトルが元の行列に対する寄与度の大きな順に並んでいるため、より重要な上位$k$個の特異値・特異ベクトルのみを残して他の成分を打ち切ることにより、簡単に低ランク近似を行うことができる。一方、通常のQR分解では寄与度の大きな成分の情報は陽にはわからない。そこで、次のような列ピボット選択付きQR分解（以下、CPQR分解¹）が用いられる。

\begin{align}
\mathbf{A} \mathbf{P} &= \mathbf{Q} \mathbf{R} \\
&= 
\begin{bmatrix}
\mathbf{Q}_{1} & \mathbf{Q}_{2}
\end{bmatrix}
\begin{bmatrix}
\mathbf{R}_{11} & \mathbf{R}_{12} \\
 & \mathbf{R}_{22} 
\end{bmatrix}
\end{align}

ここで、$\mathbf{P} \in \mathbb{R}^{m \times m}$は列の置換を表す置換行列であり、

\left| r_{11} \right| \geq \left| r_{22} \right| \geq \cdots \geq \left| r_{mm} \right|

のように上三角行列の対角成分の絶対値が降順になるような順序に対応する。なお、この場合、行列のサイズは$\mathbf{Q} \in \mathbb{R}^{n \times n}$、$\mathbf{R} \in \mathbb{R}^{n \times m}$である。また、分解後のブロック行列のサイズはそれぞれ$\mathbf{Q}_{1} \in \mathbb{R}^{n \times k}$、$\mathbf{Q}_{2} \in \mathbb{R}^{n \times (m-k)}$、$\mathbf{R}_{11} \in \mathbb{R}^{k \times k}$、$\mathbf{R}_{12} \in \mathbb{R}^{k \times (m-k)}$、$\mathbf{R}_{11} \in \mathbb{R}^{(n-k) \times (m-k)}$である。

上式を展開すると、

\mathbf{A}
= \mathbf{Q}_{1} 
\begin{bmatrix}
\mathbf{R}_{11} & \mathbf{R}_{12} 
\end{bmatrix}
\mathbf{P}^{\top}
+ \mathbf{Q}_{2} 
\begin{bmatrix}
\mathbf{0} & \mathbf{R}_{22} 
\end{bmatrix}
\mathbf{P}^{\top}

となるので、$\left\| \mathbf{R}_{22} \right\|_{F} \approx 0$であれば、

\mathbf{A}
\approx \mathbf{Q}_{1} 
\begin{bmatrix}
\mathbf{R}_{11} & \mathbf{R}_{12} 
\end{bmatrix}
\mathbf{P}^{\top}

がよい近似となる。$\left\| \mathbf{R}_{22} \right\|_{F} \approx 0$ということは、元の行列の数値的なランクが$k$であることを意味するので、CPQR分解によって行列の数値的なランクを求めることもできる。実際、特異値分解でも数値的なランクを計算でき、そのようなアルゴリズムはRank Revealingアルゴリズムと呼ばれる。

Rank Revealing QR分解

特に、以下の条件を満たすような置換行列$\mathbf{P}$が存在するとき、CPQR分解はRank Revealing QR（RRQR）分解と呼ばれる。

$\sigma_{\min}\left( \mathbf{R}_{11} \right) \approx \sigma_{k}\left( \mathbf{A} \right)$または$\sigma_{\max}\left( \mathbf{R}_{22} \right) \approx \sigma_{k+1}\left( \mathbf{A} \right)$
$\sigma_{\min}\left( \mathbf{R}_{11} \right) \gg \sigma_{\max}\left( \mathbf{R}_{22} \right) = O\left( \left\| \mathbf{R}_{22} \right\|_{2} \right)$

ここで、$\operatorname*{cond}\left( \mathbf{A} \right) = \left. \sigma_{\max}\left( \mathbf{A} \right) \right/ \sigma_{\min}\left( \mathbf{A} \right)$を行列$\mathbf{A}$の条件数とすると、最初の式は$\operatorname*{cond}\left( \mathbf{A} \right) \approx \operatorname*{cond}\left( \mathbf{R}_{11} \right)$とも表現できる。つまり、$\mathbf{R}_{11}$の上位$k$個の特異値が元の行列のそれと十分近ければよい近似になるため、Rank Revealingアルゴリズムのよさは特異値の近似精度によって判断されるということだ。
2番目の条件は数値的なランクの定義であり、マシンイプシロン程度の小さな値をしきい値として、特異値の大きさにギャップがあるところでランクを決めることを意味する。

なお、数学的にはRRQR分解は失敗する可能性があるらしい（Kahan行列という反例がある）が、相当変なデータでもなければ実際上はうまくいくので、あまり気にしなくてもよい。
また、より近似精度を高めたStrong RRQR分解というものもあるようだが、アルゴリズムがよくわからなかったのでここでは省略する。²

なお、CPQR分解を元のサイズの大きな行列に対して適用すると計算コストが大きい（列ノルムが大きい順に並べ替える処理があるため）。そのため、

\begin{align}
\mathbf{A} &= \mathbf{Q} \mathbf{R} \\
\mathbf{R} \mathbf{P} &= \mathbf{Q}' \mathbf{R}' \\
\therefore \mathbf{A}\mathbf{P} &= \left( \mathbf{Q} \mathbf{Q}' \right) \mathbf{R}'
\end{align}

のように、前処理として最初に通常のQR分解を行い、得られた上三角行列に対してCPQR分解を行うことで全体の計算コストを小さくできる。

補間分解

CPQR分解と関係が深い行列分解として次の補間分解（Interpolative Decomposition, ID）がある。

\mathbf{A} \approx \mathbf{A}_{(:,J)} \mathbf{Z}

ここで、$\mathbf{A}_{(:,J)} \in \mathbb{R}^{n \times k}$は元の行列から添字集合$J = \left[ j_{1}, j_{2},\cdots, j_{k} \right]$に対応する列のみを取り出した行列である。$\mathbf{Z}\in \mathbb{R}^{k \times m}$は$\mathbf{Z}_{(:,J)} = \mathbf{I}_{k}$を満たす行列である。
上記は列を取り出した場合の例だが、行を取り出す場合や列と行を両方を取り出す場合についても同様に考えることができる（両側補間分解はCUR分解ともいわれる）。

補間分解は元の行列の成分をそのまま使用しているため、寄与の大きい重要なデータを抽出するためなどに利用できる。また、元の行列の疎性や非負性を引き継ぐため、分解によってそれらが失われないことが保証される。

上述の通り、行列$\mathbf{A}$をCPQR分解すると、

\begin{align}
\mathbf{A} \mathbf{P} &= 
\begin{bmatrix}
\mathbf{Q}_{1} & \mathbf{Q}_{2}
\end{bmatrix}
\begin{bmatrix}
\mathbf{R}_{11} & \mathbf{R}_{12} \\
\mathbf{0} & \mathbf{R}_{22} 
\end{bmatrix}
\\
\Rightarrow \mathbf{A}
&= \mathbf{Q}_{1} 
\begin{bmatrix}
\mathbf{R}_{11} & \mathbf{R}_{12} 
\end{bmatrix}
\mathbf{P}^{\top}
+ \mathbf{Q}_{2} 
\begin{bmatrix}
\mathbf{0} & \mathbf{R}_{22} 
\end{bmatrix}
\mathbf{P}^{\top}
\end{align}

となる。右辺第一項について見てみると、

\begin{align}
\mathbf{Q}_{1} 
\begin{bmatrix}
\mathbf{R}_{11} & \mathbf{R}_{12} 
\end{bmatrix}
\mathbf{P}^{\top}
=
\mathbf{Q}_{1} \mathbf{R}_{11}
\begin{bmatrix}
\mathbf{I}_{k} & \mathbf{R}_{11}^{-1} \mathbf{R}_{12} 
\end{bmatrix}
\mathbf{P}^{\top}
\end{align}

であるから、

\begin{align}
\mathbf{A}_{(:,J)} &= \mathbf{Q}_{1} \mathbf{R}_{11} \\
\mathbf{Z} &= 
\begin{bmatrix}
\mathbf{I}_{k} & \mathbf{R}_{11}^{-1} \mathbf{R}_{12} 
\end{bmatrix}
\mathbf{P}^{\top}
\end{align}

が成り立つ。また、

\begin{align}
\mathbf{A} - \mathbf{A}_{(:,J)} \mathbf{Z} &= 
\mathbf{Q}_{2} 
\begin{bmatrix}
\mathbf{0} & \mathbf{R}_{22} 
\end{bmatrix}
\mathbf{P}^{\top}
\end{align}

であり、上式の右辺が残差に相当するため、$\left\| \mathbf{R}_{22} \right\|_{F}$が十分小さければ補間分解はよい近似となる。
実際には$\mathbf{A}_{(:,J)}$はCPQR分解から求まる置換ベクトル$P$を用いて$J = P_{(1:k)}$とする。また、$\mathbf{R}_{11}$が正則でない場合も考慮して、$\mathbf{Z}$の計算には逆行列の代わりに擬似逆行列を使うこともできる。

QLP分解

CPQR分解から派生する別の手法としてQLP分解がある。CPQR分解の結果として求まる上三角行列の対角成分の絶対値（R値）は特異値の近似になっているが、精度はそこまで高くない。そこで、近似精度を高めるために、行についても入れ替えを行うことを考える。すなわち、

\begin{align}
\mathbf{A} \mathbf{P}_{1} &= \mathbf{U} \mathbf{R} \\
\mathbf{R}^{\top} \mathbf{P}_{2} &= \mathbf{V} \mathbf{L}^{\top} \\
\end{align}

というように、上三角行列の転置行列に対して再度CPQR分解を計算する。ここで、$\mathbf{L}$は下三角行列である。その結果、次式を得る。

\begin{align}
\mathbf{R} &= \mathbf{P}_{2} \mathbf{L} \mathbf{V}^{\top} \\
\therefore \mathbf{A} &= \mathbf{U} \mathbf{P}_{2} \mathbf{L} \mathbf{V}^{\top} \mathbf{P}_{1}^{\top}
\end{align}

ここで、$\mathbf{U} \mathbf{P}_{2}$、$\mathbf{V}^{\top} \mathbf{P}_{1}^{\top}$は左右の特異ベクトルの近似であり、$\mathbf{L}$の対角成分の絶対値（L値）が特異値の近似となる。
このような特異値分解に似た形の近似的な分解は一般にUTVフレームワークと呼ばれ、QLP分解もその一種である。

特異値の近似の比較

計算例として、CPQR分解のR値とQLP分解のL値が特異値をどれだけ近似できるか試してみる。計算はJulia 1.4で行った。なお、データ行列は以前の記事でも使った円柱周りの流れのデータである。

まず、CPQR分解は組み込みのqr関数を使って簡単に計算できる。

using LinearAlgebra

# D:データ行列

F = qr(D) # QR分解
R = Matrix(F.R) # Rを行列型に変換
F = qr(R, Val(true)) # CPQR分解
R = Matrix(F.R) # Rを通常の行列に変換

内部的にはLAPACKのgeqp3ルーチンが呼ばれており、Householder変換によって直交行列が計算されるようだ。結果はメモリ節約のために圧縮された形式で保存されており、通常の行列データに変換する場合はMatrix(F.R)のように変換する必要がある。
CPQR分解にする場合はVal(true)オプションを与えると列ピボッティングが有効になる（結果の置換ベクトルはF.pで取り出せる）。

QLP分解を行う場合は上記に続けて再度CPQR分解を行えばよい。

F = qr(R', Val(true)) # CPQR分解（QLP分解）
L = Matrix(F.R)' # Lを通常の行列に変換

最後にお目当てのR値とL値を求めておく（正解としての特異値も忘れずに）。

R_val = abs.(diag(R)) # R値
L_val = abs.(diag(L)) # L値
S = svd(R).S # 特異値

結果を適当にプロットすると、以下のようになった。

using Plots

plot([S[1:100], R_val[1:100], L_val[1:100]], seriestype=:scatter, yaxis=:log, label=["singular value" "R val" "L val"])

図より、CPQR分解のR値によってデータの寄与度の大きさを判断できるようになったことがわかる。一方、R値は特異値に比べて減衰が遅く、またところどころに値のギャップが見られるので、打ち切るランクの設定によっては低ランク近似の精度が悪そうだ。また、2回CPQR分解を行ったQLP分解のL値はより特異値に近づいていることがわかる。

計算速度的にはもちろん、CPQR分解やQLP分解の方が特異値分解よりも速い。ただし、実際上は前処理によって行列のサイズを小さくしてから各種の計算を行うのが常であるため、律速要因は前処理の手法になると思われる（今の場合はQR分解）。そのため、行列のサイズと前処理手法との組み合わせ、精度と計算速度のバランスなどを考えながら最適な行列分解の手法を選択することが重要と思われる。

参考文献

A Note on QR-Based Model Reduction: Algorithm, Software, and Gravitational Wave Applications
The Interpolative Decomposition (ID)
Rank revealing factorizations, and low rank
approximations
Rank Revealing QR分解による悪条件線形方程式の解法
Golub & van Loan, Matrix Computations(4th Edition)

英語ではColumn Pivoted QR decompositionとかQR decomposition with Column Pivotingとかいうらしい。 ↩
CPQR分解のアルゴリズムは貪欲法に基づいており、最適な置換行列が得られない場合があるため、精度を高めるために置換をやり直す、という話らしい（多分）。 ↩