行列の平方根の数値計算

この記事は 数値計算 Advent Calendar 2018 7日目の記事です

本記事では，「行列の平方根」と3種類の計算方法について紹介します．
また，Juliaによる実装も行います．
「行列の平方根」は聞き慣れない言葉ですが，定義自体はシンプルです．
楽しんでいただけると幸いです．

1. はじめに

行列の平方根について

正方行列 $A$ の平方根 $A^{1/2}$ とは， $X^2=A$ を満たす $X$ のことです¹．

たとえば，$A = \begin{bmatrix}4&0\\0&9\end{bmatrix}$ のとき， $X = \begin{bmatrix}2&0\\0&3\end{bmatrix}$ は $X^2=A$ の解なので，$X$ は $A$の平方根です．
お気づきかもしれませんが，$A=ZD Z^{-1}$ と対角化できるとき，$A^{1/2}=ZD^{1/2}Z^{-1}$ です．

行列平方根がメインの計算になることは少ないですが，半正定値計画問題や境界値問題，医用イメージング²など幅広い分野に現れるため，こんな計算があるんだ！と知っておくとオトクな気がします．

計算方法の簡単な紹介

いくつかのプログラミング言語では行列平方根の計算が実装されています．
Pythonではscipy.linalg.sqrtm，Juliaではsqrt，MatlabやOctaveではsqrtmを使えます．

自力で組む場合は，以下の3種類の計算方法がおすすめです．
それぞれに長所があるので，計算用途に合わせて選びましょう．

• 行列分解 (Schur法)
  • 他の手法よりも（おそらく）高速
  • JuliaやPythonは本記事のアルゴリズムの改良版を使っているはず
• Newton法
  • 行列積と逆行列積のみで計算ができる．
  • 多倍長演算の環境などで行列分解ができなかったとしても，平方根が計算できる
• 数値積分
  • 逆行列を足していくだけなので並列計算しやすそう

本記事では上記の3手法を紹介し，プログラム例を示します．

2. 行列分解（Schur法）

はじめの例のように，$A$ が対角可能であれば，対角化によって $A^{1/2}$ を計算することができます．
実際，Juliaでは $A$ が対称行列ならば対角化によって計算を行っています (ソース)．

「対称行列ならば」と言ったのは，非対称行列の対角化は数値的に不安定になりやすく避けられているためです．
ではどうするのかというと，数値的に安定なSchur分解を用います．

Aを（複素）Schur分解すると $A=QTQ^\mathsf{H}$ のようにユニタリ行列 $Q$ と三角行列 $T$ の積に分解されます．$Q^\mathsf{H}$ は $Q$ の共役転置です．実は，三角行列の平方根は直接的に計算することが可能です．そこで，三角行列の平方根 $U=T^{1/2}$ を求めておいて， $A^{1/2}=QT^{1/2}Q^{\mathsf{H}}$ を計算します．

Schur法のアルゴリズム³は以下の通りです．

• 入力: $A$ ($n\times n$ 行列)，出力: $X \approx A^{1/2}$
• $A$ を複素Schur分解する: $A=QTQ^\mathsf{H}$
• $n\times n$ の行列 $U$ を作り，初期化する
• $u_{i,i} = \sqrt{t_{i,i}}~~(i=1,\dots,n)$
• $u_{i,j}=\cfrac{t_{i,j}-\sum_{k=i+1}^{j-1} u_{i,k}u_{k,j}}{u_{i,i}+u_{j,j}}~~(j=2,\dots,n,~i=j-1,\dots,1)$
• $X = QUQ^\mathsf{H}$

3. Newton法

Newton法は非線形方程式の解法です．スカラー方程式 $f(x) = 0$ に対するNewton法の漸化式は

$$x_{k+1}=x_k - \frac{f(x_k)}{f'(x_k)}$$

で与えられます．例えば，$f(x) = x^2 - a$ を代入して整理すると，漸化式は

$$x_{k+1} = \frac{1}{2}\left(x_k+\frac{a}{x_k}\right)$$

となり，適当な初期値のもとで $x_k$ は $\sqrt{a}$ に収束します．実際の計算では初期値を $a$ とすれば良いでしょう．

行列の平方根は $X^2-A=O$ の解なので，あたかも $X$ をスカラーと思ってNewton法を適用します．初期値が$A$のとき，行列平方根に対するNewton法は

$$
\begin{cases}
X_0 = A,\
X_{k+1} = \cfrac{1}{2}\left(X_k + X_k^{-1}A\right)
\end{cases}
$$

となり，数学的には $X_k$ は $A^{1/2}$ に収束します⁴．しかし，上記の漸化式(1)は数値的に不安定であることが知られているため，(1)と等価な
$$
\begin{cases}\tag{2}
X_0 = A,~ ~Y_0=I,~ ~M_0=A,\
X_{k+1} = \cfrac{1}{2},X_k\left(I+M_k^{-1}\right),\
Y_{k+1} = \cfrac{1}{2},Y_k\left(I+M_k^{-1}\right),\
M_{k+1} = \cfrac{1}{2}\left(I+\cfrac{M_k+M_k^{-1}}{2}\right)
\end{cases}
$$
を用いることをおすすめします．(2)は積型DB反復式と呼ばれています⁵．
式は複雑になりましたが，$\lim_{k\to\infty}X_k= A^{1/2}$ なので，適当な回数だけ解を更新して，$X_k$を取り出しましょう．

4. 数値積分近似

行列の平方根は区間 $[-1,1]$ 上の積分を用いて表せることが知られています⁶:

$$A^{1/2}=\frac{4}{\pi}A\int_{-1}^1 \left[(1+x)^2I+(1-x)^2A\right]^{-1},\mathrm{d}x.$$

被積分関数の中身が行列なのでギョッとするかもしれません．実際は，被積分関数の入力はスカラーなので，$x$ を-1から1まで動かしながら出てくる行列を足していくだけで $A^{1/2}$ が計算できます．
区間が $[-1,1]$ なので，Gauss−Legendre求積などで計算すれば良いでしょう．下の数値例ではGauss−Legendre求積の積分点と重みの計算にFastGaussQuadrature.jlを用いました．

5. 実装例

$A=\begin{bmatrix}5&4&1\\4&6&4\\1&4&5\end{bmatrix}$ に対して平方根 $A^{1/2}=\begin{bmatrix}2&1&0\\1&2&1\\0&1&2\end{bmatrix}$ を求めます．

using LinearAlgebra
using FastGaussQuadrature

# Schur法
function sqrt_schur(A)
    isreal_A = isreal(A)
    T, Q = schur(convert(Array{ComplexF64,2}, A))
    n = size(A)[1]
    U = zeros(Complex, n,n)
    for i = 1:n
        U[i,i] = sqrt(T[i,i])
    end
    for j = 2:n
        for i = j-1:-1:1
            tmp = 0.0
            for k=i+1:j-1
                tmp += U[i,k]*U[k,j]
            end
            U[i,j] = (T[i,j] - tmp) / (U[i,i]+U[j,j])
        end
    end
    return isreal_A ? real.(Q*U*Q') : Q*U*Q'
end

# Newton法
function sqrt_newton(A, TOL, MAXITER)
    n = size(A)[1]
    X = copy(A)
    M = copy(A)
    Y = diagm(0=>ones(n))
    k = 0
    relres = norm(X^2-A) / norm(A)
    while k < MAXITER && relres > TOL
        k += 1
        M_inv = inv(M)
        X = 0.5 * X * (I + M_inv)
        Y = 0.5 * Y * (I + M_inv)
        M = 0.5 * (I + 0.5*(M + M_inv))
        relres = norm(X^2-A) / norm(A)
    end
    return X
end

# 積分近似
function sqrt_integral(A, m)
    F(x) = inv((1+x)^2*I + (1-x)^2*A)
    x, w = gausslegendre(m)
    n = size(A)[1]
    G = zeros(n,n)
    for i = 1:m
        G += w[i] * F(x[i])
    end
    return 4/pi * A * G
end

A = [5 4 1; 4 6 4; 1 4 5]
Ans = [2 1 0; 1 2 1; 0 1 2]

# Juliaのsqrt関数
X = sqrt(A)
err = norm(Ans - X)/norm(Ans)

# Schur法
X_schur = sqrt_schur(A)
err_schur = norm(Ans - X_schur)/norm(Ans)

# Newton法
X_newton = sqrt_newton(A, 1e-12, 30)
err_newton = norm(Ans - X_newton)/norm(Ans)

# 数値積分
X_int = sqrt_integral(A, 32)
err_int = norm(Ans - X_int)/norm(Ans)

println("sqrt関数の誤差: $(err)")
println("Schur法の誤差: $(err_schur)")
println("Newton法の誤差: $(err_newton)")
println("数値積分近似の誤差: $(err_int)")
# sqrt関数の誤差: 2.310557899261868e-15
# Schur法の誤差: 8.860073618950564e-16
# Newton法の誤差: 1.9603600581608083e-16
# 数値積分近似の誤差: 6.010860650248393e-16

参考文献

• [1] Higham, Nicholas J. Functions of Matrices: Theory and Computation. Vol. 104. Siam, 2008.
• [2] Cardoso, João R. "Computation of the matrix $p$th root and its Fréchet derivative by integrals" Electron. Trans. Numer. Anal 39 (2012): 414-436.

1. $1$ と $-1$ がともに $1$ の平方根であるように，$X^2=A$ を満たす行列 $X$ は複数存在し得ます．本記事では，$A$ の平方根の中でも主値に相当するもの (principal square root) を考えます．すなわち，方程式の解の中でも，全ての固有値が複素平面の右半平面 $\{z\in\mathbb{C}:|\mathrm{Re}(z)|>0\}$ に属する解を考えます．なお，行列の主平方根は，$A$の全ての固有値が $\mathbb{C}\setminus (-\infty,0]$ に属せば一意に存在します [Thrm. 1.29, 1]． ↩

2. 行列の対数関数の計算の応用として医用イメージングがあります．MatlabやPythonの行列対数関数の計算では（非対称行列の時は）行列平方根の計算を用いるアルゴリズムを利用しています． ↩

3. [Algorithm 6.3, 1] など ↩

4. [Eq. 6.12, 1] ↩

5. [Eq. 6.17, 1] ↩

6. [Eq. 2.6, 2] など ↩