非線形最小二乗法

Author: Shun Suzuki

Date: 2024-01-09

件の行列方程式を次の非線形最小二乗問題として解く方法もある, $$\begin{array}{r}\min \Vert G\mathbf{q}-\mathbf{p}{\Vert }^2.\end{array}$$ ここでのパラメータは$\mathbf{p}$の位相と$\mathbf{q}$である. 反復的に上式を最小化していく方法がいくつか知られている.

Levenberg-Marquardt (LM) 法

ここでは, LM (Levenberg-Marquardt) 法 ^{1 2}を紹介する. LM法を用いた例としては ^{3 4}などがある.

後のため, 推定すべきパラメータ$\mathbf{\theta }$を $$\begin{array}{rl}\mathbf{\varphi }& =[-{\varphi }_0,\dots ,-{\varphi }_{N-1}]{}^{\mathsf{T}},\\ \mathbf{\psi }& =[-{\psi }_0,\dots ,-{\psi }_{M-1}]{}^{\mathsf{T}},\\ \mathbf{a}& =[a_0,\dots ,a_{N-1}]{}^{\mathsf{T}},\\ \mathbf{\theta }& =[\mathbf{\varphi }{}^{\mathsf{T}},\mathbf{\psi }{}^{\mathsf{T}},\mathbf{a}{}^{\mathsf{T}}]{}^{\mathsf{T}}.\end{array}$$ と置き, 最適化問題を $$\begin{array}{r}\min \Vert \mathbf{z}(\mathbf{\theta }){\Vert }^2,\end{array}$$ $$\begin{array}{rl}\Vert \mathbf{z}(\mathbf{\theta }){\Vert }^2& =\Vert G\mathbf{q}-\mathbf{p}{\Vert }^2\\ & =\Vert G\mathbf{q}-\mathrm{diag}(|p_0|,\dots ,|p_{M-1}|){\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}{\Vert }^2\\ & =\Vert G\mathbf{q}-P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}{\Vert }^2(P=\mathrm{diag}(|p_0|,\dots ,|p_{M-1}|))\end{array}$$ と書く. ここでベクトル$\mathbf{x}=[x_0,\dots ,x_k]{}^{\mathsf{T}}$に対して${\mathrm{e}}^{\mathbf{x}}$は $$\begin{array}{r}{\mathrm{e}}^{\mathbf{x}}=[{\mathrm{e}}^{x_0},\dots ,{\mathrm{e}}^{x_k}]{}^{\mathsf{T}}\end{array}$$ を意味する. まず, $\mathbf{z}$は複素関数であることに注意する. LM法を適用できるようにするため, 実部と虚部に分けると, $$\begin{array}{r}\Vert \mathbf{z}(\mathbf{\theta }){\Vert }^2={\Vert \left(\begin{array}{c}\mathcal{R}[\mathbf{z}(\mathbf{\theta })]\\ \mathcal{I}[\mathbf{z}(\mathbf{\theta })]\end{array}\right)\Vert }^2\end{array}$$ となる. ここで, $\mathcal{R}[z],\mathcal{I}[z]$はそれぞれ$z$の実部と虚部を表す. $$\begin{array}{r}\mathbf{f}(\mathbf{\theta })=\left(\begin{array}{c}\mathcal{R}[\mathbf{z}(\mathbf{\theta })]\\ \mathcal{I}[\mathbf{z}(\mathbf{\theta })]\end{array}\right)\end{array}$$ とおくと, $\mathbf{f}$は実関数であり $$\begin{array}{r}\min \Vert G\mathbf{q}-\mathbf{p}{\Vert }^2\iff \min \Vert \mathbf{f}(\mathbf{\theta }){\Vert }^2\end{array}$$ となるので, LM法が適用できる.

この時, $\mathbf{f}(\mathbf{\theta })$のJacobian $J(\mathbf{\theta })$は $$\begin{array}{rl}(\mathbf{z}(\mathbf{\theta }){)}_k& =\sum _j{G}_{kj}{a}_j{\mathrm{e}}^{\mathrm{i}{\varphi }_j}-|p_k|{\mathrm{e}}^{\mathrm{i}{\psi }_k},\\ \therefore \frac{\partial (\mathbf{z}(\mathbf{\theta }){)}_k}{\partial {\varphi }_j}& =\mathrm{i}{G}_{kj}{a}_j{\mathrm{e}}^{\mathrm{i}{\varphi }_j},\\ \frac{\partial (\mathbf{z}(\mathbf{\theta }){)}_k}{\partial {\psi }_i}& =-\mathrm{i}{\delta }_{ki}|p_k|{\mathrm{e}}^{\mathrm{i}{\psi }_k},\\ \frac{\partial (\mathbf{z}(\mathbf{\theta }){)}_k}{\partial a_j}& =G_{kj}{\mathrm{e}}^{\mathrm{i}{\varphi }_j},\end{array}$$ から, $$\begin{array}{rl}J(\mathbf{\theta })& =\frac{\partial \mathbf{f}(\mathbf{\theta })}{\partial \mathbf{\theta }}\\ & =\left(\begin{array}{ccc}\mathcal{R}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]& \mathcal{R}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]& \mathcal{R}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]\\ \mathcal{I}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]& \mathcal{I}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]& \mathcal{I}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]\end{array}\right)\end{array}$$ となる. この$J(\mathbf{\theta })$を用いて, LM法では, 以下の更新式に従いパラメータを更新する. $$\begin{array}{r}{\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-(J({\mathbf{\theta }}_t){}^{\mathsf{T}}J({\mathbf{\theta }}_t)+\lambda (t)I{)}^{-1}J({\mathbf{\theta }}_t){}^{\mathsf{T}}f(\mathbf{\theta }).\end{array}$$ $\lambda (t)$の決め方にはいくつかのバリエーションが知られている. 例えば, 文献⁵を参照.

ここで, 計算量を削減するため, $J({\mathbf{\theta }}_t){}^{\mathsf{T}}J({\mathbf{\theta }}_t)$をあらかじめ計算しておこう. $$\begin{array}{rl}& J({\mathbf{\theta }}_t){}^{\mathsf{T}}J({\mathbf{\theta }}_t)\\ & =\left(\begin{array}{cc}\mathcal{R}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]{}^{\mathsf{T}}& \mathcal{I}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]{}^{\mathsf{T}}\\ \mathcal{R}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]{}^{\mathsf{T}}& \mathcal{I}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]{}^{\mathsf{T}}\\ \mathcal{R}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]{}^{\mathsf{T}}& \mathcal{I}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]{}^{\mathsf{T}}\end{array}\right)\\ & \times \left(\begin{array}{ccc}\mathcal{R}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]& \mathcal{R}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]& \mathcal{R}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]\\ \mathcal{I}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]& \mathcal{I}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]& \mathcal{I}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]\end{array}\right)\end{array}$$ であるが, $$\begin{array}{r}\mathcal{R}[X]{}^{\mathsf{T}}\mathcal{R}[Y]+\mathcal{I}[X]{}^{\mathsf{T}}\mathcal{I}[Y]=\mathcal{R}[X^{†}Y]\end{array}$$ であるから, 結局, $$\begin{array}{rl}& J({\mathbf{\theta }}_t){}^{\mathsf{T}}J({\mathbf{\theta }}_t)\\ & =\mathcal{R}\left(\begin{array}{ccc}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}G\mathrm{diag}(\mathbf{q})& -\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})& -\mathrm{i}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})\\ \ast & \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})& \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})\\ \ast & \ast & \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})\end{array}\right)\end{array}$$ ($J({\mathbf{\theta }}_t){}^{†}J({\mathbf{\theta }}_t)$はHermite行列なので, 一部省略した. 以下同様.)

ここで, 列ベクトル$\mathbf{x},\mathbf{y}$に対して, $$\begin{array}{rl}{\left(\mathrm{diag}(\mathbf{x})M\mathrm{diag}(\mathbf{y})\right)}_{ij}& =\sum _l\left(\sum _k{\delta }_{ik}{x}_k{M}_{kl}\right){\delta }_{lj}{y}_j\\ & =\sum _l\sum _k{\delta }_{ik}{\delta }_{lj}{x}_k{M}_{kl}{y}_j\\ & =M_{ij}{x}_i{y}_j\end{array}$$ であり, したがって, $$\begin{array}{r}\mathrm{diag}(\mathbf{x})M\mathrm{diag}(\mathbf{y})=M\circ \mathbf{x}\mathbf{y}{}^{\mathsf{T}}\end{array}$$ となる. なお, $X\circ Y$は$X$と$Y$のHadamard積を意味する. ここで, $$\begin{array}{rl}B& =\left(\begin{array}{ccc}G& -P& -\mathrm{i}G\end{array}\right),\\ T& =\left(\begin{array}{c}\overline{\mathbf{q}}\\ {\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }}\\ {\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }}\end{array}\right)\end{array}$$ と定義すると, $$\begin{array}{rl}& \left(\begin{array}{ccc}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}G\mathrm{diag}(\mathbf{q})& -\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})& -\mathrm{i}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})\\ \ast & \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})& \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})\\ \ast & \ast & \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})\end{array}\right)\\ & =\left(\begin{array}{ccc}G{}^{†}G& -G{}^{†}P& -\mathrm{i}G{}^{†}G\\ \ast & P{}^{\mathsf{T}}P& \mathrm{i}P{}^{\mathsf{T}}G\\ \ast & \ast & G{}^{†}G\end{array}\right)\circ \left(\begin{array}{ccc}\overline{\mathbf{q}}\mathbf{q}{}^{\mathsf{T}}& \overline{\mathbf{q}}\left({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\right){}^{\mathsf{T}}& \overline{\mathbf{q}}\left({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }}\right){}^{\mathsf{T}}\\ \ast & {\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }}\left({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\right){}^{\mathsf{T}}& {\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }}\left({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }}\right){}^{\mathsf{T}}\\ \ast & \ast & {\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }}\left({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }}\right){}^{\mathsf{T}}\end{array}\right)\\ & =(B{}^{†}B)\circ TT{}^{†}\end{array}$$ となり, $$\begin{array}{r}J({\mathbf{\theta }}_t){}^{\mathsf{T}}J({\mathbf{\theta }}_t)=\mathcal{R}\left[(B{}^{†}B)\circ TT{}^{†}\right]\end{array}$$ と計算できる.

次に, $J({\mathbf{\theta }}_t){}^{\mathsf{T}}f(\mathbf{\theta })$を計算しよう. $$\begin{array}{rl}J({\mathbf{\theta }}_t){}^{\mathsf{T}}f(\mathbf{\theta })& =\left(\begin{array}{cc}\mathcal{R}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]{}^{\mathsf{T}}& \mathcal{I}[\mathrm{i}G\mathrm{diag}(\mathbf{q})]{}^{\mathsf{T}}\\ \mathcal{R}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]{}^{\mathsf{T}}& \mathcal{I}[-\mathrm{i}P\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})]{}^{\mathsf{T}}\\ \mathcal{R}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]{}^{\mathsf{T}}& \mathcal{I}[G\mathrm{diag}({\mathrm{e}}^{\mathrm{i}\mathbf{\varphi }})]{}^{\mathsf{T}}\end{array}\right)\left(\begin{array}{c}\mathcal{R}[\mathbf{z}(\mathbf{\theta })]\\ \mathcal{I}[\mathbf{z}(\mathbf{\theta })]\end{array}\right)\\ & =\mathcal{R}\left[\left(\begin{array}{c}-\mathrm{i}\mathrm{diag}(\mathbf{q}){}^{†}G{}^{†}\mathbf{z}\\ \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}\mathbf{z}\\ \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}\mathbf{z}\end{array}\right)\right]\end{array}$$ ここで, $\mathbf{z}=G\mathbf{q}-P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}$なので, $$\begin{array}{rl}\left(\begin{array}{c}-\mathrm{i}\mathrm{diag}(\mathbf{q}){}^{†}G{}^{†}\mathbf{z}\\ \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}\mathbf{z}\\ \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}\mathbf{z}\end{array}\right)& =\left(\begin{array}{c}-\mathrm{i}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}(G\mathbf{q}-P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})\\ \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}(G\mathbf{q}-P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})\\ \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}(G\mathbf{q}-P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }})\end{array}\right)\\ & =-\mathrm{i}\left(\begin{array}{c}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}G\mathbf{q}-\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\\ -\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}G\mathbf{q}+\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\\ \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}G\mathbf{q}-\mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\end{array}\right)\end{array}$$ となる. ここで, 列ベクトル$\mathbf{x},\mathbf{y}$に対して, $$\begin{array}{rl}{\left(\mathrm{diag}(\mathbf{x})M\mathbf{y}\right)}_i& =\sum _j\left(\sum _k{\delta }_{ik}{\mathbf{x}}_k{M}_{kj}\right){\mathbf{y}}_j\\ & =\sum _j{M}_{ij}{\mathbf{x}}_i{\mathbf{y}}_j\\ & =\sum _j{\left(\mathrm{diag}(\mathbf{x})M\mathrm{diag}(\mathbf{y})\right)}_{ij}\end{array}$$ である. したがって, $$\begin{array}{r}\left(\begin{array}{c}-\mathrm{i}\mathrm{diag}(\mathbf{q}){}^{†}G{}^{†}\mathbf{z}\\ \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}\mathbf{z}\\ \mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}\mathbf{z}\end{array}\right)=-\mathrm{i}\left(\begin{array}{c}\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}G\mathbf{q}-\mathrm{diag}(\overline{\mathbf{q}})G{}^{†}P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\\ -\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}G\mathbf{q}+\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }})P{}^{\mathsf{T}}P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\\ \mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}G\mathbf{q}-\mathrm{i}\mathrm{diag}({\mathrm{e}}^{-\mathrm{i}\mathbf{\varphi }})G{}^{†}P{\mathrm{e}}^{\mathrm{i}\mathbf{\psi }}\end{array}\right)\end{array}$$ は, $-\mathrm{i}(B{}^{†}B)\circ TT{}^{†}$の各行を最初の$N+M$列分だけ足したものに等しい. したがって, $$\begin{array}{r}J({\mathbf{\theta }}_t){}^{\mathsf{T}}f(\mathbf{\theta })=\mathrm{sumcol}(\mathcal{I}[(B{}^{†}B)\circ TT{}^{†}],N+M).\end{array}$$ ここで, $\mathrm{sumcol}(X,n)$は, $X$の各行を最初の$n$列分だけ足し合わせて得られる列ベクトルである.

なお実際にこのLM法を実装してみるとわかるが, 実は $$\begin{array}{rl}\mathbf{a}& =[a_0,\dots ,a_{N-1}]{}^{\mathsf{T}}=[a,\dots ,a]{}^{\mathsf{T}}\end{array}$$ という風に, 振動子の出力を(一律に)固定してしまうほうが, 再現性の意味でも, 計算量の意味でも性能がいい⁶. このときは, 例えば, パラメータは$\mathbf{\theta }=[\mathbf{\varphi }{}^{\mathsf{T}},\mathbf{\psi }{}^{\mathsf{T}}]{}^{\mathsf{T}}$として, $$\begin{array}{rl}B& =\left(\begin{array}{cc}G& -P\end{array}\right),\\ T& =\left(\begin{array}{c}\overline{\mathbf{q}}\\ {\mathrm{e}}^{-\mathrm{i}\mathbf{\psi }}\end{array}\right)\end{array}$$ を計算すれば良い.

勾配降下法とGauss-Newton法

なお, 更新式を $$\begin{array}{r}{\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-(J({\mathbf{\theta }}_t){}^{\mathsf{T}}J({\mathbf{\theta }}_t){)}^{-1}J({\mathbf{\theta }}_t){}^{\mathsf{T}}f(\mathbf{\theta })\end{array}$$ としたものはGauss-Newton法と呼ばれ, $$\begin{array}{r}{\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-\lambda J({\mathbf{\theta }}_t){}^{\mathsf{T}}f(\mathbf{\theta })\end{array}$$ としたものは勾配降下法 (Gradient-Descent), あるいは, 最急降下法 (Steepest-Descent) と呼ばれる. LM法はこれらの間の子である.