常微分方程式の数値解析とデータサイエンス
宮武 勇登 （大阪大学）
IBIS 2022 チュートリアル
何の図がご存知ですか？？
- Newton の Principia より引用
- Keplerの第2法則（角運動量保存）の説明
- 誤差逆伝播法 と密接な関係
（と宮武は感じている...）
自己紹介
略歴
2018 -
大阪大学サイバーメディアセンター
准教授
2015 - 2017
名古屋大学工学研究科 助教
2015
日本応用数理学会・日本数学会・SIAM
（理学部数学科・情報科学研究科情報基礎数学専攻）
所属学会
- 微分方程式の 数値解析
- 数値線形代数
- 数値計算の不確実性定量化
- 地震分野などとの協働
応用数学の一分野としての 数値解析 とは
問題を数値的に解析するための
- アルゴリズムの開発，および
- その数学解析
連立一次方程式や行列関数など
行列に関連する計算アルゴリズムの研究
応用数学の一分野としての 数値解析 とは
博士 （情報理工学） 東京大学
微分方程式
時空間に 連続的 に変化する現象のモデル方程式
ニュートンの運動方程式
ナビエ・ストークス方程式
\( \displaystyle m\frac{\mathrm d^2 \bm{r}}{\mathrm d t^2} =F\)
\( \displaystyle \frac{\partial \bm{v}}{\partial t} + (\bm{v}\cdot \nabla)\bm{v} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \bm{v} + \bm{g}\)
時間連続な現象のモデル
時空間連続な現象のモデル
常微分方程式
ODE (ordinary differential equation)
偏微分方程式
PDE (partial differential equation)
通常は 物理的な現象 や 経済 のモデル
微分方程式 と 情報学
行列
ODE Net
(Chen et al. 2018)
動的低ランク近似
最適化
画像処理
?
今日お話したいこと
- 入門書に書いてある数値計算テクニックでは多くの場合不十分
- 4つのケースの話
- 入門書に書いてある数値計算手法
- 入門書には書いていない数値計算手法（の考え方）
ODENet + 動的低ランク近似 + 最適化 （ + 画像処理）
スライド１枚で話題提供
近年、各種ライブラリは充実してきているけれど、
それはそれで解法の選択が困難に...
数値解析の入門書に書いてある内容
常微分方程式の数値計算とは
基本的には 離散変数法
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \bm{f} (\bm{u}(t)), \quad \bm{u}(0)=\bm{u}_0 \in \mathbb{R}^d\)
離散時刻 \(t=t_1,t_2,\dots\) 上の近似解を求める
（刻み幅が等間隔ならば \(t_n=nh\)）
\(\bm{u}_1\approx \bm{u}(t_1)\)
\(\bm{u}_2\approx \bm{u}(t_2)\)
\(\vdots\)
time
刻み幅：\(h\) or \(\Delta t\)
常微分方程式の数値解法の分類
離散変数法の代表的解法
一段解法
多値法
\(\bm{u}_n\) から \(\bm{u}_{n+1}\) を計算
\(\bm{u}_{n-k},\dots,\bm{u}_n\) から \(\bm{u}_{n+1}\) を計算
- Runge–Kutta法
- 分離型Runge–Kutta法
- splitting法 など
- 線形多段解法
- 一般線形法 など
- Runge–Kutta法
注意：公式を表すとき
\(\bm{u}_n\mapsto\bm{u}_{n+1}\) で表現したり
\(\bm{u}_0\mapsto \bm{u}_1\) で表現したりします
Euler法（オイラー法）
\( \displaystyle \frac{\bm{u}_{n+1}-\bm{u}_n}{h} = \bm{f}(\bm{u}_n)\)
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \bm{f} (\bm{u}(t)), \quad \bm{u}(0)=\bm{u}_0 \in \mathbb{R}^d\)
\( \displaystyle \lim_{h\to 0}\frac{\bm{u}(t_n+h)-\bm{u}(t_n)}{h} = \bm{f}(\bm{u}(t_n))\)
\( \displaystyle \frac{\bm{u}(t_n+h)-\bm{u}(t_n)}{h} \approx \bm{f}(\bm{u}(t_n))\)
微分を近似
「微分を近似する」という考え方はシンプルだけど、
この視点でできることは限られる
積分の近似
\(0\)
\(h\)
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \bm{f} (\bm{u}(t)), \quad \bm{u}(0)=\bm{u}_0 \in \mathbb{R}^d\)
\( \displaystyle \bm{u}(t_n+h) = \bm{u}(t_n) + \int_0^h \bm{f}(\bm{u}(t_n + s))\,\mathrm ds\)
\( \displaystyle \approx \bm{u}(t_n) + h \bm{f}(\bm{u}(t_n ))\)
\( \displaystyle \bm{u}_{n+1}= \bm{u}_n + h \bm{f}(\bm{u}_n)\)
\( \bm{f}(\bm{u}(t_n +s)) \)
\(s\)
厳密解の積分表現
積分を近似
右辺ベクトルの近似
\(\bm{u}(t)\)
\(\bm{f}(\bm{u}_0)\)
\(\bm{u}(0)=\bm{u}_0\)
Euler法の数値解
傾きの平均をとってみる
\(\bm{u}(t)\)
\(\bm{k}_1=\bm{f}(\bm{u}_0)\)
\(\bm{u}(0)=\bm{u}_0\)
Heun法の数値解
\(\bm{k}_2=\bm{f}(\bm{u}_{\text{Euler}})\)
\(\displaystyle \frac{\bm{k}_1+\bm{k}_2}{2}\)
（注意） Heun（ホイン）法は Heun ではなく Runge が発見した公式 (1895)
（Improved Euler法と呼ばれることも）
（入門書に載っている）Runge–Kutta法
4段4次の陽的RK法
\(\bm{k}_1 =\bm{f}(\bm{u}_0) \)
\(\bm{k}_2 =\bm{f}(\bm{u}_0+\frac{h}{2}\bm{k}_1) \)
\(\bm{k}_3 =\bm{f}(\bm{u}_0+\frac{h}{2}\bm{k}_2) \)
\(\bm{k}_4 =\bm{f}(\bm{u}_0+h\bm{k}_3) \)
\(\bm{u}_1 = \bm{u}_0 + h \big( \frac{1}{6}\bm{k}_1 + \frac{2}{6}\bm{k}_2 + \frac{2}{6}\bm{k}_3 + \frac{1}{6}\bm{k}_4 \big) \)
\begin{array}{c|cccc}
0 & & & &\\
1/2 & 1/2 & & & \\
1/2 & & 1/2 & & \\
1 & & & 1 & \\
\hline
& 1/6 & 2/6 & 2/6 & 1/6
\end{array}
Butcher配列
（こんなふうに表現することも）
- Kuttaは右の公式も発見している
- 実はこの公式の方が誤差は小さいが、
上の公式の方が覚えやすいため広く利用されている
（注意）
\begin{array}{c|cccc}
0 & & & &\\
1/3 & 1/3 & & & \\
2/3 & -1/3 & 1 & & \\
1 & 1 & -1 & 1 & \\
\hline
& 1/8 & 3/8 & 3/8 & 1/8
\end{array}
発見したのはKutta (1901)
Runge–Kutta法
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \bm{f} (\bm{u}(t),t), \quad \bm{u}(0)=\bm{u}_0 \in \mathbb{R}^d\) に対する一般の \(s\) 段 RK法
\(\displaystyle \bm{k}_i =\bm{f}\big(\bm{u}_n+h\sum_{j=1}^s a_{ij} \bm{k}_j,t_n + c_j h\big) \) (\(i=1,\dots,s\))
\(\displaystyle \bm{u}_{n+1} = \bm{u}_n + h \sum_{i=1}^s b_i \bm{k}_i \)
Butcher配列
\begin{array}{c|cccc}
c_1 & a_{11} & \cdots & a_{1s} \\
\vdots & \vdots & \ddots & \vdots \\
c_s & a_{s1} & \cdots & a_{ss} \\
\hline
& b_1 &\cdots & b_s
\end{array}
\begin{array}{c|c}
\bm{c} & A \\
\hline \\[-10pt]
& \bm{b}^\top
\end{array}
=
\(A\) が「下三角行列（対角は 0 ）」ならば陽解法 → 陽的RK法
有名な公式
陽的Euler法
\begin{array}{c|c}
0 & 0 \\
\hline
& 1
\end{array}
\begin{array}{c|cc}
0 & & \\
1/2 & 1/2 & \\
\hline
& 0 & 1
\end{array}
\begin{array}{c|cc}
0 & & \\
1 & 1 & \\
\hline
& 1/2 & 1/2
\end{array}
\begin{array}{c|c}
1 & 1 \\
\hline
& 1
\end{array}
\begin{array}{c|cc}
\displaystyle \frac12 - \frac{\sqrt{3}}{6} &\displaystyle \frac{1}{4} & \displaystyle\frac{1}{4} - \frac{\sqrt{3}}{6}\\[5pt]
\displaystyle\frac12 + \frac{\sqrt{3}}{6} & \displaystyle\frac{1}{4} +\frac{\sqrt{3}}{6} & \displaystyle\frac14 \\[5pt]
\hline
& 1/2 & 1/2
\end{array}
\begin{array}{c|c}
1/2 & 1/2 \\
\hline
& 1/2
\end{array}
Rungeの公式
Rungeの公式
陰的Euler法
中点則
Gauss公式
陽解法
陰解法
数値解法の精度
常微分方程式の数値解法の精度
以下のように評価できる数値解法を 「\({\color{magenta}p}\) 次精度解法」 と呼ぶ
\(\|\bm{u}(h)-\bm{u}_1\|= \mathrm{O}(h^{{\color{magenta}p}+1})\)
これは 1ステップ毎 に混入する「局所誤差」
大域誤差は...
とてもラフに議論すれば、目標時刻を \(T\) としたとき、
\(t_N = T(=Nh)\) における大域誤差は
\(Ch^{p+1}\times N = CNhh^p = CTh^p = \mathrm{O}(h^{\color{magenta}p})\)
次数が異なると...
- 次数 と 傾き がほぼ同じ
- 誤差の小さい数値解を得るためには
次数の高い解法の方が効率的なことが多い - 刻み幅を変えると誤差が今と比べて
相対的にどう変わるか予想できる - テスト問題を除いて，実際の誤差の大きさは
よく分からないことが多い
異なる刻み幅で計算したときの大域誤差
陽的RK法で \(p\) 次精度を達成する最小段数
- \(p=1\) : 陽的Euler法
- \(p=2\) : Runge (1895)
- \(p=3\) : Heun (1900)
- \(p=4,5\) : Kutta (1901)
ただし、\(5\) 次精度解放にはミスあり
Nyström (1925) で修正
\(5\) 段 \(5\) 次が可能かどうかは当時不明
→ 不可能 by Butcher (1963頃)
最小段数
\(p\)
\(1\)
\(1\)
\(2\)
\(2\)
\(3\)
\(3\)
\(4\)
\(4\)
\(5\)
\(6\)
\(6\)
\(7\)
\(7\)
\(9\)
\(8\)
\(11\)
\(\cdots\)
\(\cdots\)
\(10\)
\(17\)
- \(p=6\) : Butcher (1964)
- \(p=7\) : \(8\) 段では不可能 by Butcher (1965)
\(9\) 段公式は複数の研究者が1968年頃発見 - \(p=8\) : Curtis (1970)
\(10\) 段では不可能 by Butcher (1985) - \(p=10\) : Hairer (1978) ← ギネス記録
\(16\) 段で可能かどうかは未解決
陰的RK法は \(s\) 段 \(2s\) 次 が可能
数値解法の安定性
安定性の評価基準はたくさんあるが、最も基本的なものを紹介
安定性の評価基準はたくさんあるが、最も基本的なものを紹介
（複素数値のスカラーODE）
テスト方程式
\( \displaystyle \frac{\mathrm d}{\mathrm dt} u(t) = \lambda u(t),\)
\( \displaystyle u(0) = u_0 \in \mathbb{C}\)
\( \lambda\) は複素数で、その実部は 負 とする
解は \( \displaystyle u(t) = \mathrm{e}^{\lambda t} u_0 \) なので、\( |u(t)| = |\mathrm{e}^{\lambda t}u_0|<|u_0|\) \((t>0)\) が成り立つが、
RK法による近似解に対してもこの性質は成り立つだろうか？
Euler法の場合
\(u_1 = (1+h\lambda)u_0\) より \(|u_1|<|u_0|\) となるには
陽的Euler法
陰的Euler法
\(|1+h\lambda|<1\)
\(u_1 = u_0 + h\lambda u_1\) より \(\displaystyle u_1 = \frac{1}{1-h\lambda}u_0\)
\(h>0\) ならば必ず \(|u_1|<|u_0|\)
無条件に安定
\(h\) が大きいとこの不等式は成り立たない
Runge–Kutta法の場合
RK法の時間発展の一表現
\(u_1 = \underbrace{R(z)}_{}u_0\) (\(z=h\lambda\))
\(\displaystyle R(z) = \frac{\mathrm{det}(I-zA + z\bm{e}\bm{b}^\top)}{\mathrm{det}(I-zA)}\) \((\bm{e}=[1,\dots,1]^\top)\)
RK法の安定性因子
A安定性
- 安定性領域： \(\mathcal{R} := \{z\in\mathbb{C} \mid |R(z)|<1\}\)
- \(\mathcal{R} \supset \mathbb{C}^{-} =\{z\in\mathbb{C} \mid \mathrm{Re}(z)<0\}\) を満たすRK法を
「A安定」な解法という
（\(\mathbb{C}^{-}\)：複素平面の左半分）
陽的Runge–Kutta法の安定性領域
- \(p\) 段 \(p \) 次の陽的RK法の安定性領域は同一
（例：4段4次の陽的RK法は複数存在するが、
その安定性領域はどれも同じ） - \(p=5\) （6段） の図は一例
- どれも \(\mathbb{C}^{-}\) （=左半平面） からは程遠い
\(\lambda\) の実部が負の方に大きいと
\(h\) はものすごく小さくとらないといけない - Gauss–Runge–Kutta法 （陰的RK法） は A安定
硬い系（stiffness）
硬い系
陽的RK法で計算する場合に、
刻み幅を極めて小さくとらないと、数値的に不安定になってしまう系
RK法の場合
- 陰解法を用いる
- 安定性領域ができる限り実軸の左の方を覆うように設計された
RK法を用いる
例題
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \begin{bmatrix} -1& 1\\2 & -3\end{bmatrix} \bm{u}(t) + \begin{bmatrix}\cos t \\ 3\cos t - \sin t \end{bmatrix} ,\quad \bm{u}(0) = \begin{bmatrix} 1\\ 2\end{bmatrix}\)
例 1
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \begin{bmatrix} -2& 1\\1998 & -1999\end{bmatrix} \bm{u}(t) + \begin{bmatrix}-\cos t \\ 1999\cos t - \sin t \end{bmatrix} ,\quad \bm{u}(0) = \begin{bmatrix} 1\\ 2\end{bmatrix}\)
例 2
実はどちらも厳密解は \( \displaystyle \bm{u}(t) = \mathrm{e}^{-t} \begin{bmatrix} 1 \\ 1\end{bmatrix} + \begin{bmatrix}0 \\ \cos t\end{bmatrix}\)
微分方程式もなんとなく似通っているし、
RK法で計算した近似解も似ているはず...？
（「三井, 常微分方程式の数値解法」 の例）
陽的Runge–Kutta法 と 中点則
4段4次の陽的RK法
中点則（陰解法）
例 1
例 2
例 1 と 例 2 はほとんど同じ
（刻み幅はどちらも \(h=0.01\)）
例 2 は 2 ステップ目 で発散...
違いは何か？
固有値：\(-1, -4\)
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \underbrace{\begin{bmatrix} -1& 1\\2 & -3\end{bmatrix}}_{} \bm{u}(t) + \begin{bmatrix}\cos t \\ 3\cos t - \sin t \end{bmatrix} \)
例 1
\( \displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u}(t) = \underbrace{\begin{bmatrix} -2& 1\\1998 & -1999\end{bmatrix}}_{} \bm{u}(t) + \begin{bmatrix}-\cos t \\ 1999\cos t - \sin t \end{bmatrix} \)
例 2
固有値：\(-1, -2000\)
だいたい、\(-3 < -4h \, (h<3/4)\) くらいで安定に計算できる！
- RK法を使うなら \(h<3/2000\) くらい小さくとらないといけない
- 中点則なら特に制限なし （ただし精度の問題はある）
改めて「硬い系」とは？
- 行列 \(A\) について絶対値が大きな固有値を持つ問題は
「硬い系」になる - 硬度比が大きな問題も「硬い系」になりうる
- 硬度比が大きくても硬い系ではないかもしれないが、
初期値の変動に対しては不安定
線形系：\( \displaystyle \frac{\mathrm d}{\mathrm dt}\bm{u}(t) = A\bm{u}(t) + \bm{\phi}(t) \) に対して...
硬度比 \(\displaystyle =\)
| \(A\) の絶対値最大固有値 |
| \(A\) の絶対値最小固有値 |
例 2 は「硬い系」とまでは言えず、「硬い系に近い」という程度の問題
（注） 非線形問題の場合は、線形化してあらわれるJacobi行列の固有値に着目する
硬い系に対する数値解法
安定性領域が、（特に左方向に） 広い解法を使う！
- A安定な解法
- 「硬い系専用」の解法
- 陽解法ならば、安定性領域が実軸左方向に広がっている解法
（引用：Hairer, Wanner, Solving Ordinary Differential Equations II, 1996）
数値解析の入門書に書いてない現代的な内容
Störmer法
Newtonの運動方程式： \(\displaystyle \frac{\mathrm d^2}{\mathrm dt^2} \bm{q} = -\nabla V(\bm{q}) \) \(\Leftrightarrow\) \(\displaystyle \frac{\mathrm d}{\mathrm dt} \begin{bmatrix} \bm{q} \\ \bm{p} \end{bmatrix} = \begin{bmatrix} \bm{p} \\ -\nabla V(\bm{q}) \end{bmatrix}\)
Störmer法 (1907)
\( \displaystyle \bm{q}_{n+1/2} = \bm{q}_n + \frac{h}{2}\bm{p}_n\)
\( \displaystyle \bm{p}_{n+1} = \bm{p}_n - h \nabla V(\bm{q}_{n+1/2})\)
\( \displaystyle \bm{q}_{n+1} = \bm{q}_{n+1/2} + \frac{h}{2}\bm{p}_{n+1}\)
- Verlet法 （1967; 分子動力学） も本質的に同じ
- 数値解析学者は Störmer–Verlet法 と呼んだりする
あるいは、Newton–Störmer–Verlet法
位置 \(\bm{q}\) を陽的Euler法で半ステップ更新
運動量 \(\bm{p}\) を1ステップ更新
位置 \(\bm{q}\) を残りの半ステップ更新
Keplerの二体問題の数値計算
約1周期分の計算
約10周期分の計算
計算する時間が長くなるとStörmer法が有利
エネルギー保存則
太陽系のシミュレーション
陽的Euler法
陰的Euler法
Störmer法
Symplectic Euler法
Störmer法の
親戚（こども）解法
（岩瀬さん （阪大） 提供）
Störmer–Verlet法はなぜ優れている？
Newtonの運動方程式は Hamilton系
→ 解の時間発展写像は シンプレクティック (symplectic)
SV法による近似解の離散的な時間発展写像もシンプレクティック
→ もとの系に非常によく似たHamilton系を厳密に解いている
\(\mathbb{R}^2\) の場合：
シンプレクティック性 = 位相空間内の面積保存性
（引用：Hairer, Lubich, Wanner, Geometric Numerical Integration, 2006）
幾何学的数値積分法 ・ 構造保存数値解法
（Geometric Numerical Integration）
（Structure-preserving numerical methods）
問題 （微分方程式） の数理構造を離散化後も厳密に再現する数値解法
- シンプレクティック解法
- エネルギー保存解法 （シンプレクティック解法とは両立しない）
- エネルギー散逸解法
- 多様体上の微分方程式に対する数値解法
- 高振動系やマルチスケール問題に対する数値解法
研究者の主なタスク
- どのように構造保存アルゴリズムを設計すればよいか？
- 恩恵は何か？
Störmer法はシンプレクティックだが、さらに...
Keplerの第２法則 （角運動量保存則）
双線形形式の保存！
実は Störmer法 を含む「分離型Runge–Kutta法」という解法クラスに対して
Hamilton系に対してシンプレクティック \(\Leftrightarrow\) 任意の双線形形式を保存 \({}^\ast\)
\(\displaystyle \frac{\mathrm d}{\mathrm dt} \begin{bmatrix} \bm{q} \\ \bm{p} \end{bmatrix} = \begin{bmatrix} \bm{p} \\ -\nabla V(\bm{q}) \end{bmatrix}\), \(\displaystyle V(\bm{q}) = \frac{1}{q_1^2 + q_2^2} \)
\(q_1 p_2 - q_2 p_1=\text{const.}\)
（\({}^\ast\) 双線形形式の保存量があるとき自動的に保存性を再現）
なお、RK法に関しては
Hamilton系に対してシンプレクティック \(\Leftrightarrow\) 任意の2次の保存量を保存
ODE Net の誤差逆伝播法
問題設定
\(\displaystyle \frac{\mathrm d}{\mathrm dt} \bm{u} = \bm{f} (\bm{u},t,\bm{\theta})\), \(\bm{u}(0)=\bm{u}_0\)
Loss: \(\displaystyle \mathcal{L}(\bm{u}(T;\bm{\theta})) \)
小さくしたい
現実には... \( \mathcal{L}(\bm{u}_N (\bm{\theta}))\) を小さくしたい
\(\displaystyle \nabla_{\bm{\theta}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \) や \(\displaystyle \nabla_{\bm{\theta}}\mathcal{L}(\bm{u}_N(\bm{\theta})) \) の計算が重要！
（注意） \(\mathcal{L}\) の具体形などを限定しなければ、
データ同化 （特に4次元変分法） と本質的に同じ設定
特にこれが問題
連続の場合：アジョイント （随伴） 法
\(\dot{\bm{\lambda}} = -\big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big)^\top \bm{\lambda}\)
\(\displaystyle \dot{\bm{u}} = \bm{f} (\bm{u})\), \(\bm{u}(0)=\bm{\theta}\)
対象の方程式 ：
アジョイント方程式 ：
\(\displaystyle \nabla_{\bm{\theta}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \) の評価
\(\displaystyle \bm{\lambda}(T)=\nabla_{\bm{u}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \) を満たす解 \(\bm{\lambda}(t)\) に対して \(\displaystyle \bm{\lambda}(0)=\nabla_{\bm{\theta}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \)
簡単のため初期値を \(\bm{\theta}\) とする
アジョイント方程式は時間逆向きに解く！
上の方程式の解
変分方程式を用いた証明
\(\dot{\bm{\lambda}} = -\big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big)^\top \bm{\lambda}\)
\(\displaystyle \dot{\bm{u}} = \bm{f} (\bm{u})\), \(\bm{u}(0)=\bm{\theta}\)
対象の方程式 ：
アジョイント方程式 ：
\(\dot{\bm{\delta}} = \big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big) \bm{\delta}\)
変分方程式 ：
補題 ： \(\displaystyle \frac{\mathrm d}{\mathrm dt} \big( \bm{\lambda}^\top \bm{\delta}\big) = 0\)
双線形形式の保存！
\(\displaystyle \because \frac{\mathrm d}{\mathrm dt} \big( \bm{\lambda}^\top \bm{\delta}\big) =\dot{\bm{\lambda}}^\top \bm{\delta} + \bm{\lambda}^\top \dot{\bm{\delta}}= 0\)
特に \(\bm{\lambda}(T)^\top \bm{\delta}(T)=\bm{\lambda}(0)^\top \bm{\delta}(0)\)
連鎖律 \(\displaystyle \nabla_{\bm{\theta}}\mathcal{L}(\bm{u}(T;\bm{\theta})) = (\nabla_{\bm{\theta}}\bm{u}(T;\bm{\theta}))^\top \nabla_{\bm{u}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \) と
変分方程式の性質 \(\displaystyle \bm{\delta}(T)= (\nabla_{\bm{\theta}}\bm{u}(T;\bm{\theta})) \bm{\delta}(0) \) より
\(\displaystyle \big(\nabla_{\bm{u}}\mathcal{L}(\bm{u}(T;\bm{\theta})))^\top \bm{\delta}(T) = \big(\nabla_{\bm{\theta}}\mathcal{L}(\bm{u}(T;\bm{\theta})))^\top \bm{\delta}(0)\)
(\(\ast\))
(\(\ast\)) と比べると、\(\displaystyle \bm{\lambda}(T)=\nabla_{\bm{u}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \) を満たす解 \(\bm{\lambda}(t)\) に対して \(\displaystyle \bm{\lambda}(0)=\nabla_{\bm{\theta}}\mathcal{L}(\bm{u}(T;\bm{\theta})) \)
双線形形式の保存！
離散の場合
\(\displaystyle \dot{\bm{u}} = \bm{f} (\bm{u})\), \(\bm{u}(0)=\bm{\theta}\)
対象の方程式 ：
\(\displaystyle \nabla_{\bm{\theta}}\mathcal{L}(\bm{u}_N(\bm{\theta})) \) の評価
近似計算して \(\bm{u}_N(\bm{\theta}) \approx \bm{u}(T;\bm{\theta})\)
\(\dot{\bm{\lambda}} = -\big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big)^\top \bm{\lambda}\)
アジョイント方程式 ：
\(\dot{\bm{\delta}} = \big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big) \bm{\delta}\)
変分方程式 ：
近似計算して \(\bm{u}_N(\bm{\theta}) \approx \bm{u}(T;\bm{\theta})\)
対象の方程式と同じ手法で離散化
（実際に計算する必要はない）
\(\bm{\lambda}_n^\top \bm{\delta}_n =\text{const.}\) となるように離散化
双線形形式の保存！
終端条件 \(\displaystyle \bm{\lambda}_N=\nabla_{\bm{u}}\mathcal{L}(\bm{u}_N(\bm{\theta})) \) を課し、時間逆向きに計算すると \(\displaystyle \bm{\lambda}_0=\nabla_{\bm{\theta}}\mathcal{L}(\bm{u}_N(\bm{\theta})) \)
アジョイント方程式の離散化 （RK法の場合）
\(\dot{\bm{\lambda}} = -\big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big)^\top \bm{\lambda}\)
\(\displaystyle \dot{\bm{u}} = \bm{f} (\bm{u})\), \(\bm{u}(0)=\bm{\theta}\)
対象の方程式 ：
アジョイント方程式 ：
\(\dot{\bm{\delta}} = \big(\nabla_{\bm{u}} \bm{f}(\bm{u})\big) \bm{\delta}\)
変分方程式 ：
「対象の方程式」と「変分方程式」を \(s\) 段RK法 \((A,\bm{b},\bm{c})\) で離散化した場合
「アジョイント方程式」は \(s\) 段RK法 \((\widehat{A},\widehat{\bm{b}},\widehat{\bm{c}})\) で離散化すればよい
ただし、
\(b_i \widehat{a}_{ij} + \widehat{b}_ja_{ji} = b_i \widehat{b}_j,\)
\(i,j=1,\dots,s\)
\(\widehat{b}_i=b_i,\) \(\widehat{c}_i = c_i,\)
\(i=1,\dots,s\)
- 分離型RK法が双線形形式の保存量を自動的に再現する条件
- Newton–Störmer–Verlet法は一例
Newton–Störmer–Verlet法 再考
分離型ODE系
\begin{cases}
\dot{\bm{y}} = \bm{f}(\bm{y},\bm{z})\\
\dot{\bm{z}} = \bm{g}(\bm{y},\bm{z})
\end{cases}
に対する 分離型RK法
\(\displaystyle \bm{k}_i = \bm{f}\big(\bm{y}_0 + h\sum_{j=1}^s a_{ij}\bm{k}_j, \bm{z}_0 + h\sum_{j=1}^s \widehat{a}_{ij} \bm{l}_j\big) \)
\(\displaystyle \bm{l}_i = \bm{g}\big(\bm{y}_0 + h\sum_{j=1}^s a_{ij}\bm{k}_j, \bm{z}_0 + h\sum_{j=1}^s \widehat{a}_{ij} \bm{l}_j\big) \)
\(\displaystyle \bm{y}_1 = \bm{y}_0 + h\sum_{i=1}^s b_i \bm{k}_i\)
\(\displaystyle \bm{z}_1 = \bm{y}_0 + h\sum_{i=1}^s \widehat{b}_i \bm{l}_i\)
\begin{array}{c|c}
\bm{c} & A \\
\hline \\[-10pt]
& \bm{b}^\top
\end{array},
\begin{array}{c|c}
\widehat{\bm{c}} & \widehat{A} \\
\hline \\[-10pt]
& \widehat{\bm{b}}^\top
\end{array}
NSV法の場合
\begin{array}{c|cc}
0 & & \\
1 & 1/2 & 1/2\\
\hline
& 1/2 & 1/2
\end{array},
\begin{array}{c|cc}
1/2 & 1/2 & \\
1/2 & 1/2 & \\
\hline
& 1/2 & 1/2
\end{array}
（注意） 一般には陰解法
前頁の条件を満たす
補足 ： 2階微分も計算できます！
2階微分 = ヘッシアン \(H\)
\(H \bm{x} = \bm{b}\) のような連立一次方程式を解きたいことがほとんど
- \(H\) は通常、大規模・密行列なので、\(H\) の計算は非現実的
- CG法 （共役勾配法 ： Krylov部分空間法の一種） など を用いるなら
「ヘッシアン \(\times\) ベクトル」 の計算ができれば十分
通常のネットワークの場合
動的低ランク近似
行列の低ランク近似
\(A\in\mathbb{R}^{m\times n}\) に対し、ランクが \(k\) 以下の \(m\times n\) 型行列の中で
\(\displaystyle \min_{\text{rank}(X)\leq k}\| X - A\|_{\mathrm{F}}^2\)
フロベニウスノルムの意味で \(A\) との距離が最小になるものを求めよ
（注） \(\displaystyle \|A\|_{\mathrm{F}}^2 := \sum_{i=1}^m \sum_{j=1}^n a_{ij}^2\)
特異値分解で解ける！
特異値分解
ランク \( r\) の実行列 \(A\in\mathbb{R}^{m\times n}\) は
適当な \(m\) 次直交行列 \(V\) と \(n\) 次直交行列 \(U\) と用いて
\(A = U\Sigma V^\top\)
の形に分解される
ここで、\(\Sigma\) は
\Sigma =
\begin{pmatrix}
\begin{matrix}
\sigma_1 & & \\
& \ddots & \\
& & \sigma_r
\end{matrix}
& \Large{O_{r,n-r}} \\
\Large{O_{m-r,r}} &
\Large{O_{m-r,n-r}}
\end{pmatrix}
\quad
(\sigma_1 \geq \sigma_2\geq \cdots \geq \sigma_r > 0)
の形の \(m\times n\) 型行列
\(\sigma_1,\dots,\sigma_r\) : 「特異値」という
左上から大きい順に並んでいるとする
低ランク近似問題の解
\(A=\)
\(U\)
\(\Sigma\)
\(V^\top\)
\begin{matrix}
\sigma_1 & & & & \\
& \ddots & & &\\
& & \sigma_k & & \\
& & & \ddots & \\
& & & & \sigma_r
\end{matrix}
\(A_k=\)
\begin{matrix}
\sigma_1 & & \\
& \ddots & \\
& & \sigma_k
\end{matrix}
\(U_k\)
\(\Sigma_k\)
\(V_k^\top\)
最適解！
注意 ： \(U_k^\top U_k = V_k^\top V_k = I_k\) (単位行列)
Eckart–Young–Mirskyの定理
最適解 \(=A_k\)
この事実は非常に有名だが、「初等的な証明」を与えている文献は少ない
スペクトルノルム (\(l^2\)ノルム) とフロベニウスノルムの場合に限定して
初等的証明
ユニタリ不変ノルム （フロベニウスノルムも含む） に対して
初等的証明
問題設定
\(A(t)\in\mathbb{R}^{m\times n}\) : 変数 \(t\) に関して連続な行列 に対して
\(\| X(t) - A(t)\|_{\mathrm{F}}\) を最小にする ランク \(k\) 行列 \(X(t)\) を求めよ
あるいは、連続的に変化する行列の列 \(A_1,A_2,\dots\) に対して
\(\| X_i - A_i\|_{\mathrm{F}}\) を最小にする ランク \(k\) 行列 \(X_1,X_2,\dots\) を求めよ
当然、各 \(t\) や 各 \(i\) で低ランク近似したものが最適解だが
それなりに高コスト...
∵ 大規模行列の特異値分解の計算は大変...
問題を近似
\(\mathcal{M}_k\) : ランク \(k\) の \(m\times n\) 行列の集合 （多様体）
\(T_Y \mathcal{M}\) : \(Y\) における \(\mathcal{M}_k\) の接空間
\(\| X(t) - A(t)\|_{\mathrm{F}}\) を最小にする 行列 \(X(t)\in\mathcal{M}_k\) を求めよ
\(\| \dot{Y}(t) - \dot{A}(t)\|_{\mathrm{F}}\) を最小にする 行列 \(\dot{Y}(t) \in T_{Y(t)}\mathcal{M}_k\) を求めよ
- \(Y(t)\) は多様体 \(\mathcal{M}_k\) 上のある種の微分方程式の解
- \(X(0)\) を \(Y(t)\) の初期値とし、
初期値問題としての微分方程式を近似的に解く
準最適性
ある条件のもとで
\(\displaystyle \| Y(t) - X(t)\|_\mathrm{F} \leq C(t) \int_0^t \| X(s) - A(s)\|_\mathrm{F}\,\mathrm ds\)
近似解と最適解の差 \(\leq \) 定数 \(\times\) 最適解の誤差の蓄積
(Koch, Lubich, 2007)
どの程度よく近似するか？
(引用 ： Koch, Lubich, 2007)
各時刻で特異値分解したときの誤差
動的低ランク近似の解の誤差
動的低ランク近似と最適近似の差
悪くなさそう！
失敗する例
\(\displaystyle A(t) = \begin{bmatrix} \mathrm{e}^{-t} & 0 \\ 0 & \mathrm{e}^t\end{bmatrix}\) に対して
- \(\displaystyle X(t) = \begin{bmatrix} \mathrm{e}^{-t} & 0 \\ 0 & 0\end{bmatrix}\) (\(t<0\)), \(\displaystyle X(t) = \begin{bmatrix} 0 & 0 \\ 0 & \mathrm{e}^{t} \end{bmatrix}\) (\(t>0\))
- \(\displaystyle Y(t) = \begin{bmatrix} \mathrm{e}^{-t} & 0 \\ 0 & 0\end{bmatrix}\) (\(t>t_0\))
\(t_0<0\) のときに \(Y(t_0)=X(t_0)\) として
近似問題を解いた解
このとき、\(t>0\) で \(X(t)\) と \(Y(t)\) の乖離は指数的に大きくなる
\(A(t)\) の特異値が \(t=0\) で交差するのが原因
微分方程式表現
(\(X(0) = U_0 S_0 V_0^\top\))
\(\dot{S} = U^\top \dot{A} V\)
\(\dot{U} = (I-UU^\top) \dot{A} V S^{-1}\)
\(\dot{V} = (I-VV^\top) \dot{A}^\top U S^{-\top}\)
\(S(0)=S_0\)
\(U(0) = U_0\)
\(V(0) = V_0\)
\(Y(t) = U(t) S(t) V(t)^\top\) は、近似問題の解を与える
- \(U(t)\) と \(V(t)\) は自動的に \(U(t)^\top U(t) = I\), \(V(t)^\top V(t) = I\)、
\(S(t)\) はランク \(k\) 行列となる - その帰結として、\(Y(t)=U(t)S(t)V(t)^\top\) もランク \(k\)
離散化の際の要点
(\(X(0) = U_0 S_0 V_0^\top\))
\(\dot{S} = U^\top \dot{A} V\)
\(\dot{U} = (I-UU^\top) \dot{A} V S^{-1}\)
\(\dot{V} = (I-VV^\top) \dot{A}^\top U S^{-\top}\)
\(S(0)=S_0\)
\(U(0) = U_0\)
\(V(0) = V_0\)
- \(U_n\) と \(V_n\) が \(U_n^\top U_n = I\), \(V_n^\top V_n = I\) をみたすように
離散化しなければならない - いろいろな方法が提案されている
最適化
最急降下法？
\(\bm{x}_{n+1} = \bm{x}_n - h_n \nabla f(\bm{x}_n)\)
最適化
数値解析
\(\displaystyle \min_{\bm{x}\in\mathbb{R}^d} f(\bm{x})\)
に対する最急降下法
\(\displaystyle \frac{\mathrm d}{\mathrm dt} \bm{x}(t) = - \nabla f(\bm{x}(t))\)
に対する陽的Euler法
\(h_n\) の決め方
目的
最適解の発見
関数値を基準に選ぶ
誤差を基準に選ぶ
厳密解フローの再現
連続最適化と微分方程式の数値解析
最適化問題
最適化問題
例：無制約最適化
\(\displaystyle \min_{\bm{x}\in\mathbb{R}^d} f(\bm{x})\)
最適化問題
微分方程式
例：勾配流
\(\dot{\bm{x}} = -\nabla f(\bm{x})\)
最適化問題
最適化手法
例：最急降下法
数値解法
例：陽的Euler法
\(\bm{x}_{n+1} = \bm{x}_n - h_n \nabla f(\bm{x}_n)\)
\(\bm{x}_{n+1} = \bm{x}_n - h_n \nabla f(\bm{x}_n)\)
他の例： Runge–Kutta法、線形多段解法
他の例： Nesterovの加速法、Newton法
最適解の探索
離散化
連 続 極 限
最適化に適した離散化手法？
適切なダイナミクス？
安定な陽解法の活用 1/5
目的関数 \(f\) についての仮定
- \(f\) が \(L\)-平滑 ： 勾配 \(\nabla f\) が \(L\)-Lipschitz 連続
- \(\text{argmin}\ f \neq \emptyset\) ： 最適解 \( \bm{x}^*\) が存在
- \(\mu\)-強凸 ： \(\displaystyle f(\bm{x}) - \frac{\mu}{2}\|\bm{x}\|^2\) が凸関数
最急降下法の収束条件 (cf. Nesterov 2004)
\(h_n < 2/L\) のとき \( f(\bm{x}_n) \to f^*\)
安定な陽解法の活用 2/5
- 勾配流 \(\dot{\bm{x}} = -\nabla f(\bm{x})\) の線形化 ： \(\dot{\bm{\delta}}=-\nabla^2 f(\bm{x})\bm{\delta}\)
- \(f\) が \(L\)-平滑な凸関数 \(\Rightarrow\) \(-\nabla^2 f(\bm{x})\) の固有値 \(\in [-L,2]\)
\(-2 < -Lh < 0\) ならば 陽的Euler法は安定
\(\Leftrightarrow 0<h<2/L\) ： 最急降下法の収束条件と同じ！
安定性領域が実軸左方向に広がっている陽的RK法を使うと
より大きな刻み幅を使って安定に最適化問題を解けそう！
安定な陽解法の活用 3/5
（引用：Hairer, Wanner, Solving Ordinary Differential Equations II）
Text
こういう点は避けたい
左端 \(\approx 2s^2\)
安定な陽解法の活用 4, 5/5
離散勾配法の活用 1/4
ところで、刻み幅の制約を受けない陽的な解法は作れないだろうか？
A. 可能！
観察：勾配系 \(\dot{\bm{x}} = -\nabla f(\bm{x})\) に対して \(f(\bm{x}(t))\) は単調減少
\( \displaystyle \frac{\mathrm d}{\mathrm dt} f(\bm{x}) = \nabla f(\bm{x}) ^\top \dot{\bm{x}} = - \| \nabla f(\bm{x})\|^2 \leq 0\)
方程式の代入 (自明)
連鎖律 （本質）
離散勾配法の活用 2/4
連鎖律 \( \displaystyle \frac{\mathrm d}{\mathrm dt} f(\bm{x}) = \nabla f(\bm{x}) ^\top \dot{\bm{x}}\) を模倣して
- \( \displaystyle f(\bm{y})-f(\bm{x}) = \nabla_\text{d} f(\bm{x},\bm{y}) ^\top ( \bm{y} - \bm{x})\)
- \( \nabla_\text{d}f(\bm{x},\bm{x}) = \nabla f(\bm{x})\)
を満たす離散量 \( \nabla_\text{d} f(\bm{x},\bm{y})\) を用いて
勾配流 \(\dot{\bm{x}} = -\nabla f(\bm{x})\) を
\(\displaystyle \frac{\bm{x}_{n+1}-\bm{x}_n}{h} = -\nabla_\text{d} f(\bm{x}_n,\bm{x}_{n+1})\)
と離散化すると
\(\displaystyle f(\bm{x}_{n+1})-f(\bm{x}_n)=\big(\nabla_\text{d} f(\bm{x}_n,\bm{x}_{n+1})\big)^\top (\bm{x}_{n+1}-\bm{x}_n) = -h \| \nabla_\text{d}\bm{f}\|^2 \leq 0\)
\(0<h\) でありさえすれば、\(f\) は単調減少
離散勾配法
離散勾配法の活用 3/4
- 一般に離散勾配は無数に存在する
ほとんどの場合、陰解法になる
脱線？？
SOR法 （Successive Over-Relaxation法、逐次加速緩和法）
- 連立一次方程式 \( A \bm{x} = \bm{b}, \ \ A \in \mathbb{R}^{m\times m}, \ \bm{b} \in \mathbb{R}^m\) に対する代表的な定常反復法
- \(\bm{x}^{(n+1)} \mapsto \bm{x}^{(n)}\) の定義：
\(\displaystyle \bm{x}_i^{(m+1)}=(1-\omega)\bm{x}_i^{(n)} + \omega \frac{1}{a_{ii}} \big( b_i - \sum_{j=1}^{i-1}a_{ij}\bm{x}_j^{(k+1)} - \sum_{j=i}^{m}a_{ij}\bm{x}_j^{(k)} \big),\) \(i=1,\dots,m\)
- \(\omega=1\) のとき、Gauss–Seidel法； \(\omega\)：緩和パラメータ
- \(A\) が正定値対称行列のとき、
SOR法が任意の初期ベクトル \(\bm{x}^{(0)}\) に対して収束 \(\Leftrightarrow\) \(\omega\in (0,2)\)
離散勾配法の活用 4/4
- \(\dot{\bm{x}} = - D^{-1} \nabla f(\bm{x}) (= -D^{-1} (A\bm{x}-\bm{b}) ) \)
を Itoh–Abeの離散勾配 を用いて離散化すると
\(A\) が正定値対称ならば、\(0<h_n\) でありさえすれば収束する - \(\displaystyle h=\frac{2\omega}{2-\omega}\) と変換すれば、SOR法と一致！
SOR法の連続極限
画像処理
（たぶん時間切れなのでちょっとだけ）
画像処理に微分方程式？
- ODE Net を使うと自然に微分方程式が登場
- 離散勾配法を用いるなど、微分方程式を介して
最適化アルゴリズムを導出
- 画像処理の過程そのものを微分方程式でモデル化
例 ： Inpainting （画像修復）
最後に...重要そうな文献をいくつか
常微分方程式の数値解析の基礎
和書で、入門書よりは深く、
一方で専門的過ぎない
定番の教科書：1990年代半ばくらいまでの研究内容が
網羅的に整理されている
常微分方程式の数値解析の基礎
（第3版は） より最先端の内容も含んでいる
PDEも扱っている
構造保存数値解法
対象はHamilton系が中心
（やや古く、対象が限定的）
数学的にはHLWほど高度ではない
定番中の定番
数学的難易度はやや高め
一番新しい
対象はやや限定的だが読みやすい
まとめ
- 微分方程式を扱う際は、ほとんどの場合、離散化が必須
数値解析のいろいろな手法が （数値解析で忘れられているような手法も） 役立つ！
数値解析学者として嬉しい！
- これまでの数値解析学の知見では不十分なことも多い
（例）
多くの物理シミュレーションと違い、「低精度演算」の需要が大きかったりする
新しい数値解析学的研究テーマ
数値解析学者として嬉しい！
今日の内容を「すべて」含んだ本を書いています（きっと来年度中には...）
阪大での講義資料
Thank You!
IBIS2022（公開版）
By Yuto Miyatake
IBIS2022（公開版）
IBIS2022 https://ibisml.org/ibis2022/ のチュートリアルの資料です（部分的にカットしています）
