行列方程式 - シルベスター方程式とリャプノフ方程式

行列を未知数とする方程式を行列方程式と呼びます。$AX=B$ のような単純な形であれば $X=A^{-1}B$ で解けますが、未知行列が複数の行列に挟まれた形になると話が複雑になります。ここではシルベスター方程式とリャプノフ方程式という代表的な 2 つの行列方程式を扱います。

シルベスター方程式

シルベスター方程式は

$$AX+XB=C$$

の形をした行列方程式です。$A$ は $m\times m$、$B$ は $n\times n$、$C$ は $m\times n$ の既知の行列であり、$X$ は $m\times n$ の未知行列です。未知行列 $X$ が $A$ と $B$ の両方に挟まれている点が特徴的で、$X$ について単純に因数分解できません。

この方程式は 1884 年にジェームズ・ジョゼフ・シルベスターによって研究されました。制御理論や安定性解析など、工学的な応用が広いことで知られています。

解の存在条件

シルベスター方程式が任意の $C$ に対して一意な解を持つための必要十分条件は、$A$ と $-B$ が共通の固有値を持たないことです。すなわち

$$\sigma (A)\cap \sigma (-B)=\varnothing$$

が成り立つことです。ここで $\sigma (\cdot )$ は固有値の集合（スペクトル）を表します。同じことですが、$A$ の固有値 ${\lambda }_i$ と $B$ の固有値 ${\mu }_j$ に対して

$${\lambda }_i+{\mu }_j\ne 0(\text{すべての }i,j)$$

が成り立てば一意な解が存在します。

なぜこの条件になるのかを理解するために、シルベスター方程式をベクトル化して考えます。$\mathrm{vec}(X)$ を行列 $X$ の列を縦に積んだベクトルとすると、クロネッカー積の性質

$$\mathrm{vec}(AXB)=(B^T\otimes A)\mathrm{vec}(X)$$

を使って $AX+XB=C$ は

$$(I_n\otimes A+B^T\otimes I_m)\mathrm{vec}(X)=\mathrm{vec}(C)$$

と書き直せます。これは $mn\times mn$ の通常の連立一次方程式です。この係数行列 $I_n\otimes A+B^T\otimes I_m$ が正則であるための条件が、$A$ と $-B$ の固有値が重ならないことに対応しています。$I_n\otimes A$ の固有値は $A$ の固有値を $n$ 回繰り返したもの、$B^T\otimes I_m$ の固有値は $B$ の固有値を $m$ 回繰り返したものであり、和の固有値は ${\lambda }_i+{\mu }_j$ の形になるためです。

小さな例

$2\times 2$ の具体例で計算してみます。

$$A=\left(\begin{array}{cc}1& 0\\ 0& 2\end{array}\right),B=\left(\begin{array}{cc}3& 0\\ 0& 4\end{array}\right),C=\left(\begin{array}{cc}4& 5\\ 6& 12\end{array}\right)$$

$A$ の固有値は $\{1,2\}$、$-B$ の固有値は $\{-3,-4\}$ であり、共通部分はないため一意な解が存在します。

$X=6$ とおいて $AX+XB=C$ を成分ごとに書き下すと

$$\left(\begin{array}{cc}{x}_{11}& x_{12}\\ 2x_{21}& 2x_{22}\end{array}\right)+\left(\begin{array}{cc}3x_{11}& 4x_{12}\\ 3x_{21}& 4x_{22}\end{array}\right)=\left(\begin{array}{cc}4& 5\\ 6& 12\end{array}\right)$$

各成分を比較すると

$$4x_{11}=4,5x_{12}=5,5x_{21}=6,6x_{22}=12$$

となり

$$X=\left(\begin{array}{cc}1& 1\\ 6/5& 2\end{array}\right)$$

が得られます。$A$ と $B$ が対角行列の場合は各成分が独立に解けるため、$({\lambda }_i+{\mu }_j)x_{ij}=c_{ij}$ という形になります。${\lambda }_i+{\mu }_j\ne 0$ であれば各成分を割り算で求められるわけです。

リャプノフ方程式

リャプノフ方程式はシルベスター方程式の特殊な場合であり、次の形をしています。

$$AX+X{A}^T=C$$

シルベスター方程式で $B=A^T$ としたものです。$C$ が対称行列であれば解 $X$ も対称行列になります。

リャプノフ方程式の解の存在条件はシルベスター方程式の条件から導かれます。$B=A^T$ なので $-B=-A^T$ であり、$A^T$ の固有値は $A$ の固有値と同じです。したがって条件は

$${\lambda }_i+\overline{{\lambda }_j}\ne 0(\text{すべての }i,j)$$

です（実行列の場合は ${\lambda }_i+{\lambda }_j\ne 0$）。特に、$A$ のすべての固有値の実部が負（安定行列）であれば、この条件は自動的に満たされます。

連続リャプノフ方程式と安定性

制御理論では、線形システム $\dot{\mathbf{x}}=A\mathbf{x}$ の安定性を判定するためにリャプノフ方程式が使われます。具体的には次の定理が成り立ちます。

$A$ のすべての固有値が負の実部を持つことと、任意の正定値対称行列 $Q$ に対して

$$A^{T}P+PA=-Q$$

を満たす正定値対称行列 $P$ が存在することは同値です。

この定理が有用な理由は、大きな行列の固有値を計算するのは数値的に不安定になりうるのに対して、リャプノフ方程式を解くことは数値的に安定なアルゴリズムが存在するからです。

数値解法

シルベスター方程式を数値的に解く標準的な方法はバートレス・スチュアート（Bartels–Stewart）アルゴリズムです。手順は次のとおりです。

$A=UR{U}^T$, $B=VS{V}^T$（$R$, $S$ は上三角行列、$U$, $V$ は直交行列）とシュール分解すると、元の方程式は

$$RY+YS=\tilde{C}$$

の形に帰着されます（$Y=U^{T}XV$, $\tilde{C}=U^{T}CV$）。$R$ と $S$ が上三角行列であるため、この方程式は後退代入で列ごとに解けます。計算量は $O(m^3+n^3)$ です。

Python では SciPy の scipy.linalg.solve_sylvester で直接解くことができます。

import numpy as np
from scipy.linalg import solve_sylvester

A = np.array([[1, 0], [0, 2]])
B = np.array([[3, 0], [0, 4]])
C = np.array([[4, 5], [6, 12]])

X = solve_sylvester(A, B, C)
print(X)

リャプノフ方程式には scipy.linalg.solve_continuous_lyapunov が用意されています。

from scipy.linalg import solve_continuous_lyapunov

A = np.array([[-1, 0], [0, -2]])
Q = np.eye(2)

P = solve_continuous_lyapunov(A, -Q)
print(P)

離散リャプノフ方程式

連続時間のリャプノフ方程式とは別に、離散時間システム ${\mathbf{x}}_{k+1}=A{\mathbf{x}}_k$ に対応する離散リャプノフ方程式

$$AX{A}^T-X+Q=O$$

も存在します。この方程式の解が正定値になるための条件は、$A$ のすべての固有値の絶対値が 1 未満（$|{\lambda }_i|<1$）であることです。連続時間では「実部が負」、離散時間では「絶対値が 1 未満」という安定性条件の違いに対応しています。