大阪大学医学部 Python会

Now is better than never.

行列方程式のクロネッカー積による解法

2021-02-09(Tue) - Posted by 山本 in 技術ブログ    tag:Python

Contents

    運動制御モデルの論文を読んでいた際に，行列方程式\(AXB + CXD = E\)の解法を忘却していたので，メモとして記事を書いた．結論から言えばクロネッカー積を用いて解くことができる．以下ではクロネッカー積の定義について触れ，簡単な行列方程式\(AXB=C\)に対してクロネッカー積を用いた解法を用いた後，\(AXB + CXD = E\)の解き方について説明する．

    なお，これに基づいて筆者が執筆した記事が無限時間最適制御モデル (infinite-horizon optimal feedback control model)[外部リンク] である．

    クロネッカー積

    定義と例

    以下，Wikipediaの説明から引用する．

    \(A = (a_{ij})\) を \(m \times n\) 行列，\(B = (b_{kl})\) を \(p \times q\) 行列とすると、それらのクロネッカー積 \(A \otimes B\) は

    $$ A\otimes B={\begin{bmatrix}a_{11}B&\cdots &a_{1n}B\\\vdots &\ddots &\vdots \\a_{m1}B&\cdots &a_{mn}B\end{bmatrix}} $$

    で与えられる \(mp \times nq\) 区分行列である．例えば，

    $$ \begin{bmatrix}1&2\\3&4\end{bmatrix} \otimes {\begin{bmatrix}0&5\\6&7\end{bmatrix}}={\begin{bmatrix}1\cdot 0&1\cdot 5&2\cdot 0&2\cdot 5\\1\cdot 6&1\cdot 7&2\cdot 6&2\cdot 7\\3\cdot 0&3\cdot 5&4\cdot 0&4\cdot 5\\3\cdot 6&3\cdot 7&4\cdot 6&4\cdot 7\\\end{bmatrix}}={\begin{bmatrix}0&5&0&10\\6&7&12&14\\0&15&0&20\\18&21&24&28 \end{bmatrix}} $$

    のような計算が成り立つ．

    実装

    Pythonではnumpy.kron()，JuliaではKronecker.jlを用いる．以下ではPythonを使用する．

    In [1]:
    import numpy as np
np.set_printoptions(precision=3, suppress=True) # print時の設定
    In [2]:
    A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 5], [6, 7]])
print(np.kron(A, B))

    [[ 0  5  0 10]
 [ 6  7 12 14]
 [ 0 15  0 20]
 [18 21 24 28]]

    行列方程式の解法

    1. AXB = C

    \(A, B, C\)が与えられていて，\(X\)を未知とするときの方程式 \(AXB = C\)を考える．この場合，\(A, B\)は正則であると仮定すると，\(X=A^{-1}CB^{-1}\)と求まる．ただし，以下では後の問題のためにRoth's column lemma (vec-trick)を用いた方法について考えよう．Roth's column lemmaは，先程と同様の方程式 \(AXB = C\)を考えると，この方程式は

    $$ (B^\top \otimes A)\text{vec}(X) = \text{vec}(AXB)=\text{vec}(C) $$

    の形に書き下すことができる，というものである．\(\text{vec}(\cdot)\)はvec作用素（行列を列ベクトル化する作用素）である．numpyだとflatten()で実現できる．Roth's column lemmaを用いれば，\(AXB = C\)の解は

    $$ X = \text{vec}^{-1}\left((B^\top \otimes A)^{-1}\text{vec}(C)\right) $$

    として得られる．ただし，\(\text{vec}(\cdot)^{-1}\)は列ベクトルを行列に戻す作用素(inverse of the vectorization operator)である．numpyだとreshape()で実現できる (numpyのベクトルはデフォルトで行ベクトルであることに注意)．2つの作用素をまとめると，

    $$ \begin{align} \text{vec} &: R^{m\times n}\to R^{mn}\\ \text{vec}^{−1} &: R^{mn}\to R^{m×n} \end{align} $$

    であり，\(\text{vec}^{−1}\left(\text{vec}(X)\right)=X\ (\text{for all}\ X\in R^{m\times n})，\text{vec}\left(\text{vec}^{−1}(x)\right)=x\ (\text{for all}\ x \in R^{mn})\)となる．

    以下では\(P, Q, R\)が与えられている際に，方程式\(PYQ=R\)の解\(Y\)を求める．

    In [3]:
    m = 4
P = np.random.randn(m, m)
Q = np.random.randn(m, m)
R = np.arange(m**2).reshape(m, m)
print(R)

    [[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

    上記に基づいて方程式を解く．

    In [4]:
    Y = (np.linalg.inv(np.kron(P, Q.T)) @ R.flatten()).reshape(m, m)
print(P @ Y @ Q)

    [[ 0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]
 [12. 13. 14. 15.]]

    少々元の式と実装が異なるが，計算はできているので良しとしておく．また，\(Y=P^{-1}RQ^{-1}\)としても解けることも示す．

    In [5]:
    Y2 = np.linalg.inv(P) @ R @ np.linalg.inv(Q)
print(P @ Y2 @ Q)

    [[-0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]
 [12. 13. 14. 15.]]

    2. AXB + CXD = E

    この記事の本題．\(A, B, C, D, E\)が与えられていて，\(X\)を未知とするときの方程式 \(AXB + CXD = E\)を考えると，この方程式は

    $$ \begin{align} & (B^\top \otimes A)\text{vec}(X) + (D^\top \otimes C)\text{vec}(X) =\text{vec}(E)\\ & X = \text{vec}^{-1}\left((B^\top \otimes A + D^\top \otimes C)^{-1}\text{vec}(E)\right) \end{align} $$

    として解ける．例えば連続的リアプノフ方程式 (continuous Lyapunov equation) \(AX+XA^{H}+Q=0\)は\(AXB + CXD = E\)の特別な場合であるので，

    $$ \begin{align} & (I^\top \otimes A)\text{vec}(X) + (\bar{A} \otimes I)\text{vec}(X) =-\text{vec}(Q)\\ & X = -\text{vec}^{-1}\left((I^\top \otimes A + \bar{A} \otimes I)^{-1}\text{vec}(Q)\right) \end{align} $$

    と解ける．ただし，\(A^H\)は\(A\)の随伴行列，\(\bar{A}\)は\(A\)の各成分を複素共役で置き換えた行列，\(I\)は単位行列である．他にはシルベスター方程式 (Sylvester equation) \(AX+XB=C\)も同様に解ける．

    以下では\(PZQ + SZT = R\)の解\(Z\)を求める．

    In [6]:
    S = np.random.randn(m, m)
T = np.random.randn(m, m)
    In [7]:
    Z = (np.linalg.inv(np.kron(P, Q.T) + np.kron(S, T.T)) @ R.flatten()).reshape(m, m)
print(P @ Z @ Q + S @ Z @ T)

    [[ 0.  1.  2.  3.]
 [ 4.  5.  6.  7.]
 [ 8.  9. 10. 11.]
 [12. 13. 14. 15.]]

    3. XAX = B

    クロネッカー積は使わないが，\(XAX=B\)の解法を記載しておく．\(A\)が正則であると仮定すると，行列平方根（の一つ）\(A^{1/2}\)が求まり，両辺の両側から\(A^{1/2}\)をかけることで，\(\left(A^{1/2}X A^{1/2}\right)^2 = A^{1/2}B A^{1/2}\)が得られる．よって，

    $$ X = A^{-1/2}\left(A^{1/2}B A^{1/2}\right)^{1/2}A^{-1/2} $$

    と求まる．\(A^{1/2}\)の計算にはscipy.linalg.sqrtm()を用いる．なお，一般に\(n\)次正則行列の場合に行列平方根は（よって\(XAX = B\)の解も）\(2^n\)個存在する．以下では\(XPX=R\)を解く．

    In [8]:
    from scipy.linalg import sqrtm
    In [9]:
    Psqrt = sqrtm(P)
Pinvsqrt = np.linalg.inv(Psqrt)
X = Pinvsqrt @ sqrtm(Psqrt @ R @ Psqrt) @ Pinvsqrt
print(X @ P @ X)

    [[-0.-0.j  1.-0.j  2.+0.j  3.-0.j]
 [ 4.-0.j  5.-0.j  6.+0.j  7.-0.j]
 [ 8.+0.j  9.+0.j 10.-0.j 11.-0.j]
 [12.-0.j 13.+0.j 14.+0.j 15.+0.j]]

    参考文献

    Tags Python