📂行列代数

特異値分解を用いたティホノフ正則化の計算量削減トリック

定理

リッジ回帰の最適解: $$ L \left( \beta \right) = \left\| Y - X \beta \right\|_{2}^{2} + \lambda \left\| \beta \right\|_{2}^{2} $$ $\lambda$が定数として与えられているとき、リッジ回帰の目的関数$L$を上のように表そう。リッジ回帰の最適解$\hat{\beta} = \argmin_{\beta} L \left( \beta \right)$は次の通りである。 $$ \hat{\beta} = \left( X^{\top} X + \lambda I \right)^{-1} X^{\top} Y $$ ここで$A^{\top}$は$A$の転置行列、$I$は単位行列、$A^{-1}$は$A$の逆行列である。

行列$A \in \mathbb{R}^{m \times n}$とベクトル$\mathbf{b} \in \mathbb{R}^{m}$が与えられたとき、$A \mathbf{x} = \mathbf{b}$の最小二乗解を求める際に、$\lambda \ge 0$に対するティホノフ正則化^{Tikhonov regularization}を適用すると、次の最適化問題を解くことになる。 $$ \argmin_{\mathbf{x}} \left\| A \mathbf{x} - \mathbf{b} \right\|^2 + \lambda \left\| \mathbf{x} \right\|^2 $$ 通常この問題の解はクローズドフォームである$\mathbf{x} = \left( A^{\top} A + \lambda I \right)^{-1} A^{\top} \mathbf{b}$として知られているが、$A$の特異値分解$A = U \Sigma V^{\top}$を用いれば、行列分解を一度だけ行って$\lambda$を変えながら計算できる次の公式でも表すことができる。 $$ \mathbf{x}(\lambda) = V \left( \Sigma^{2} + \lambda I \right)^{-1} \Sigma U^{\top} \mathbf{b} $$

説明

さまざまな行列分解による最小二乗法の実装について速度を比較し、安定性を比較してみると、次のように非常に明確なトレードオフが見える。

方法	速度	安定性
LU	速い	不安定
QR	普通	安定的
SVD	遅い	非常に安定的

実際、QRとSVDそのものの計算量だけを見るとSVDは遅いが、たとえば行列が非常に大きく、リッジ回帰のために$\lambda$を変えながら計算を繰り返して最適な値を探す場合には、SVDによるトリックを通じて計算量を減らすことができる。提示された公式において$\left( \Sigma^{2} + \lambda I \right)$は対角行列なので逆行列を求めるのは非常に容易であり、$V$と$U^{\top} \mathbf{b}$は一度だけ計算すればすべての$\lambda$に対して再利用可能である。

このようなトリックは単に計算量を削減するだけでなく、SVDによる安定性まで確保できるので、莫大な規模と悪条件の危険性が高まるデータに近づくほど有用な方法となりうる。

もちろん上で述べた問題がない場合、たいていはQR分解に勝てないか勝つ必要がない。たとえば1秒かかっていたタスクが0.5秒に減ったからといって、わざわざSVDベースのリッジ回帰を実装するのはあまりにも過剰であり、検証されていないコードを自分で書く時点で実装の完成度も疑いを招く。SVDを適用すべき状況かどうかを合理的に判断し、実装に成功したとしても実際に性能向上があったかについての検証が必ず必要である。

証明

$A = U \Sigma V^{\top}$を$A$の特異値分解 (SVD)としよう。ここで$U$と$V$は直交行列であり、$\Sigma$は特異値を含んでいる対角行列である。$V$が直交行列なので$\mathbf{y} = V^{\top} \mathbf{x}$と置換すると次を得る。 $$ \lambda \left\| \mathbf{x} \right\|^2 = \lambda \left\| V \mathbf{y} \right\|^2 = \lambda \mathbf{y}^{\top} V^{\top} V \mathbf{y} = \mathbf{y}^{\top} \lambda I \mathbf{y} $$ そして $$ \begin{align*} \left\| A \mathbf{x} - \mathbf{b} \right\|^2 =& \left\| U \Sigma V^{\top} V \mathbf{y} - \mathbf{b} \right\|^2 \\ =& \left\| U \Sigma \mathbf{y} - \mathbf{b} \right\|^2 \\ =& \left( U \Sigma \mathbf{y} - \mathbf{b} \right)^{\top} \left( U \Sigma \mathbf{y} - \mathbf{b} \right) \\ =& \mathbf{y}^{\top} \Sigma^{\top} U^{\top} U \Sigma \mathbf{y} - 2 \mathbf{y}^{\top} \Sigma U^{\top} \mathbf{b} + \left\| \mathbf{b} \right\|^{2} \\ =& \mathbf{y}^{\top} \Sigma^{2} \mathbf{y} - 2 \mathbf{y}^{\top} \Sigma U^{\top} \mathbf{b} + \left\| \mathbf{b} \right\|^{2} \end{align*} $$ であるので、$\mathbf{y}$に対する目的関数とそのグラディエントは次のように表すことができる。 $$ \begin{align*} J \left( \mathbf{y} \right) =& \mathbf{y}^{\top} \left( \Sigma^{2} + \lambda I \right) \mathbf{y} - 2 \mathbf{y}^{\top} \Sigma U^{\top} \mathbf{b} + \left\| \mathbf{b} \right\|^{2} \\ \partial J \left( \mathbf{y} \right) / \partial \mathbf{y} =& 2 \left( \Sigma^{2} + \lambda I \right) \mathbf{y} - 2 \Sigma U^{\top} \mathbf{b} . \end{align*} $$ 最後に、与えられた$\lambda$に対する最適解$\mathbf{x}(\lambda)$は次の通りである。 $$ \mathbf{x}(\lambda) = V \left( \Sigma^{2} + \lambda I \right)^{-1} \Sigma U^{\top} \mathbf{b} $$

■

コード

次は筆者が研究に実際に用いているSVDベースのティホノフ正則化である。

using LinearAlgebra

function tikhonov_λ(A, b)
    U, s, V = svd(A)
    Uᵀb = U'b
    return λ -> V * (Diagonal(s ./ (s.^2 .+ λ)) * Uᵀb)
end

筆者の場合、Uᵀb = U'bを書き忘れて無名関数にそのまま入れ、毎回$\lambda$を変えながら$U^{\top} \mathbf{b}$を計算してしまい、計算量が基本実装(QR分解)に劣るというミスをしたことがある。数式だけを信じず、実際のベンチマークを必ずやってみよう。