大阪大学医学部 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