固有値問題と行列方程式

Author: Shun Suzuki

Date: 2024-01-09

ここでは, 件の最適化問題の解法の1つとして, Longらによる方法 ¹を紹介する. 正直, この論文はお世辞にも理解しやすいとは言い難く, この節はおそらくこうだろうという自己解釈を多分に含むので注意されたい.

NOTE: この手法は現在SDKでは実装されていない.

Longらの方法は, 単一焦点解の重ね合わせに基づく. 単一焦点解の重ね合わせとして駆動した場合, 各焦点の位相は何が最適か, を固有値問題として解く.

まずはじめに, 単一焦点を生成する解について考えよう. ある点${\mathbf{r}}_i$での複素音圧は $$\begin{array}{r}p({\mathbf{r}}_i)=\sum _j{G}_{ij}{a}_j{\mathrm{e}}^{\mathrm{i}(k\Vert {\mathbf{r}}_i-{\mathbf{r}}_j\Vert -{\varphi }_j)},\end{array}$$ であるため, ${\mathbf{r}}_i$に焦点を生成する解は明らかに $$\begin{array}{r}{\varphi }_j=k\Vert {\mathbf{r}}_i-{\mathbf{r}}_j\Vert ,\end{array}$$ である. これは, 伝搬による位相遅れを補償するように, 遠くの振動子が早めに音を出すことを意味する.

これには, 時間対称性を用いたもう1つの解釈がある. (損失のない) 波動方程式は時間に対して対称であるため, 焦点に音源を置いたときの各振動子における音圧を記録し, これを逆再生することで焦点が生成できる. 逆再生するため,焦点から振動子への伝搬行列は$\overline{G_{ij}}$と表される. ここで, $\overline{\cdot }$は複素共役を表す. すなわち, $a_j{\mathrm{e}}^{-\mathrm{i}{\varphi }_j}=\overline{G_{ij}}|p({\mathbf{r}}_i)|$²とすれば, 逆に焦点が生成されるはずである. これより, $$\begin{array}{r}\sum _j{G}_{ij}{\mathrm{e}}^{\mathrm{i}(k\Vert {\mathbf{r}}_i-{\mathbf{r}}_j\Vert )}|p({\mathbf{r}}_i)|\overline{G_{ij}}{\mathrm{e}}^{-\mathrm{i}k\Vert {\mathbf{r}}_i-{\mathbf{r}}_j\Vert }=|p({\mathbf{r}}_i)|\sum _j|G_{ij}{|}^2,\end{array}$$ となり, 焦点での音圧$|p({\mathbf{r}}_i)|$からずれる. これは, 振動子の配置が離散的であることに起因する. しかし, 今考えているのはフェーズドアレイなので, この分を補正しなくてはいけない.

したがって, 実際の振動子の駆動は $$\begin{array}{r}{a}_j{\mathrm{e}}^{-\mathrm{i}{\varphi }_j}=\frac{|p_i|\overline{G_{ij}}}{{\sum }_j|G_{ij}{|}^2},\end{array}$$ とすべきである(論文 ¹の式(8)).

ここで, 次のような行列$X$を考える, $$\begin{array}{r}X=\left(\begin{array}{ccc}\frac{|p_0|\overline{G_{00}}}{{\sum }_j|G_{0j}{|}^2}& \cdots & \frac{|p_{M-1}|\overline{G_{M-1,0}}}{{\sum }_j|G_{M-1,j}{|}^2}\\ ⋮& \ddots & ⋮\\ \frac{|p_0|\overline{G_{0,N-1}}}{{\sum }_j|G_{0j}{|}^2}& \cdots & \frac{|p_{M-1}|\overline{G_{M-1,N-1}}}{{\sum }_j|G_{M-1,j}{|}^2}\end{array}\right).\end{array}$$ 少しわかりにくいが, 例えば, 第0列は${\mathbf{r}}_0$の位置に音圧$|p_0|$の焦点を生成しようとしたときの振動子の駆動ベクトル$\mathbf{q}$である. ここで, 焦点音圧の位相 ${\psi }_i$から, 位相ベクトル$\mathbf{x}$を, $$\begin{array}{r}\mathbf{x}=[{\mathrm{e}}^{-\mathrm{i}{\psi }_0},\dots ,{\mathrm{e}}^{-\mathrm{i}{\psi }_{M-1}}]{}^{\mathsf{T}},\end{array}$$ $$\begin{array}{r}X\mathbf{x},\end{array}$$ とすると, これは焦点に (位相込の) 複素音圧$\mathbf{p}$の焦点を生成しようとしたときの振動子の各駆動ベクトルの線形重ね合わせをあらわす. したがって, $GX\mathbf{x}$は, そのような線形重ね合わせ駆動ベクトルが, 実際に焦点に生成する音圧を表す. (行列$GX$が論文 ¹の行列$R$に等しい.) これが, もともとの焦点音圧$\mathbf{p}=\mathrm{diag}(|p_0|,\dots ,|p_{M-1}|)\mathbf{x}$ (の定数倍) に等しくなってほしいので $$\begin{array}{r}GX\mathbf{x}=\lambda \mathrm{diag}(|p_0|,\dots ,|p_{M-1}|)\mathbf{x},\end{array}$$ が解くべき問題になり, これはまさに固有値問題である. 定数倍の自由度は, 最終的に振動子の出力を$\lambda$分の1倍すればいい話なので問題にならない.

さて, 上の固有値問題の解のうち, もっとも望ましいのは$\lambda$のもっとも大きいものである. なぜなら, $\lambda$が大きいなら, その分だけ振動子の出力を抑えることができるからである. これは, 省電力化とアーティファクトの抑制につながる. 以上により, 焦点の位相の最適化が終わる.

焦点の位相が求まったので, 各焦点を生成する駆動ベクトルの重ね合わせである, $X\mathbf{x}$を$\mathbf{q}$とすれば良いような気もするが, この部分は振動子の出力可能な最大音圧を一切考慮していない. すでに述べたように, 一部の振動子のみが突出して強い音を出すような解は望ましくない. そこで, 振動子の出力のばらつきを抑えるためにTikhonov正則化を導入し, 最適化問題を次のように拡張したものを解く. $$\begin{array}{r}\left(\begin{array}{c}\mathbf{p}\\ 0\\ ⋮\\ 0\end{array}\right)=\left(\begin{array}{ccc}& G& \\ {\sigma }_0^{\gamma }& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\sigma }_{N-1}^{\gamma }\end{array}\right)\mathbf{q}.\end{array}$$ ここで, $$\begin{array}{r}{\sigma }_j=\sqrt{\left|\sum _{i=0}^{M-1}\frac{G_{ij}|p_i|}{M}\right|},\end{array}$$ である.

正直に言うと, ここの解釈はいまだにきちんと理解できていない. おおよそは以下のような感じであろう. まず, ${\sigma }_j$は$j$番目の振動子が持つ各焦点への“影響力“のようなものを表している. 連立方程式に${\sigma }_j{q}_j=0$のような式を追加することにより, 影響力の強い振動子の出力は弱めることにより, 振動子の出力のバランスを良くする. ただし, なぜ${\sigma }_j$がこの表記になるのかはわからない.

なお, Longらの論文 ¹では, 固有値問題や逆問題の具体的な解き方は明記されていない. また, 出力は振動子の限界を超えた場合には切り捨てるという方針を取る.

計算量は, 固有値問題を解くのに$O(M^3)$, 逆問題を解くのに, 例えばLU分解を使うと$O(N^3)$となる. したがって, $O(\max (M,N{)}^3)$である³.