Help us understand the problem. What is going on with this article?

級数 S = X+AXA*+A^2XA^2*+... を求める方法?(離散リアプノフ方程式の数値解導出)

Qiitaだとちょっとおかしい表記になるので気になる人はGithub上で見て。

追記:「離散のリヤプノフ方程式」です。完全に恥さらしですが反省のために残しておきます。

目的

以下の級数(名前知らない)の無限和を求めたい。

$$
\begin{align}
S &= X + AXA^\top + A^2 X (A^2)^\top + A^3 X (A^3)^\top +...\\
&= \sum^{\infty}_{n=0} A^n X (A^n)^\top
\end{align}
$$

条件として,

  • $A$はフルランクである。
  • $A$の全ての固有値の絶対値は1より小さい:$0<|eig(A)|<1$
  • $A,X \in \mathcal{R}^{n\times n}$
  • $X$は対称行列

とする。
なんか足りなそうな条件あったら教えてください。

数式の変換

ここで$S$の両辺から$A$と$A^\top$をかけた式とを見比べる。

$$
\begin{align}
S &= X + AXA^\top + A^2 X (A^2)^\top + ...\\
ASA^\top &= AXA^\top + A^2 X (A^2)^\top + A^3 X (A^3)^\top + ...\\
\end{align}
$$

この時。$\lim_{n\rightarrow \infty}A^nX(A^n)^\top=0$になるため,引き算をして

$$
\begin{align}
S = ASA^\top + X
\end{align}
$$

を得る。これを方程式(*)と呼ぼう。(Markdownの表記都合上)

これを$S$について解ければ良さそう。

解法1:Sylvestar 方程式 への変換

Sylvestar方程式とは以下の方程式のことをいう。

$$
AX + XB = C
$$

なお,この方程式が一意的な解 $X$ を持つのは、行列 $A$ と $−B$ が共通する固有値を持たないときのみに限られる。

Jameson, 1968の論文に解法がある(非常にめんどい!)。

余談だが$A=B$の時,この式はご存知Lyapnov(リアプノフ)方程式となる。
こういった名前を知らないと調べる時に非常に不便するので,勉強とはこのためにあると言っても過言ではない。

先程の方程式(*)で$A$がフルランクなので,求めたい$S$について,

$$
A^{-1}S + S(-A^\top) = A^{-1}X
$$

とできる。
解がある条件は,$A^{-1}$と$A^\top$が固有値を共有しない時,である。

  • $A$がフルランクより逆行列の固有値は元の行列の逆数になる。
  • 転置行列の固有値は元の行列と変わらない

ということを考えると,$A$が固有値に1を持たない or 逆数の関係にある固有値を持たないという条件になり,これは前提にて満たされている。

ということでscipyやmatlabのソルバを使えば一応この値が計算できるということになる。

解法2:固有ベクトルの分解を用いて解く(名前知らない)

名前は知らないが$A$の固有ベクトルを用いて方程式(*)を解くことができる。

今,線形写像$L(\Phi)$を次のように定める。

$$
L(\Phi)=\Phi-A\Phi A^\top
$$

この写像に$\Phi$に$S$を代入した時,方程式(*)より

$$
L(S)=S-ASA^\top=X
$$

となる。

また,$A$の固有値と固有ベクトルをそれぞれ$\lambda_i, u_i (i=1,...,n)$と置くと,

$$
\begin{align}
L(u_iu_j^\top)&=u_iu_j^\top-Au_iu_j^\top A^\top\\
&=u_iu_j^\top-\lambda_iu_iu_j^\top\lambda_j\\
&=(1-\lambda_i\bar{\lambda_j}) u_iu_j^\top
\end{align}
$$

という風に写像される。$u_iu_j^\top$はより一般化して$u_iu_j^*$(エルミート)とした方が良いのかな…?その場合バーをつけている固有値も複素共役である必要がある。

次もちょっと自信がないのだが,

  • $A$がフルランクなので固有ベクトルの積で$\mathcal{R}^{n\times n}$の空間を貼れる?

ため,それぞれ

$$
X=\sum_{ij} \tilde{x}{ij} u_iu_j^\top\\
S=\sum
{ij} \tilde{s}_{ij} u_iu_j^\top
$$

と分解できる。$\tilde{x}{ij},\tilde{s}{ij}$は各要素のゲインを表すスカラである。

先程の線形写像の関係式を用いれば,求めたいSの要素は

$$
\tilde{s}{ij} = \frac{1}{(1-\lambda_i\bar{\lambdaj})}\tilde{x}_{ij}
$$

となる。カンのいい人は$A$が固有値に1を持たない or 逆数の関係にある固有値を持たないという条件がこの分母を0にしない条件と等価であることを察せるだろうか。

したがって,$\tilde{x}_{ij}$を求めればいいが,固有ベクトルの直交条件より

$A=U\Sigma U^\top$の$U=[u_1 u_2 ... u_n]$を用いて

$$
\begin{align}
X=U \tilde{X}U^\top = U\begin{bmatrix}
\tilde{x}{11}&\tilde{x}{12}& \cdots \
\tilde{x}{21}&\tilde{x}{22}& \cdots \
\vdots & \vdots & \cdots
\end{bmatrix}
U^\top
\end{align}
$$

となるため,
同様の考察をすれば,$\Lambda=[\lambda_1 \lambda_2 ... \lambda_n]^\top[\lambda_1 \lambda_2 ... \lambda_n]$を用いて$S=U\tilde{S}U^\top$の$\tilde{S}$は

$$
\tilde{S}=G.*\tilde{X}^\top\
\mbox{while, } \ G=1./(ones-\Lambda)
$$

と計算できる。ここで1で引かれているので任意の

onesは全ての要素が1の行列,.*と./は要素ごとの積と除算を表す。(matlabコマンドを使う人は想像付きやすいかと)

実際にやってみる

実際にやってみよう。

A = [1.5 1 ; -0.7 0 ]
X = [1 0.5 ; 0.5 0.25 ]

の2次元行列として,

  1. 級数和を地道に計算(数値解?)
  2. シルベスタ方程式を用いて求める(解1)
  3. 固有値ベクトルを用いた方法で求める(解2)

として結果を比較してみよう

% 変数設定
A = [1.5 1 ; -0.7 0 ]
X = [1 0.5 ; 0.5 0.25 ]

eig(A)
A =

    1.5000    1.0000
   -0.7000         0


X =

    1.0000    0.5000
    0.5000    0.2500


ans =

   0.7500 + 0.3708i
   0.7500 - 0.3708i

1.愚直に計算する

愚直に計算するとだいたい,

S1 =

18.8694  -11.3596
  -11.3596    9.4952

くらいの値になる。

S1 = X;

for i = 1:21
S1 = S1 + A^i*X*A^i';
    if mod(i,5)==1
        disp(['sum of ',num2str(i),'th term'])
        display(S1)
    end
end
sum of 1th term

S1 =

    5.0000   -0.9000
   -0.9000    0.7400

sum of 6th term

S1 =

   17.2341  -10.3678
  -10.3678    8.6879

sum of 11th term

S1 =

   18.6998  -11.2240
  -11.2240    9.3343

sum of 16th term

S1 =

   18.8266  -11.3207
  -11.3207    9.4558

sum of 21th term

S1 =

   18.8694  -11.3596
  -11.3596    9.4952

2.シルベスタ方程式で解く

これはコマンドを叩くだけ。

値は

S2 =

   18.8802  -11.3672
  -11.3672    9.5013

で同じ感じ。

S2 = sylvester(inv(A),-A',A\X)
S2 =

   18.8802  -11.3672
  -11.3672    9.5013

3.固有ベクトルを用いて解く

めんどくさいようでコマンド自体はそこまで複雑でもない。

S3 =

  18.8802 - 0.0000i -11.3672 + 0.0000i
 -11.3672 + 0.0000i   9.5013 - 0.0000i

複素数が混じっているものの同じ様な値を取れた。

[U,D] = eig(A);

Sigma = diag(D);

G = 1./(ones(2,2)-Sigma*Sigma');

X_tilde = inv(U)*X*inv(U');
S_tilde = G.*X_tilde;

S3 = U*S_tilde*U'


S3 =

  18.8802 - 0.0000i -11.3672 + 0.0000i
 -11.3672 + 0.0000i   9.5013 - 0.0000i

まとめ

一応この方程式を解くことはできたが,この級数の導出についてまともな説明のある論文や教科書があったら教えてほしい。

絶対あると思うんだよなぁ。

→ 離散リアプノフ方程式でした

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away