2018-08-24

ベクトル場

この記事では、多様体上のベクトル場について解説します。なお、この記事を理解するためには多様体や接空間などの概念を理解しておくことが必要です。それらについては

を参考にしてください。

多様体上のベクトル場
ベクトル場全体の集合の代数構造
積分曲線
二つの多様体上のベクトル場の関係性
参考文献

多様体上のベクトル場

多様体 $M$ 上のベクトル場 $X$ とは、 $X: M\ni p\mapsto X(p) \in T_p M$ という対応のことです。 $X(p)$ をよく $X_p$ と書きます。この記事では基本的には $X(p)$ と書くことにして、 $X(p)$ だと読みづらくなると思われるところで $X_p$ を使うことにします。また、 $m$ 次元多様体 $M$ の座標近傍系が $\{ U_{\alpha}; x_1^{\alpha}, x_2^{\alpha}, \ldots, x_m^{\alpha} \}$ であるとき、 $M$ 上のベクトル場 $X$ が滑らか( $C^{\infty}$ 級 )であるとは $X$ を $U_{\alpha}$ 上で $X= \sum_{i=1}^m \xi_i^{\alpha} \frac{\partial}{\partial x_i^{\alpha}}$ と局所座標表示したときに、成分 $\xi_1^{\alpha}, \xi_2^{\alpha}, \ldots, \xi_m^{\alpha}$ が $U_{\alpha}$ 上で滑らかであるときに言います。以下では単にベクトル場と言ったら、滑らかなベクトル場のことを言ってることにします。

ベクトル場を可視化した例を2次元ユークリッド空間 ${\bf R}^2$ という馴染みのある多様体で与えておきます。多様体 ${\bf R}^2$ の(局所)座標は $(x,y)$ だとします。このとき、多様体 ${\bf R}^2$ 上のベクトル場 $X$ を $X(p) = (x^2+y)\left( \frac{\partial}{\partial x} \right)_p + xy \left( \frac{\partial}{\partial y}\right)_p$ としますと、以下のような感じでこのベクトル場 $X$ を可視化できます。

f:id:ogyahogya:20180824090812p:plain

つまり、 ${\bf R}^2$ 上のベクトル場が与えられると各点に馴染みのある矢印のベクトルがくっついてるというように考えることができます。したがって、天気予報でよく見かける以下のような図は地球の表面の一部を ${\bf R}^2$ で表現して風というベクトル場を可視化していると考えることができます。

f:id:ogyahogya:20180824090255p:plain

ベクトル場全体の集合の代数構造

多様体 $M$ 上のベクトル場全体の集合を $\mathfrak{X}(M)$ と書くことにして、任意の $X, Y\in \mathfrak{X}(M)$ の和 $X+Y$ を任意の $p\in M$ に対して\begin{align} (X+Y)(p):= X(p)+ Y(p) \end{align}と定義します。また、任意の $f\in C^{\infty}(M)$ と $X\in \mathfrak{X}(M)$ に対し、 $X$ の $f$ 倍を\begin{align} (fX)(p):=f(p)X(p) \end{align} と定義します。このような和と積を $\mathfrak{X}(M)$ に導入すると、多様体 $M$ 上の滑らかな関数の集合 $C^{\infty}(M)$ は代数学の言葉では環の性質を満たすので、 $\mathfrak{X}(M)$ は $C^{\infty}(M)$ 加群ということになります。また関数 $f$ が定数関数、つまり任意の $p\in M$ に対して $f(p)$ が一定なものも $C^{\infty}$ 級の関数だと考えられますので、 $\mathfrak{X}(M)$ は上の和と積で ${\bf R}$ 上のベクトル空間の構造も入っていることになります。

さらに、任意の $f\in C^{\infty}(M)$ と $X\in \mathfrak{X}(M)$ に対し、\begin{align} (Xf)(p):= X_p f \end{align} と定義します。ここで $fX$ は $\mathfrak{X}(M)$ の要素だったのに対して $Xf$ は $C^{\infty}(M)$ の要素であることに注意してください。

ベクトル場 $X$ は局所座標表示すると、 $X= \sum_{i=1}^m \xi_i^{\alpha} \frac{\partial}{\partial x_i^{\alpha}}$ というように書けたので1階の微分作用素だとみなせます。次に、 $X, Y\in \mathfrak{X}(M)$ , $f\in C^{\infty}(M)$ として $(XY)(f):=X(Y(f))$ を定義して $XY$ を局所座標表示すると2階の微分作用素だとみなせます。同様に $YX$ も2階の微分作用素だとみなせます。ところが、 $XY$ と $YX$ の差 $XY-YX$ は1階の微分作用素となります。つまり、 $XY-YX\in \mathfrak{X}(M)$ となります。

この $XY-YX$ は重要で、 $[X,Y]:=XY-YX\in \mathfrak{X}(M)$ は $X$ と $Y$ のLie(リー)括弧積と呼ばれています。リー括弧積は以下の性質があります。

f:id:ogyahogya:20180824103811p:plain

ただし、上の $X, Y, Z$ は $\mathfrak{X}(M)$ の任意の要素です。一般に、 ${\bf R}$ 上のベクトル空間 $V$ と $[\cdot, \cdot]:V\times V\rightarrow V$ が与えられときに、リー括弧積の性質(i), (ii), (iii) を満たすならベクトル空間 $V$ はリー代数と呼ばれています。したがって、多様体 $M$ 上のベクトル場全体の集合 $\mathfrak{X}(M)$ はリー代数という代数構造も持ちます。

積分曲線

上で紹介した天気予報でよく出てくる風のベクトル場を可視化した矢印たちを眺めていると曲線が見えてくるかもしれません。以下では、そのような曲線を数学的に考えます。

で説明しているように、多様体 $M$ 上の曲線 $c:(a,b)\rightarrow M$ が与えられると曲線 $c$ の速度ベクトルを考えることができて、速度ベクトルは曲線上の点の接空間の元、つまり接ベクトルになっていました。曲線上の各点の速度ベクトルを集めると曲線上のベクトル場が定義できたことになります。この曲線 $c$ の速度ベクトルを集めてできるベクトル場を $\dot{c}$ と書くことにします。このとき、曲線 $c$ の積分曲線は以下のように定義されます。

f:id:ogyahogya:20180824112030p:plain

上の $\dot{c}(t)=X(c(t))$ が具体的には何ものなのかを見るために、曲線 $c$ の一部が多様体 $M$ の座標近傍 $(U; x_1,x_2,\ldots, x_m)$ 上にあるとして、 $\dot{c}(t)=X(c(t))$ を局所座標系を使って表してみます。曲線 $c$ を $U$ 上で $c(t) = (x_1(t),x_2(t),\ldots, x_m(t))$ と書くと、\begin{align} \dot{c}(t) = \sum_{i=1}^m \frac{dx_i}{dt}(t) \left( \frac{\partial}{\partial x_i} \right)_{c(t)} \end{align} となります。また、 $X= \sum_{i=1}^m f_i\frac{\partial}{\partial x_i}$ を $X$ の居所座標表示とすると、\begin{align} X(c(t))=\sum_{i=1}^m f_i(c(t))\left(\frac{\partial}{\partial x_i}\right)_{c(t)} \end{align} と書けます。よって、 $\dot{c}(t) = X(c(t))$ は局所的には\begin{align} \frac{dx_i}{dt}(t) = f_i(x_1(t),x_2(t),\ldots,x_m(t))\quad (i=1,2,\ldots, m) \end{align} という微分方程式のことであることが分かりました。

微分方程式の解の存在定理より、初期条件 $(x_1(0),x_2(0),\ldots, x_m(0))$ を与えると、\begin{align} \frac{dx_i}{dt}(t) = f_i(x_1(t),x_2(t),\ldots,x_m(t))\quad (i=1,2,\ldots, m) \end{align} を満たす解 $(x_1(t),x_2(t),\ldots, x_m(t))$ が $-\epsilon＜t＜\epsilon$ で存在することが分かります。ただし、 $\epsilon$ は十分小さな正の実数です。また微分方程式の解の一意性定理より、 $a＜t＜b$ で $(x_1(t),x_2(t),\ldots, x_m(t))$ が解で $c＜t＜d$ で $(y_1(t), y_2(t),\ldots, y_m(t))$ が解で $(x_1(0),x_2(0),\ldots, x_m(0)) = (y_1(0),y_2(0),\ldots, y_m(0))$ なら $(a,b)\cap (c,d)$ で $(x_1(t),x_2(t),\ldots, x_m(t)) = (y_1(t), y_2(t),\ldots, y_m(t))$ となります。以上は局所座標を使った議論ですが、局所座標を使わないで同じことを述べると以下のようになります。

f:id:ogyahogya:20180824125617p:plain

二つの多様体上のベクトル場の関係性

集合 $M$ と $N$ を多様体だとして、滑らかな写像 $F: M\rightarrow N$ とベクトル場 $X\in \mathfrak{X}(M)$ が与えられているとします。このとき、 $F$ の点 $p\in M$ における微分写像 $(dF)_p: T_p M \rightarrow T_{F(p)} N$ によって、ベクトル $(dF)_p(X_p)\in T_{F(p)} N$ を得ることができます。しかし、そのような $(dF)_p(X_p)$ が $N$ 上のベクトル場を定義するとは限りません。実際に、写像 $F$ が全射でないとすると、点 $q\in N\backslash F(M)$ には $F$ の微分によってベクトルを定めることはできません。また、写像 $F$ が単射でないとすると、 $F$ は $M$ の異なる2点を $N$ の同一の点に移すことがあり、そのような $M$ の2点での接ベクトルは $F$ の微分写像によって異なる接ベクトルに移されることがあるからです（接ベクトルの始点は同じなのに）。したがって、滑らかな写像 $F: M\rightarrow N$ とベクトル場 $X\in \mathfrak{X}(M)$ 、 $Y\in \mathfrak{X}(N)$ が与えられたときに、 $X$ と $Y$ が $F$ で関係付けられるとは一般には言えません。ベクトル場 $X\in \mathfrak{X}(M)$ と $Y\in \mathfrak{X}(N)$ が任意の $p\in M$ に対して以下のように滑らかな写像 $F$ で関係付けられるとき、ベクトル場 $X$ と $Y$ は $F$ -関係を持つと言います。

f:id:ogyahogya:20190412141134p:plain

現段階では、 $X\in \mathfrak{X}(M)$ に $F$ -関係な $N$ 上のベクトル場が存在することは言えていないのですが、実は $F:M\rightarrow N$ が微分同相写像なら任意の $X\in \mathfrak{X}(M)$ に $F$ -関係な $N$ 上のベクトル場が存在することが言えます。実際に、任意の $q\in N$ に対して\begin{align} Y_q = (dF)_{F^{-1}(q)} (X_{F^{-1}(q)}) \end{align} とベクトル場 $Y$ を定めると、 $X$ と $Y$ は $F$ -関係を持ち、 $Y: N\rightarrow TN$ が滑らかであることは、 $N\rightarrow M\rightarrow TM\rightarrow TN$ という滑らかな写像の合成だからです。また、 $Y$ の定め方より $X$ に $F$ -関係な $N$ 上の滑らかなベクトル場は一意だということも分かります。このような $X$ に $F$ -関係な $N$ 上の滑らかなベクトル場を $F_* X$ と書いて $F$ による $X$ の押し出しと言います。上の議論から、任意の $q\in N$ に対して\begin{align} (F_*X)_q = (dF)_{F^{-1}(q)} (X_{F^{-1}(q)}) \end{align}です。

参考文献

(1) ベクトル場の説明を参考にしました。

多様体の基礎 (基礎数学5)

作者: 松本幸夫
出版社/メーカー: 東京大学出版会
発売日: 1988/09/22
メディア: 単行本
購入: 7人クリック: 36回
この商品を含むブログ (33件) を見る

(2) 微分方程式の解の存在と一意性などについて詳しく書いています。

常微分方程式の安定性 (実教理工学全書)

作者: 山本稔
出版社/メーカー: 実教出版
発売日: 1979/01
メディア: 単行本
この商品を含むブログを見る

(3) ベクトル場の関係性についての議論を参考にしました。

Introduction to Smooth Manifolds (Graduate Texts in Mathematics)

作者: John Lee
出版社/メーカー: Springer
発売日: 2012/08/24
メディア: ハードカバー
この商品を含むブログを見る

2018-05-14

可制御可観測な線形システム全体の集合は多様体になる

制御理論多様体とか集合と位相

この記事では、可制御・可観測な線形システム全体の集合は多様体になることを説明します。ここで言う線形システムとは、

で解説した

f:id:ogyahogya:20180511094841p:plain

のことですが、状態方程式表現は行列の三つ組 $(A,B,C)\in {\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ で決定するので、一つの線形システムは ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ の一点である\begin{align} (A,B,C)\end{align}

のことだ考えられます。そうすると、線形システム全体の集合とは\begin{align} {\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}\end{align}

のこととなります。この線形システム全体の集合の中には、

で解説した可制御な線形システムや可観測な線形システムが含まれることになり、上の解説記事を注意深く読むと可制御な線形システム全体の集合とは、 \begin{align} S_c:=\{ (A,B,C)\, |\, {\rm rank} \begin{pmatrix} B & AB & \cdots & A^{n-1}B \end{pmatrix}=n\} \end{align}

のことであり、可観測な線形システム全体の集合とは、 \begin{align} S_o:=\{ (A,B,C)\, |\, {\rm rank} \begin{pmatrix} C \\ CA \\ \vdots \\ CA^{n-1} \end{pmatrix}=n\} \end{align}

のことだということになります。つまり、この記事では、

\begin{align} S_c\cap S_o\end{align}

という可制御可観測な線形システム全体の集合が多様体になるということを説明します。多様体については、

を参照してください。

多様体の開集合は多様体になる。

結論から先に言うと、可制御可観測な線形システム全体の集合 $S_c\cap S_o$ が多様体であることは、 $S_c\cap S_o$ が多様体である ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ の開集合であることから分かります。そこで、多様体の開集合は多様体になるという事実を解説しておきます。あとで、 $S_c\cap S_o$ が ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ の開集合だということを説明します。

そもそも開集合とは何かというと、

で説明したように、位相空間の要素のことでした。位相空間 $(X,\mathcal{O})$ が与えられると、 $X$ の任意の部分集合 $Y$ に対して、

\begin{align} \mathcal{O}_Y:= \{ U\cap Y\,|\, U\in \mathcal{O} \} \end{align}

とすることで、 $\mathcal{O}_Y$ が $Y$ の位相になります。このような $\mathcal{O}_Y$ を $X$ の位相 $\mathcal{O}$ から導かれた $Y$ の相対位相と呼んで、位相空間 $(Y,\mathcal{O}_Y)$ を $(X,\mathcal{O})$ の部分空間と呼びます。多様体は上のリーマン多様体の記事で定義したように任意の異なる2点を開集合で分離できるという性質を持ったハウスドルフ空間でした。実は、ハウスドルフ空間 $X$ の部分空間 $Y$ もハウスドルフ空間になります。なぜなら、任意の $a, b\in Y$ , $a\neq b$ に対して、 $X$ はハウスドルフ空間なので、 $a\in U, b\in V, U\cap V=\emptyset$ となる $X$ の開集合 $U, V$ が存在して、 $U':= U\cap Y, V':=V\cap Y$ は $Y$ の開集合で $a\in U', b\in V', U'\cap V'=\emptyset$ となるからです。

上のことから、多様体 $(X,\mathcal{O})$ の開集合 $Y$ に $\mathcal{O}$ から導かれた相対位相を導入すると開集合 $Y$ はハウスドルフ空間となります。さらに、多様体 $X$ の $C^{\infty}$ 級の座標近傍系が $\{ (U_{\lambda}, \phi_{\lambda})\}_{\lambda\in \Lambda}$ だとすると、これを $Y$ に制限した $\{ (U_{\lambda}\cap Y, \phi_{\lambda}|U_{\lambda}\cap Y)\}_{\lambda\in \Lambda}$ は $Y$ の $C^{\infty}$ 級の座標近傍系になります。したがって、多様体 $X$ の任意の開集合は $X$ の位相から導かれた相対位相を導入することで多様体になります。

距離空間は位相空間になる

応用上もっとも重要な空間は距離が定められた空間であり、その距離を使って位相を導入できる、つまり、空間の異なる点たちの近さを議論できるようになります。あとで線形システム全体の集合にも距離を導入して位相空間としますので、ここでは簡単に距離空間について説明しときます。

距離空間の定義は以下の通りです。

f:id:ogyahogya:20180513132711p:plain

ベクトル空間にノルムが導入されたノルム空間は距離空間になります。ここで、ノルム空間の定義は以下の通りです。

f:id:ogyahogya:20180513132822p:plain

実際に、ベクトル空間 $X$ にノルム $||\cdot ||$ が導入されているとすると、 $X$ 上の任意の2点 $x,y$ に対して

\begin{align} d_{\rm norm}(x,y):= ||x-y|| \end{align}

と定めると、 $X$ 上に距離 $d_{\rm norm}:X\times X\rightarrow {\bf R}$ が導入できます。この $d_{\rm norm}$ をノルム $||\cdot||$ から定められる距離と言います。

距離はノルムを使わなくても定義できます。例えば、集合 $X$ の任意の2点 $x,y$ に対して

\begin{align} d(x,y) := \begin{cases} 1\quad (x\neq y\, のとき) \\ 0 \quad (x=y\, のとき) \end{cases} \end{align}

と定めると、 $X$ 上に距離 $d$ が導入できますが、これはあるノルムから定められた距離ではありません。というのは、ノルムから定められる距離 $d_{\rm norm}$ はノルムの定義の(3)から $d_{\rm norm}(\lambda x, \lambda y) = |\lambda| d_{\rm norm}(x,y)$ を満たしますが、上の距離 $d$ はこれを満たさないからです。

距離空間 $(X,d)$ に対しては次のように距離 $d$ を用いて位相を導入できます。

f:id:ogyahogya:20180513140632p:plain

f:id:ogyahogya:20180513140643p:plain

f:id:ogyahogya:20180513140820p:plain

工学系の研究論文は、距離空間を扱うことが多く、位相について何も言わなくても暗黙のうちに上の意味での距離 $d$ から定められる位相 $\mathcal{O}_d$ が導入されていることがほとんどです。この記事でも距離空間の位相は上のように距離 $d$ から定義された位相を導入するとします。

線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ は多様体になる

線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ が多様体になることを確認しておきます。

まず、線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ は次の和(+)とスカラー倍( $\cdot$ )でベクトル空間になります。

\begin{align} (A_1,B_1,C_1)+(A_2,B_2,C_2) &:=(A_1+A_2,B_1+B_2,C_1+C_2) \\ \lambda \cdot (A,B,C)&:= (\lambda A,\lambda B, \lambda C) \end{align}

次に、ベクトル空間である線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ に以下のノルムを導入します。任意の $(A,B,C)\in {\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ に対して、

\begin{align} ||(A,B,C)||:= ||A||_F +||B||_F + ||C||_F \end{align}

この $||\cdot ||_F$ はフロベニウスノルムと呼ばれており、任意のサイズの行列 $M$ に対して

\begin{align} ||M||_F:= \sqrt{\sum_{i,j} m_{ij}^2} \end{align}

というように各成分の２乗和の平方根で定義され、行列の集合のノルムとしてよく使われています。上で説明したように、線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ には

\begin{align} d( (A_1,B_1,C_1), (A_2,B_2,C_2) ) &:= ||(A_1,B_1,C_1)- (A_2,B_2,C_2)|| \\ &=||A_1-A_2||_F + ||B_1-B_2||_F+||C_1-C_2||_F\end{align}

というように上のノルムを用いて距離が定義できて距離空間となります。したがって、上で説明したように、この距離 $d$ を用いて線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ には位相 $\mathcal{O}_d$ が定義でき位相空間になります。

この位相空間 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ がハウスドルフ空間 (任意の異なる2点が開集合で分離できる位相空間) になることは任意の $p_1:=(A_1,B_1,C_1), p_2:=(A_2,B_2,C_2)\in {\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ とすると、 $B(p_1; d(p_1,p_2)/2), B(p_2; d(p_1,p_2)/2) \in \mathcal{O}_d$ で

\begin{align} B(p_1; d(p_1,p_2)/2)\cap B(p_2; d(p_1,p_2)/2) =\emptyset \end{align}

となることから分かります。さらに、 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ の座標近傍は ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ 自身と行列の3つ組をベクトル化することと同値な微分同相写像 $\psi: {\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}\rightarrow {\bf R}^{n(n+m+p)}$ の組だけの集合 $\{({\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}, \psi)\}$ であり、これは $C^{\infty}$ 級の座標近傍系なので線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ は多様体になります。

可制御可観測な線形システム全体の集合 $S_c\cap S_o$ は多様体になる

可制御可観測な線形システム全体の集合 $S_c\cap S_o$ が多様体になることは $S_c$ と $S_o$ が線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ の開集合であることから従います。なぜなら、 $S_c$ と $S_o$ が多様体 ${\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}$ の開集合なら位相空間の定義から $S_c\cap S_o$ も開集合となり、多様体の開集合は上で説明したように多様体になるからです。

可制御な線形システム全体の集合 $S_c$ が開集合であることを示します。可観測な線形システム全体の集合 $S_o$ が開集合であることも同様に示せます。 $S_c$ が開集合であるとは、任意の $(A,B,C)\in S_c$ とすると、ある $\delta ＞ 0$ が存在して、

\begin{align} B( (A,B,C); \delta ) \subset S_c \end{align}

となることを意味しているので、これを示します。 $(A,B,C)\in S_c$ だとすると

\begin{align} \begin{pmatrix} B & AB & \cdots & A^{n-1} B \end{pmatrix} \end{align}

の行ベクトル $v_1, v_2, \ldots, v_n\in {\bf R}^{1\times nm}$ は1次独立です。つまり、

\begin{align} (\lambda_1, \lambda_2, \ldots, \lambda_n)\neq (0, 0, \ldots, 0)\quad \Rightarrow \quad \lambda_1 v_1+\lambda_2v_2+\cdots +\lambda_n v_n \neq 0 \end{align}

が成り立ちます。行ベクトル $v_1,v_2\ldots, v_n$ たちは行列 $A, B$ の各成分の和と積から構成されているので、

\begin{align} f_i: {\bf R}^{n\times n}\times {\bf R}^{n\times m}\times {\bf R}^{p\times n}\rightarrow {\bf R}^{1\times nm}\quad {\rm s.t.}\quad v_i= f_i(A,B,C) \end{align}

となる連続関数 $f_1,f_2,\ldots, f_n$ が存在します。連続関数のスカラー倍や連続関数同士の和も連続関数となるので、

\begin{align} f(A,B,C):= \lambda_1 f_1(A,B,C) + \lambda_2 f_2(A,B,C)+\cdots +\lambda_n f_n(A,B,C) \end{align}

で定義された関数 $f$ も連続関数となります。よって、任意の $\epsilon＞0$ に対して、ある $\delta＞0$ が存在して、 \begin{align} d( (A,B,C), (A',B',C') ) < \delta \Rightarrow \| f(A,B,C)-f(A',B',C') \|_F < \epsilon \end{align}

となります。これより、 $\epsilon$ を適切に選ぶことで

\begin{align} \begin{pmatrix} B' & A'B' & \cdots & A'^{n-1} B' \end{pmatrix} \end{align}

の行ベクトルが１次独立にできることを示せれば、ある $\delta ＞ 0$ が存在して、

\begin{align} B( (A,B,C); \delta ) \subset S_c \end{align}

となることが示せたことになります。いま $v_1,v_2,\ldots, v_n$ が1次独立だったことに注意すると、

\begin{align} \exists \alpha＞0\quad {\rm s.t.}\quad ||f(A,B,C)||_F=\alpha \end{align}

が成り立ちます。さらに、\begin{align} | \|f(A',B',C') \|_F - \|f(A,B,C)\|_F | < \| f(A,B,C)-f(A',B',C') \|_F < \epsilon \end{align}

が成り立つので、\begin{align} \| f(A,B,C)\|_F -\epsilon< \| f(A',B',C') \|_F<\| f(A,B,C)\|_F +\epsilon \end{align}

が成り立ちます。よって、\begin{align} \alpha -\epsilon< \| f(A',B',C') \|_F<\alpha +\epsilon \end{align} が成り立つので、 $\epsilon$ を $0＜\epsilon＜\alpha$ が成り立つように選ぶと、\begin{align} \| f(A',B',C')\|_F \neq 0 \end{align} が成り立ちます。これは、関数 $f$ の作り方から、\begin{align} \begin{pmatrix} B' & A'B' & \cdots & A'^{n-1} B' \end{pmatrix} \end{align} の行ベクトルが１次独立であることを意味しているので、証明完了です。

このように可制御可観測な線形システムの全体の集合 $S_c\cap S_o$ は多様体になって、特に開集合であることが分かりました。 $S_c\cap S_o$ は開集合なので、ある一つの可制御可観測な線形システムの十分小さな近傍内のすべての線形システムは可制御可観測だということになります。実は、可制御可観測な線形システムの全体の集合 $S_c\cap S_o$ は線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m} \times {\bf R}^{p\times n}$ の中で稠密であることも示すことができます。つまり、線形システム全体の集合 ${\bf R}^{n\times n}\times {\bf R}^{n\times m} \times {\bf R}^{p\times n}$ の任意の開集合を $U$ としたとき、 $(S_c\cap S_o)\cap U\neq \emptyset$ が成り立ちます。このことより、一様にランダムに線形システム $(A,B,C)$ を生成すると、可制御可観測なシステムが得られる可能性が極めて高いということや、可制御可観測でないシステムに微小な摂動が加わっただけで可制御可観測なシステムになるということが言えます。言い換えると、可制御可観測なシステムであっても、実際上はほとんど可制御可観測とは言えないようなシステムもあるということであり、

で紹介したような指標を使って可制御性や可観測性の大きさを測ることが応用上は重要だということになります。

参考文献

記事を書くにあたって以下の本を参考にしました。

(1) 距離空間とかの部分を参考にしました。

集合・位相入門

作者: 松坂和夫
出版社/メーカー: 岩波書店
発売日: 1968/06/10
メディア: 単行本
購入: 12人クリック: 62回
この商品を含むブログ (56件) を見る

(2) 第2章の定理11の可制御な線形システム全体の集合は開集合であることの証明を参考にしました。

現代の最適制御理論〈上〉 (1974年) (数学叢書〈23〉)

作者: L.マーカス,占部実,古川長太,伊藤久美子
出版社/メーカー: 吉岡書店東京丸善
発売日: 1974
メディア: ?
この商品を含むブログを見る

2018-02-20

arXivについて

雑記

この記事はarXivについてのメモです。

arXivとは

arXivとは理数系の英語の原稿を無料で読むことができる以下のウェブサイトです。

https://arxiv.org/

arXivは以下のような特徴を持ちます。

1) arXivは査読がないので誰でもすぐに掲載できます。一方で、専門誌は査読があるので投稿から掲載まで1年以上かかったりします。

2) 多くの場合、専門誌に投稿した論文もarXivに投稿できます。

arXivへの投稿は研究成果の優先権の確保になるので、似たようなことをやってる人がいそうな場合はできるだけarXivに投稿した方が良さそうです。最新の成果も速く広まりますし。

arXivに画像入りの論文を投稿するときの注意

arXivには、基本的にtexのソースファイルを投稿することになりますが、画像を入れるときは、画像のファイルも投稿することになります。texや画像のファイルを投稿すると、arXivがpdfを作ってくれますが、texで

\usepackage[dvipdfmx]{graphicx}

を使っていたらダメ（arXivが生成したpdfは、本来あるはずの画像が消滅した内容になっている）です。

画像を載せるためには、\usepackage[dvipdfmx]{graphicx}ではなく

\usepackage[pdftex]{graphicx}

を使うと良いです（TeXworksだと、これ使うとエラーになんだけど、arXivでは問題なくpdfを作ってくれる）。また、elsevierが提供しているクラスファイルelsarticle.clsを利用している場合、texの中に\usepackage[dvipdfmx]{graphicx}等の図を載せるための命令を書かなくても図は載るので、tex中に\usepackage[dvipdfmx]{graphicx}等の命令を書かないでarXivにelsarticle.clsもアップロードすれば画像を載せることができます。

arXivに投稿する画像の入ったepsファイルのサイズを小さくする方法

arXivにアップロードするファイルは一つあたり400KB、合計で1MBの制限があるらしく、制限に引っかかるためepsファイルをアップロードできないことが多いです。そのため、epsファイルのサイズを小さくする方法を知っていると良いです。小さくする方法はたくさんありますが、一例をあげると、epsファイルを「GSview」で開き、Fileからconvertを選びResolution(解像度)を300に設定すると、見た目はほとんど変わらないで良い感じに小さくなったepsファイルが出来上がります。また、epsファイルをpdfファイルに変換してもサイズが小さくなります。epsファイルのpdfファイルへの変換はオンライン上でもできて、例えば以下でもできます。

convertio.co

例

arXivに投稿した自分の原稿はこんな感じです。

1) 与えられた行列に最も近いグラフラプラシアンを構成する方法を提案した論文。

[1802.06482] Optimal Graph Laplacian

2) 多数のリーダーがいるマルチエージェントシステムのリーダーの数を減らしたときに、リーダーを入れ替えた方がもとのシステムのパフォーマンスに近くなるか否かを論じた論文。結論は、入れ替えない方が近い。

[1802.06479] Optimal leader selection and demotion in leader-follower multi-agent systems

3) 安定な線形システムの $H^2$ ノルムの意味では最強のモデル低次元化法を提案した論文。 $H^2$ ノルムの意味ではこれより良い結果になる方法はない。

[1803.01523] Riemannian optimal model reduction of stable linear systems

4) IEEE Transactions on Automatic Control に掲載された論文。このように掲載済みの論文も出版社がarXivにアップすることを許している場合がある（出版社によっては掲載してから２年たったら許すとかがある）。ここで説明している正定値対称行列全体の集合のリーマン幾何学は制御理論だけでなく、信号処理、情報幾何などの研究でも役立つと思われる。

[1803.04097] Structure-preserving $H^2$ optimal model reduction based on Riemannian trust-region method

5) リーマン多様体上の最適化手法を用いた状態遷移行列に対称性をもつ連続時間線形システムの同定法

[1804.10379] Riemannian optimal system identification method of linear continuous-time systems with symmetry

2018-01-26

グラフラプラシアン

グラフ理論とか

この記事では、電力網のネットワーク、交通網のネットワーク、人間関係のネットワーク、神経ネットワーク、遺伝子ネットワークのようなネットワークシステムの性質を解析する際に重要なグラフラプラシアンについて解説します（枝に向きのないネットワークだけ解説します）。

グラフ、隣接行列、次数行列

下図のように節点と枝から構成されるネットワークを数学的に表現するには、グラフという概念が役立ちます。

グラフとは、節点の集合 $V$ と枝の集合 $E\subset V\times V$ の組 $(V, E)$ のことです。例えば上のネットワークだと節点集合が $\{1, 2, 3, 4, 5, 6\}$ で枝集合が $\{(1,2), (1,3), (1,5), (2,1),(2,3), (3,1), (3,2), (3,6), (4,5), (4,6), (5,1), (5,4), (6,3),(6,4)\}$ です。このように $(i,j)\in E$ なら $(j,i)\in E$ が成り立つグラフを正確には無向グラフといいます（この記事では、無向グラフだけを説明します）。節点は頂点、枝は辺とも呼ばれます。

上のグラフはそれぞれの枝が同等の重要度を持っているとすると、次のような行列で表現できます。

つまり、上のグラフは枝 $(1,2)$ が存在するので行列の $(1,2)$ 成分と $(2,1)$ 成分のところを1として、枝 $(1,4)$ は存在しないので行列の $(1,4)$ 成分と $(4,1)$ 成分を0としています。他の成分も同様に定義すると、上のような対称行列が得られます。逆に上の対称行列から上のグラフが構成できます。もしグラフのそれぞれの枝で重要度が異なる（重み付きグラフという）なら、重要度に応じて1以外の正の実数を行列の成分に定義することで重み付きグラフを対称行列によって表現できます。

より数学的には、グラフ $G=(V,E)$ の枝と重み（重要度）との対応関係を表す関数 $w:E\rightarrow {\bf R}_+$ をグラフ $G$ 上の重み関数といいます( ${\bf R}_+$ は正の実数の集合)。ただし、 $(i,j)\in E$ なら $w ( (i,j) )=w( (j,i) )$ を満たします。この重み付きグラフを $(G,w)$ と書きます。このような重み付きグラフ $(G,w)$ を行列によって表現した対称行列 $A=(a_{ij})$ を隣接行列といい、次のように定義されます。

f:id:ogyahogya:20180123162009p:plain

通常のグラフ $G$ はすべての $(i,j)\in E$ に対して $w( (i,j) )=1$ となる重み付きグラフのことです。

重み付きグラフ $(G,w)$ が与えられたとき、節点 $i\in V$ の次数を $d_i:= \sum_{j=1}^n a_{ij}$ で定義します。ただし、 $n$ はそのグラフの節点の数です。次のように各節点の次数を対角上に並べた対角行列を重み付きグラフ $(G,w)$ の次数行列と呼びます。\begin{align}
D={\rm diag}( d_1, d_2,\ldots, d_n)
\end{align}

例えば上のグラフが与えられたとすると、次数行列は $D={\rm diag} (3, 2, 3, 2, 2, 2)$ です。

グラフラプラシアン、接続行列

重み付きグラフ $(G,w)$ が与えられたとします。このとき、隣接行列 $A$ と次数行列 $D$ を用いて定義される次の行列 $L\in {\bf R}^{n\times n}$ を重み付きグラフ $(G,w)$ のグラフラプラシアンといいます。\begin{align}
L:=D-A
\end{align} グラフラプラシアンの定義からグラフラプラシアン $L\in {\bf R}^{n\times n}$ は固有値 0 を少なくとも一つもち、1をすべての成分に並べたベクトル ${\bf 1}\in {\bf R}^n$ は固有値 0 に対応する固有ベクトルとなります。なぜなら、 $L{\bf 1} = 0$ となるからです。

グラフラプラシアンは接続行列と呼ばれる行列によっても特徴付けられます。接続行列を定義するために、重み付き無向グラフ $(G,w)$ の枝に一意な番号と任意に向きを付けます。例えば、はじめに与えた例のグラフの枝に次のように番号と向きをつけてみましょう。

f:id:ogyahogya:20180124100158p:plain

この枝に番号と向きの付いたグラフは次の行列と対応付けることができます。

つまり、枝が節点から出るときに 1 を、枝が節点に入るときに -1 を、どちらでもないときに 0 を割り当てるというルールで枝に番号と向きの付いたグラフを行列表現できます。このような行列を接続行列といい、正確には次のように定義されます。

重み付き無向グラフの接続行列を定義したとき、各枝には番号が付いています。したがって、枝の上に定義されていた重み関数 $w:E\rightarrow {\bf R}_+$ から各枝に付いている番号の上の関数 $w_e$ が $w_e = w( (i,j) ) = w( (j,i) )$ ( $e= (i,j) = (j,i)$ ) というように定義されます。

重み付き無向グラフ $(G,w)$ のグラフラプラシアン $L\in {\bf R}^{n\times n}$ と接続行列 $B\in {\bf R}^{n\times m}$ は、枝の番号の上に定義された関数 $w_e$ を用いて次の関係で結ばれます。

この等式は枝に任意に向きを付けた接続行列から、グラフラプラシアンという枝の向きに独立な行列を構築できることを意味しています。また、この等式から重み付き無向グラフ $(G,w)$ のグラフラプラシアンは固有値 0 を少なくとも一つ持つ半正定値対称行列であることが即座に分かります。この等式の証明は以下の通りです。

f:id:ogyahogya:20180124164228p:plain

f:id:ogyahogya:20180124164354p:plain

f:id:ogyahogya:20180124164514p:plain

グラフラプラシアンという名前の由来

グラフラプラシアンという名前はベクトル解析で登場するラプラシアンという作用素に由来します。ここでは重みのない無向グラフ $(G,w)$ 、つまり任意の $(i,j)\in E$ に対して $w( (i,j) ) =1$ となる無向グラフのグラフラプラシアンがどのようにラプラシアンに関係しているかを説明します。

まず、グラフラプラシアン $L\in {\bf R}^{n\times n}$ のベクトル $\psi\in {\bf R}^n$ への作用は次のようになることに注意します。

次に、以下のような向きと重みのない格子グラフを考えます。

この格子グラフの節点を格子座標 $(x,y)$ で表し、節点の集合を $V$ とすると、上図のようにグラフの境界から離れた節点の次数は4なので、グラフラプラシアン $L$ の節点集合上の関数 $\phi: V\rightarrow {\bf R}$ への作用は、グラフラプラシアンのベクトルへの作用のアナロジーで、\begin{align}
(L\phi)(x,y) = 4\phi (x,y) -\phi(x+1,y) -\phi(x-1,y)-\phi(x,y+1)-\phi(x,y-1)
\end{align}となると考えられます。

さらに今度は同じ記号を使ってユークリッド空間 ${\bf R}^2$ 上のテイラー展開可能な関数 $\phi: {\bf R}^2\rightarrow {\bf R}$ を考えます。このとき、 $h$ を非常に小さいとして、 $\phi(x+h,y)$ と $\phi(x,y+h)$ を点 $(x,y)\in {\bf R}^2$ の周りでテイラー展開すると、\begin{align}
\phi(x+h,y) &= \phi(x,y) + h\frac{\partial \phi}{\partial x}(x,y) + \frac{h^2}{2!}\frac{\partial^2 \phi}{\partial x^2}(x,y)+\cdots \\ \phi(x,y+h) &= \phi(x,y) + h\frac{\partial \phi}{\partial y}(x,y) + \frac{h^2}{2!}\frac{\partial^2 \phi}{\partial y^2}(x,y)+\cdots
\end{align}となるので、\begin{align}
\frac{\partial^2 \phi}{\partial x^2}(x,y) &\approx \frac{\phi(x+h,y)+\phi(x-h,y)-2\phi(x,y)}{h^2} \\ \frac{\partial^2 \phi}{\partial y^2}(x,y) &\approx \frac{\phi(x,y+h)+\phi(x,y-h)-2\phi(x,y)}{h^2}
\end{align}となります。よって、\begin{align}
\frac{\partial^2 \phi}{\partial x^2}(x,y)+ \frac{\partial^2 \phi}{\partial y^2}(x,y) \approx \frac{\phi(x+h,y)+\phi(x-h,y)+\phi(x,y+h)+\phi(x,y-h)-4\phi(x,y)}{h^2}
\end{align}となりますが、 ${\bf R}^2$ 上のラプラシアン $\Delta$ は $\Delta= \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ なので、\begin{align}
(\Delta \phi)(x,y) \approx \frac{\phi(x+h,y)+\phi(x-h,y)+\phi(x,y+h)+\phi(x,y-h)-4\phi(x,y)}{h^2}
\end{align}となります。したがって、グラフラプラシアン $L$ はラプラシアン $\Delta$ の(-1)倍の離散化とみなせます。

グラフラプラシアンの名前の由来のもう少し深いところ

上で説明したように、グラフ理論のグラフラプラシアンとベクトル解析のラプラシアンには密接な関係があります。ここでは先ほどと同様に重みのない無向グラフを対象として、もう少し深いつながりを見たいと思います。

ベクトル解析のラプラシアン $\Delta$ は\begin{align} \Delta = {\rm div} ({\rm grad}) \end{align}のことです。ここで、 ${\rm grad}$ は勾配(gradient)の作用素、 ${\rm div}$ は発散(divergence)の作用素を表します。つまり、ラプラシアンは勾配と発散の合成によって表現されます。

このことと同様なことがグラフラプラシアンでも言えて、接続行列の転置が勾配に対応し、接続行列が発散に対応します。つまり、\begin{align} L = BB^T \end{align} が成り立ちます。この関係式は上で既に証明しています（今は重みなしで考えているので $w_1=w_2=\cdots =w_m=1$ であることに注意）。つまり、上で証明したグラフラプラシアンと接続行列の関係式は、ベクトル解析のラプラシアンが勾配と発散によって関係づけられていることから自然な関係式なのだということが分かります。

まず、接続行列の転置 $B^T\in {\bf R}^{m\times n}$ が勾配 ${\rm grad}$ に対応することを示します。 $B^T$ を $\phi = (\phi_1,\phi_2,\ldots, \phi_n)\in {\bf R}^n$ に作用させると、第k成分は\begin{align} (B^T\phi)_k = \sum_{i=1}^n B_{ik} \psi_i \end{align} となります。接続行列 $B$ の行は節点、列は枝に対応していたので、 $k=(i,j)\in E$ で節点iから節点jに行く向きが付いてると $B$ の第k列の第i成分は1、第j成分は-1になります。よって、\begin{align} (B^T\phi)_k = \phi_i-\phi_j \end{align}となり、節点集合上に定義された関数の差、つまり、勾配として $(B^T\phi)_k$ を見ることができます。さらに、 $B^T\phi$ 自身は枝上に定義された関数として解釈できます。したがって、ベクトル解析の勾配はスカラー場（関数）をベクトル場（関数）にしましたが、グラフ上の勾配（接続行列の転置）は節点集合上の関数を枝集合上の関数にすると解釈できます。

次に、接続行列 $B\in {\bf R}^{n\times m}$ が発散 ${\rm div}$ に対応することを示します。 $\psi = (\psi_1,\psi_2,\ldots, \psi_m)\in {\bf R}^m$ として $B$ を作用させると、\begin{align} (B\psi)_i = \sum_{k=1}^m B_{ik}\psi_k \end{align} となりますが、接続行列の定義から、これは\begin{align} (B\psi)_i = \sum_{節点iから出る枝k} \psi_k -\sum_{節点iに入る枝k}\psi_k \end{align}を意味しています。例えば、以下のような感じです。

これは、ベクトル解析の発散のイメージそのものです。ベクトル解析の発散のイメージについては、

ＥＭＡＮの物理学・電磁気学・ガウスの定理の証明

あたりを参考にしてください。

以上をまとめると、以下の対応関係が得られます。

f:id:ogyahogya:20180125165128p:plain

先ほどまで、勾配、発散、ラプラシアンはユークリッド空間上の関数に作用させていましたが、リーマン多様体上の関数に対しても勾配、発散、ラプラシアンが定義でき、上のような対応関係が得られます。リーマン多様体については、以下を参考にしてください。

ogyahogya.hatenablog.co

参考文献

記事を書くにあたって参考にした文献です。

(1) 重み付きグラフとグラフラプラシアンの定義のところを参考にしました。グラフ理論の初歩的な内容の解説と、初歩的なグラフ理論の制御工学の問題への応用が分かりやすく説明してあります。

マルチエージェントシステムの制御 (システム制御工学シリーズ)

作者: 東　俊一,永原　正章,石井　秀明,林　直樹,桜間　一徳,畑中　健志
出版社/メーカー: コロナ社
発売日: 2015/08/28
メディア: 単行本（ソフトカバー）
この商品を含むブログを見る

(2) ベクトル解析の勾配、発散、ラプラシアンがグラフ理論の接続行列の転置、接続行列、グラフラプラシアンにそれぞれ対応しているという解説を参考にしました。この本は電気回路のネットワークの性質をグラフ理論の概念を使って色々解析しています。面白いです。

線形代数とネットワーク

作者: ??金久
出版社/メーカー: 日本評論社
発売日: 2017/03/08
メディア: 単行本
この商品を含むブログを見る

(3) リーマン多様体上の関数への勾配、発散、ラプラシアンの作用素が詳しく説明されています。

ラプラシアンの幾何と有限要素法 (朝倉数学大系)

作者: 浦川肇
出版社/メーカー: 朝倉書店
発売日: 2009/10/01
メディア: 単行本
クリック: 5回
この商品を含むブログを見る

2017-11-09

最小平均二乗誤差推定値は条件付き期待値

確率とか

この記事では、下図のように $y$ を観測してパラメータ $x$ を推定しようとしたとき、推定したいパラメータ $x$ と推定値 $\hat{x}$ の二乗誤差 $||x-\hat{x}||^2$ の期待値が最小となる推定値（最小平均二乗誤差推定値） $\hat{x}$ は条件付き期待値 $E(x|y)$ で与えられ、この推定値は不偏推定値になることを説明します。 f:id:ogyahogya:20171108130548p:plain

ここでは、 $x\in {\bf R}^n$ , $y\in {\bf R}^p$ , $\hat{x}\in {\bf R}^n$ として、 $||\cdot||$ をユークリッドノルムとします。

推定したいパラメータ $x$ と観測値 $y$ は確率変数だと考える。そうすると、推定値 $\hat{x}$ も確率変数だと考えるのが自然になる。

この記事の中では、推定したいパラメータ $x$ と観測値 $y$ は確率変数だと考えます。そのとき、推定値 $\hat{x}$ も確率変数だと考えるのが自然になります。これらの理由を以下で説明します。

推定したいパラメータ $x$ は全知全能の神様にとっては確定値かもしれませんが、「 $x$ は $\hat{x}$ か？」と思っている人にとっては確定値ではなく、確率的にしか決められないものかもしれません。推定したいパラメータ $x$ は確率的に定まるものだと考えることは、数学的には、 $x$ は確率変数だと考えることに相当します。確率変数は下の記事で紹介しています。

ogyahogya.hatenablog.com

推定したいパラメータ $x$ が確率変数だと考えると、観測値 $y$ も確率変数だと考えることが自然なことが多いです。実際に、観測値は $y=x+v$ というように、確率変数である $x$ に確率変数である雑音 $v$ が加算されたものかもしれず、もしそうなら観測値 $y$ は確率変数となるためです。他にも $y$ が $x$ の関数として表現される多くの場合に、 $y$ は確率変数ということになります。

推定値 $\hat{x}$ は観測値 $y$ の関数だと考えて、関数 $f$ を用いて $\hat{x} = f(y)$ という関係があると考えるのが自然です。なぜなら、推定値 $\hat{x}$ は観測値 $y$ を参考にして決定するはずだからです。したがって、観測値 $y$ が確率変数なら $\hat{x}$ も確率変数だと考えることができます(関数 $f$ が可測関数という連続関数を含む応用上とても広いクラスの関数なら、 $y$ が確率変数のときに $f(y)$ も確率変数となる)。

なぜ $E(||x-\hat{x}||^2)$ を最小化するのか

この記事の冒頭で、二乗誤差 $||x-\hat{x}||^2$ の「期待値」を最小にする $\hat{x}$ は条件付き期待値 $E(x|y)$ だということを紹介すると書きましたが、なぜ二乗誤差 $||x-\hat{x}||^2$ そのものではなく、その「期待値」を最小にする $\hat{x}$ を考えるのでしょうか？その理由は、二乗誤差 $||x-\hat{x}||^2$ は $\hat{x}=x$ のときに最小値0となりますが、ここでは $x$ を確率変数だと考えるため $\hat{x}=x$ とすることができない上に、 $x$ が不明なので、どうやって $||x-\hat{x}||^2$ を小さくできるかも分からないためです。そのため、代わりに二乗誤差の期待値 $E(||x-\hat{x}||^2)$ の最小化を考えるのです。

$E(||x-\hat{x}||^2)$ を最小化するとは、もっと正確にはどういうことなのか

二乗誤差の期待値 $E(||x-\hat{x}||^2)$ は関係式 $\hat{x}=f(y)$ の関数 $f$ が与えられたら、次のように計算できます。

\begin{align*} E(||x-\hat{x}||^2) &= E(||x-f(y)||^2) \\ &= \int_{ {\bf R}^p } \int_{ {\bf R}^n} ||x-f(y)||^2 p(x,y) dx dy \end{align*}

この式を見れば分かるように、二乗誤差の期待値 $E(||x-\hat{x}||^2)$ は同時確率分布 $p(x,y)$ が与えられたら関数 $f$ の関数なので、 $E(||x-\hat{x}||^2)$ を

\begin{align*} R(f) := E(||x-\hat{x}||^2) \end{align*}

と書けます。このように、二乗誤差の期待値 $E(||x-\hat{x}||^2)$ を最小化するということは、同時確率分布 $p(x,y)$ が任意に与えられたときに関数 $R(f)$ を最小化する関数 $f$ を求めることを意味しています。

※関数 $R$ は関数 $f$ の関数だということから、汎関数(関数の関数)だということになります。

汎関数 $R(f)$ を最小にする関数 $f(y)$ は条件付き期待値 $E(x|y)$ になる

二乗誤差の期待値 $E(||x-\hat{x}||^2)$ 、つまり汎関数 $R(f)$ の最小化を考えてみましょう。まず、 $R(f)$ は

\begin{align*} R(f) &= \int_{ {\bf R}^p } \int_{ {\bf R}^n} ||x-f(y)||^2 p(x,y) dx dy \\ &=\int_{ {\bf R}^p } R_c(f) p(y) dy \end{align*}

となることに注意します。ここで、

\begin{align*} R_c(f):= \int_{ {\bf R}^n} ||x-f(y)||^2 p(x|y) dx =E(||x-f(y)||^2 | y) \end{align*}

です。

実は、 $R(f)$ を最小にする $f$ と $R_c(f)$ を最小にする $f$ は一致します。したがって、 $R(f)$ を最小化する $f$ の代わりに、 $R_c(f)$ を最小化する $f$ を求めれば良いということになります。

実際に、 $R_c(f)$ を計算してみると以下のようになります。\begin{align*} R_c(f) &= E(||x-f(y)||^2 | y) \\ &= E(|| (x-E(x|y)) + (E(x|y)-f(y))||^2 | y) \\ &= E(||x-E(x|y)||^2 | y) + ||E(x|y)-f(y)||^2 \\ &\geq E(||x-E(x|y)||^2 | y) \end{align*}

最後の不等式で等号が成り立つのは $f(y) = E(x|y)$ のときなので、 $R_c(f)$ を最小にする関数 $f(y)$ は条件付き期待値 $E(x|y):= \int_{{\bf R}^n} x p(x|y) dx$ だということが分かります。よって、最小平均二乗誤差推定値は $\hat{x} = E(x|y)$ となります。

※上の証明の中で突然 $E(x|y)$ が出てきました。これは、最小平均二乗誤差推定値は $\hat{x} = E(x|y)$ となることを知っていないと思いつかないかもしれません。このことを知らないとして、関数 $R_c$ を最小にする関数 $f$ を求めるためには微分 $\frac{\partial R_c}{\partial f}=0$ を満たす $f$ を求めれば良さそうです。実際に、このへんの話はカルマンフィルタ関係の本によく載っているのですが、多くの本のなかで微分 $\frac{\partial R_c}{\partial f}=0$ を満たす $f$ を求めています。が、実はこの微分といっているものは通常の微分だと考えることができません。なぜなら $f$ は数値ではなく関数だからです。ということで、微分 $\frac{\partial R_c}{\partial f}$ を計算するといっても通常の意味の微分は計算できません。では、どうするかというと汎関数 $R_c$ の変分(汎関数微分)を計算し、その変分がゼロになる関数 $f$ が答えだということになります。つまり、力学でオイラー・ラグランジュ方程式を導出する際にラグランジアンの作用積分の変分を計算したように、汎関数 $R_c$ の変分を計算すれば良いということになります。では実際に汎関数 $R_c$ の変分を計算してみましょう。汎関数 $R_c$ の点 $f$ における $h$ 方向への微分は

\begin{align*} \delta R_c(f)[h] := \lim_{\epsilon\rightarrow 0} \frac{R_c (f+\epsilon h) -R_c(f)}{\epsilon} \end{align*} となります（上記の $\delta R_c(f)$ がゼロとなる関数 $f$ を見つけることを変分問題 $\delta R_c(f)=0$ と書く）。汎関数 $R_c(f)$ の定義から、

\begin{align*} R_c(f+\epsilon h) -R_c(f) = \int_{{\bf R}^n} \left( 2\epsilon (x-f(y))^T h(y) + \epsilon^2 ||h(y)||^2 \right) p(x|y) dx \end{align*}

となります。よって、

\begin{align*} \delta R_c(f)[h] &= 2 \int_{ {\bf R}^n} (x-f(y))^T h(y) p(x|y) dx \\ &= 2(E(x|y)-f(y))^T h(y) \end{align*}

ということになります。したがって、 $f(y)=E(x|y)$ のとき $\delta R_c(f)=0$ となることが分かりました。

最小平均二乗誤差推定値 $\hat{x} = E(x|y)$ は不偏推定値

一般に、推定したいパラメータ $x$ と推定値 $\hat{x}$ の誤差の期待値 $E(x-\hat{x})$ がゼロになるとき推定値 $\hat{x}$ は推定したいパラメータ $x$ の不偏推定値だと言います。つまり、推定したいパラメータの期待値と推定値の期待値が一致しているときに、その推定値は不偏だと言うのです。

最小平均二乗誤差推定値 $\hat{x}=E(x|y)$ は次のように不偏推定値だということが分かります。

\begin{align*} E(\hat{x}) &= E( E(x|y)) \\ &= \int_{{\bf R}^p} \left( \int_{{\bf R}^n} x p(x|y) dx \right) p(y) dy \\ &= \int_{{\bf R}^p} \int_{{\bf R}^n} x p(x,y) dx dy \\ &= \int_{ {\bf R}^n} x \int_{ {\bf R}^p} p(x,y) dy dx\\ &= \int_{{\bf R}^n} x p(x) dx = E(x) \end{align*}

参考文献

記事を書くにあたって参考にした文献です。

(1) 理論的なことがしっかり書いてる。

応用カルマンフィルタ

作者: 片山徹
出版社/メーカー: 朝倉書店
発売日: 2000/02/01
メディア: 単行本
クリック: 6回
この商品を含むブログを見る

この記事の中では、二乗誤差 $||x-\hat{x}||^2$ の期待値の最小化を考えましたが、この本の中では、他にも絶対誤差の期待値や一様誤差の期待値の最小化も考えられており、条件付き確率分布 $p(x|y)$ がガウス分布に従うなら、絶対誤差と一様誤差の期待値の最小値を与える関数は二乗誤差の期待値の最小値を与える条件付き期待値 $E(x|y)$ に一致することが説明されています。ということで、条件付き期待値 $E(x|y)$ は他の評価指標を使っても最も良い推定値になることがあるということが分かります。

ガウス分布は、以下の記事で説明したように中心極限定理と密接に関連した確率分布で応用上よくでてきます。

ogyahogya.hatenablog.com(2) 上の参考文献と同じ著者の本。汎関数の変分を計算すべきところを普通の微分のように計算しているのは誤りであることに注意。この本のように英語の本でも普通の微分のように計算しているものが多々あります。。。

非線形カルマンフィルタ

作者: 片山　徹
出版社/メーカー: 朝倉書店
発売日: 2011/11/30
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

2017-07-04

可制御性グラミアンと可観測性グラミアン

制御理論

対象とする線形システム

f:id:ogyahogya:20170630161131p:plain

の可制御性と可観測性の定義は

ogyahogya.hatenablog.com

で紹介しましたが、この記事では可制御性と可観測性の「大きさ」を定量的に測る手段を紹介します。線形システムについては

ogyahogya.hatenablog.com

で紹介しています。この記事の中では行列 $A$ を安定、つまり、 $A$ のすべての固有値の実部は負であると仮定します。

入力エネルギー最小化問題

まず、線形写像 $\Psi_c: L^2(-\infty,0] \rightarrow {\bf R}^n$ を

f:id:ogyahogya:20170630161918p:plain

と定義します。ここで、 $L^2(-\infty,0]$ は区間 $(-\infty,0]$ 上で定義された2乗可積分な関数全体の線形空間です。このような空間については、

ogyahogya.hatenablog.com

でも紹介しています。この $\Psi_c(u)$ は無限の過去 $t=-\infty$ に状態が0であるときの現在 $t=0$ の状態 $x(0)$ を表していると解釈できます。実際に、 $\lim_{t\rightarrow -\infty} x(t) =0$ とすると、行列 $A$ が安定なので、\begin{align}
x(0) &= \lim_{t_0\rightarrow -\infty}\left(e^{-At_0}x(t_0)\right) + \int_{-\infty}^0 e^{A(0-t)}Bu(t)dt\\&=\int_{-\infty}^0 e^{-At}Bu(t)dt = \Psi_c(u)
\end{align}となるからです。

$v_1, v_2\in L^2(-\infty,0]$ とするとき、 $\langle v_1, v_2 \rangle_2:= \int_{-\infty}^0 v_1^T(t) v_2(t) dt$ によって $\langle \cdot, \cdot \rangle_2$ を定義し、 $||v_1||_2:=\sqrt{\langle v_1, v_1 \rangle_2}$ によって $||\cdot ||_2$ を定義します。

このとき、次の入力エネルギー最小化問題を考えます。

f:id:ogyahogya:20170630164230p:plain

ただし、 $||\cdot||$ はユークリッドノルムを表します。この問題は、「状態を原点から単位球面上に到達させることが可能な2乗可積分な入力の中で、エネルギーが最小になるものを求めよ。」ということを意味しています。

可制御性グラミアン

上の入力エネルギー最小化問題と次の可制御性グラミアン

f:id:ogyahogya:20170703161119p:plain

は密接な関係があります。ここでは、その関係を見るための準備をします。

線形作用素 $\Psi_c$ と行列 $W_c$ の間には次の関係が成り立ちます。

f:id:ogyahogya:20170703161504p:plain

ここで、 $\Psi_c^*$ は $\Psi_c$ の共役作用素です。共役作用素はこちらで説明しましたので、興味がありましたら参照してください。

ogyahogya.hatenablog.com

上の関係が成り立つことは次のように確認できます。

f:id:ogyahogya:20170703163148p:plain

f:id:ogyahogya:20170703163648p:plain

任意の $\xi\in {\bf R}^n$ に対して、上記の式が成り立つことから $W_c=\Psi_c\Psi_c^*$ が成り立つことが証明できました。

可制御性グラミアンの表式から $W_c$ は半正定値対称行列だということが分かります。実は、線形システムが可制御であることと、 $W_c$ が正定値対称行列であることは等価です。ここでは簡単のために、線形システムは可制御であるとして話を進めます。

直交射影

入力エネルギー最小化問題の解と可制御性グラミアンの関係を、すっきりと理解するために、直交射影の概念を紹介します。

f:id:ogyahogya:20170704093559p:plain

直交射影には、次の性質があります。

f:id:ogyahogya:20170704100020p:plain

入力エネルギー最小化問題の解と可制御性グラミアンの関係

入力エネルギー最小化問題の解は、次のように可制御性グラミアン $W_c$ を用いて表現できることが分かります。

f:id:ogyahogya:20170703164814p:plain

実は、ある直交射影 $P: L^2(-\infty,0] \rightarrow L^2(-\infty,0]$ が存在し、(※)を満たす任意の $u\in L^2(\infty,0]$ に対して、 $Pu=u_{最適}$ となって $Pu$ も(※)を満たすことが示せます。直交射影の定義の下で示したように、 $||Pu||_2 \leq ||u||_2$ なので、上の主張が成り立つということになります。

そのような直交射影 $P$ は

f:id:ogyahogya:20170704101524p:plain

で与えられます。 $Pu$ が(※)を満たすことは、

f:id:ogyahogya:20170704101724p:plain

から分かります。さらに、★が直交射影であることも少し計算すれば示せます。

さらに、入力エネルギーと可制御性グラミアンの間には、次の関係が成り立つことが分かります。

f:id:ogyahogya:20170704102951p:plain

なお、 $W_c$ の逆行列 $W_c^{-1}$ は、線形システムを可制御であると考えているので存在します。

可制御性グラミアンの幾何学的解釈

ここでは、入力 $u\in L^2(-\infty,0]$ ( $||u||_2\leq 1$ ) で到達可能な状態を可制御性グラミアン $W_c$ を用いて特徴付けます。この特徴付けは次の事実に基づきます。

f:id:ogyahogya:20170704111337p:plain

f:id:ogyahogya:20170704111350p:plain

上で示したことから、入力 $u\in L^2(-\infty,0]$ ( $||u||_2\leq 1$ ) で到達可能な状態は $W_c^{1/2} x$ ( $||x||\leq 1$ )だということが分かりました。

以上のことから「可制御性の大きさ」、すなわち、「制御のしやすさ」と可制御性グラミアン $W_c$ の関係を調べることができます。このことを見るために、単位エネルギー $||u||_2= 1$ の入力 $u\in L^2(-\infty,0]$ で到達可能な状態の集合は

f:id:ogyahogya:20170704122935p:plain

であることに注意しましょう。今、 $W_c^{1/2}$ の固有値を $\lambda_1,\lambda_2,\ldots, \lambda_n$ 、対応する正規直交固有ベクトルを $v_1,v_2,\ldots,v_n$ とすると、

f:id:ogyahogya:20170704123212p:plain

となり、 $\mathcal{E}_c$ は固有ベクトル $v_1,v_2,\ldots,v_n$ の方向を主軸として固有値 $\lambda_1,\lambda_2,\ldots,\lambda_n$ を主値とする楕円であることが分かります。イメージ図は下のような感じです。

f:id:ogyahogya:20170704123700p:plain

つまり、単位エネルギー $||u||_2= 1$ の入力 $u\in L^2(-\infty,0]$ で原点から最も遠くまで到達可能な方向が可制御性グラミアン $W_c$ の最大固有値に対応する固有ベクトルの方向で、2番目に遠くまで到達可能な方向が２番目に大きい固有値に対応する固有ベクトルの方向、３番目に・・・ということになります。言い換えると、最も制御しやすい方向が可制御性グラミアン $W_c$ の最大固有値に対応する固有ベクトルの方向で、最も制御しにくい方向が $W_c$ の最小固有値に対応する固有ベクトルの方向だということです。

以上から、線系システムが可制御であっても、制御のしやすい方向としにくい方向があり、それらを定量的に調べるためには可制御性グラミアン $W_c$ の固有値と固有ベクトルを調べれば良いということが分かりました。

可観測性グラミアン

次に、観測のしやすさを定量的に評価する可観測性グラミアンについて説明します。

今、現在の状態 $x(0)$ を $x_0$ とします。このとき、 $\Psi_o:{\bf R}^n\rightarrow [0,\infty)$ を

f:id:ogyahogya:20170704142459p:plain

で定義します。線形写像 $\Psi_o$ は $[0,t]$ 上で入力を $u=0$ としたときの線形システムの出力値と関係付きます。なぜなら、 $y(t) = C\exp (At) x_0 = \Psi_o x_0$ が成り立つからです。可制御性グラミアンの時と同様の計算で、 $\Psi_o$ は可観測性グラミアン

f:id:ogyahogya:20170704142729p:plain

と

f:id:ogyahogya:20170704142927p:plain

の関係があることが分かります。ここでは、 $L^2 [0,\infty)$ 上に内積を $\langle f_1, f_2 \rangle_2 := \int_0^{\infty} f_1(t)^Tf_2(t) dt$ で定め、ノルムを $||f||_2:= \sqrt{\langle f, f \rangle_2}$ で定めることにします。すると、出力は $y=\Psi_o x_0$ のため、出力エネルギーが

f:id:ogyahogya:20170704143554p:plain

と可観測性グラミアン $W_0$ を用いて表せることが分かります。可制御性グラミアンの時と同様に、可観測性グラミアン $W_o$ を用いて定義された集合

f:id:ogyahogya:20170704143756p:plain

は楕円であり、主軸は $W_o$ の固有ベクトルの方向で、主値は固有値の平方根ということになります。これは球面上の状態 $x_0$ にも「観測のしやすさ」があり、可観測性グラミアン $W_o$ の最大固有値に対応する固有ベクトルの方向の状態が最も観測しやすく、最小固有値に対応する固有ベクトルの方向の状態が最も観測しにくいという意味になります。（ここで、観測のしやすさを出力エネルギー $||y||_2^2$ の大小で解釈しています。）

参考文献

記事を書くにあたって次の本を参考にしました。

A Course in Robust Control Theory: A Convex Approach (Texts in Applied Mathematics)

作者: Geir E. Dullerud,Fernando Paganini
出版社/メーカー: Springer
発売日: 2000/02/02
メディア: ハードカバー
この商品を含むブログを見る

2017-02-06

行列の指数関数と対数関数

多様体とか

この記事では、行列の指数関数と対数関数について解説します。Lie群とLie環という概念を理解するための準備に相当します。

●行列の内積とノルム

${\bf R}^{n\times n}$ をn行n列の実数を成分に持つ行列全体の集合とします。このとき、 ${\bf R}^{n\times n}$ に以下の内積とノルムを導入できます。

f:id:ogyahogya:20170202134934p:plain

任意の $A, B\in {\bf R}^{n\times n}$ に対して $||AB||\leq ||A||\cdot ||B||$ が成り立つことがコーシー・シュワルツの不等式を使うことで確認できます。

●行列の指数関数

実数 $x$ の指数関数 $e^x$ を $x=0$ の周りでテイラー展開すると $e^x = \sum_{k=0}^{\infty} \frac{1}{k!} x^k$ でした。これと形式的に同じになるように行列の指数関数が次のように定義されます。

f:id:ogyahogya:20170202140703p:plain

右辺の無限級数が収束することは以下のように分かります。

f:id:ogyahogya:20170202140816p:plain

●行列の指数関数の性質

行列の指数関数の性質は以下の結果を利用して明らかにすることができます。

f:id:ogyahogya:20170202144600p:plain

解であることは $\gamma(t)=\exp(tX)$ を左辺と右辺に代入して等しくなることを見れば良いです。一意性は、微分方程式の解の一意性定理から従います。

指数関数の性質を明らかにする前に、二つの行列が可換であるか否かを判定するカッコ積を導入しておきます。

f:id:ogyahogya:20170202145056p:plain

つまり、 $[X,Y]=0$ なら $X$ と $Y$ の積は可換です。

次の結果が重要です。

f:id:ogyahogya:20170202151228p:plain

これを応用して、以下の結果が得られます。

f:id:ogyahogya:20170202151802p:plain

これは指数関数の値域は一般線形群 ${\rm GL}(n,{\bf R})$ であることを意味しています。

f:id:ogyahogya:20170202152621p:plain

さらに、以下の結果も成り立ちます。

f:id:ogyahogya:20170202152800p:plain

これらより、 $X\in {\bf R}^{n\times n}$ に対して、 $\gamma(t)=\exp(tX)$ は加法群 ${\bf R}$ から ${\rm GL}(n,{\bf R})$ への微分可能な準同型写像 $\gamma:{\bf R}\rightarrow {\rm GL}(n,{\bf R})$ を与えていることが分かります。この曲線 $\gamma$ は実は次のように $\dot{\gamma}(0)$ を用いて表されます。