離散化に由来する不確実性定量化の近年の動向

宮武勇登 (大阪大学)

常微分方程式の初期値問題を陰的な数値解法で数値計算する場合

2009年頃（学部4年生）の個人的な疑問

誤差は時間発展に伴い（通常）増加する傾向

（実際の計算の誤差を知るのは困難だが...）

(非)線形ソルバの収束の閾値、

時間発展に伴って変更できると効率的？

しかし

当時は構造保存数値解法に強い関心

踏み込んで考察することはなかった

（閾値は十分小さくしないと構造保存性が失われる）

数値計算の不確実性定量化って何？

誤差を「定量的」に評価

さらに、「確率95%で誤差はこれくらい」というようなことがいえればベター

既存の「決定論的」な数値解法を確率統計（特にベイズ統計）の視点で理解

（精度保証とは似て非なる...イメージ）

様々なことに役立てたい

計算の効率化（無駄の削減）
逆問題に対してより適切な不確実性定量化

「これよりは良い！悪い！」というより、「だいたいこれくらい！」

UQ: Uncertainty Quantification

精度がすごく良いわけでもすごく悪いわけでもない計算が主なターゲット

なぜ、数値計算のUQが重要なのか？

「超高精度な計算は必要ないが，ラフ過ぎてもダメ」というシチュエーションは案外多い

例：PDEの計算，画像処理，逆問題，...

実際の計算の精度について，一定の信頼度を保って定量的評価ができれば有用

例えば逆問題の場合...

通常の理論で導いた95%信頼区間には，数値計算の影響は加味されていない
数値計算の誤差の定量化ができれば，バイアスの評価に繋がり，
より適切な信頼区間が構築できる

結局のところ数値解析との関係は？

（従来型の）数値解析

数値計算のUQ

計算前に定性的評価

計算後に定量的評価

関数解析など

確率・統計

いつ評価？

主要な道具

【私見】　数値解析の一トピックとみなしてもよい気がします

常微分方程式

今日、お話すること、しないこと

確率的手法

その他のトピック

数値線形代数（CG法など）
数値積分
連続最適化
偏微分方程式

Probabilistic Numerics

A free electronic version
for personal use only is available

昨年出版！

確率的でない（陽には用いない）手法

ODEフィルタ・スムーザー
摂動解法

データ活用手法（逆問題）
parareal法を活用した手法

ベイズ的再解釈

様々な動機で様々な手法が開発されていて

それぞれに利点と（大きな）欠点がある

計算の過程で確率を扱う

あとで比較したりはしますが，基本的に独立した内容です

特に微分方程式の数値計算に関連するUQの研究は...

計算のUQ研究の動向

数値解析学者はほぼ不在
欧米でも統計・データ同化・機械学習などで著名な研究者も参入しているが
コミュニティの広がりは限定的（結構仲良くしてくれる）

日本では...

「計算のUQ」と謳っている研究はほとんどないように思われる
（実は融合的な研究がありうるのではと思っている）精度保証分野も強力

確率的手法

確率的でない（陽には用いない）手法

ODEフィルタ・スムーザー
摂動解法

データ活用手法（逆問題）

ベイズ的再解釈

parareal法を活用した手法

ODEフィルタ（とスムーザー）

（数式が煩雑になるのでアルゴリズムの詳細は示しません）

何らかの方法で，1ステップを計算（予測分布を出力）
何らかの方法で，予測分布を修正

基本的な考え方

どのように実現するか

カルマンフィルタ（およびその拡張版）の理論・手法を応用

例えば，各時刻で修正された分布（フィルタ分布）の

平均：数値解

分散：誤差の定量化

とみなす

カルマンフィルタ

状態空間モデル

\bm{x}_{j} = F_j \bm{x}_{j-1} + \bm{\xi}_j % y_{j+1} &= H v_{j+1} + \eta_{j+1}

\bm{y}_{j} = H _j\bm{x}_{j} + \bm{\eta}_{j}

一期先予測：

フィルタ：

（非線形モデルの場合に，線形化した行列を用いて予測・フィルタを行う手法を「拡張カルマンフィルタ」という）

\(p(\bm{x}_j\mid \bm{y}_{1:j-1})\)

\(p(\bm{x}_j\mid \bm{y}_{1:j})\)

\(\bm{x}_0 \sim \mathrm{N} (\bm{x}_{0\mid 0},V_{0\mid 0})\) : given に対して

\(\bm{x}_j \mid \bm{y}_{1:j} \sim \mathrm{N} (\bm{x}_{j\mid j},V_{j\mid j})\)

\(\bm{x}_j \mid \bm{y}_{1:j-1} \sim \mathrm{N} (\bm{x}_{j\mid j-1},V_{j\mid j-1})\)

を交互に計算

\(\bm{x}_{j\mid j-1} = F_j \bm{x}_{j-1\mid j-1}\)

\(V_{j\mid j-1} = F_j V_{j-1\mid j-1} F_j^\top + G_j Q_j G_j^\top\)

\( K_t = V_{j\mid j-1} H_j^\top (H_j V_{j\mid j-1} H_j^\top + R_j)^{-1}\)

\(\bm{x}_{j\mid j} = \bm{x}_{j\mid j-1}＋K_j (\bm{y}_j - H_j \bm{x}_{j\mid j-1})\)

\(V_{j\mid j} = V_{j\mid j-1}- K_j H_j V_{j\mid j-1}\)

システムモデル

観測モデル

\(\sim \mathrm N (\bm{0},Q_j)\)

\(\sim \mathrm N (\bm{0},R_j)\)

システムモデル

ODEに対する状態空間モデル

観測モデル

\mathbf{X}(t) : \begin{bmatrix} \bm{x}(t) \\ \bm{x}^\prime(t)\\ \vdots\\ \bm{x}^{(\nu-1)}(t) \end{bmatrix}

SDE: \( \mathrm d \mathbf{X}(t)=F\mathbf{X}(t) \mathrm dt +L \mathrm d \mathbf{\omega}_t\) を離散化して

\mathbf{x}_{j} = A(h) \mathbf{x}_{j-1} + \bm{\xi}_j % y_{j+1} &= H v_{j+1} + \eta_{j+1}

F = \begin{bmatrix} O & I & O & \cdots & O \\ O & O & I & O & \vdots \\ \vdots & & \ddots & \ddots & O \\ O & & & O & I \\ O & O & \cdots & O & O \end{bmatrix}, \quad L = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \vdots \\ \mathbf{0}\\ \sigma \mathbf{1} \end{bmatrix}

A(h) = \exp(hF), \quad \bm{\xi}_j \sim \mathrm N (\mathbf{0},Q(h)), % & Q(h) = \int_0^h \mathrm{e}^{Fs} LL^\top \mathrm{e}^{F^\top s} \,\mathrm ds

\mathbf{y}_{j} = E_1 \mathbf{x}_j - \bm{f}(E_0 \mathbf{x}_j) + \bm{\eta}_j % y_{j+1} &= H v_{j+1} + \eta_{j+1}

E_0 = \begin{bmatrix} I & O & O & \cdots \end{bmatrix}, \quad E_1 = \begin{bmatrix} O & I & O & \cdots \end{bmatrix}

観測は常に \(\mathbf{y}_j = \bm{0}\) であるとする

拡張カルマンフィルタにおける線形化の例

\( H = E_1\), \( H = E_1 - \nabla \bm{f} (E_0 \mathbf{x}_{j\mid j-1}) E_0\)

\( \displaystyle \frac{\mathrm d\bm{x}}{\mathrm dt} = \bm{f}(\bm{x}), \quad \bm{x}(0) = \bm{x}_0 \in \mathbb{R}^d\)

Q(h) = \int_0^h \mathrm{e}^{Fs} LL^\top \mathrm{e}^{F^\top s} \,\mathrm ds

相当の \(d\nu\) 次元ベクトル

\(\bm{x}^\prime \) 相当

\(\bm{x} \) 相当

ODEフィルタのポイント

高階微分相当もまとめた状態ベクトル
- 初期時刻における高階微分の計算は
  「Taylorモード自動微分」
ベクトル場 \(\bm{f}\) の情報は使わなくてもよい
かなり自由度がある

観測は常に \( \bm{0}\) ：

システムモデル

観測モデル

SDE: \( \mathrm d \mathbf{X}(t)=F\mathbf{X}(t) \mathrm dt +L \mathrm d \mathbf{\omega}_t\) を離散化して

\mathbf{x}_{j} = A(h) \mathbf{x}_{j-1} + \bm{\xi}_j % y_{j+1} &= H v_{j+1} + \eta_{j+1}

A(h) = \exp(hF), \quad \bm{\xi}_j \sim \mathrm N (\mathbf{0},Q(h)), % & Q(h) = \int_0^h \mathrm{e}^{Fs} LL^\top \mathrm{e}^{F^\top s} \,\mathrm ds

\mathbf{y}_{j} = E_1 \mathbf{x}_j - \bm{f}(E_0 \mathbf{x}_j) + \bm{\eta}_j % y_{j+1} &= H v_{j+1} + \eta_{j+1}

E_0 = \begin{bmatrix} I & O & O & \cdots \end{bmatrix}, \quad E_1 = \begin{bmatrix} O & I & O & \cdots \end{bmatrix}

観測は常に \(\mathbf{y}_j = \bm{0}\) であるとする

Q(h) = \int_0^h \mathrm{e}^{Fs} LL^\top \mathrm{e}^{F^\top s} \,\mathrm ds

\(\bm{x}^\prime \) 相当

\(\bm{x} \) 相当

残差の計算により \(\bm{f}\) の情報を取り込む

FitzHugh—Nagumo model

数値実験（「驚異的」にうまく動作するわけではない）

\displaystyle \begin{cases} \frac{\mathrm dV}{\mathrm dt} =c \big(V - \frac{V^3}{3} + R \big) \\ \frac{\mathrm dR}{\mathrm dt} = -\frac{1}{c} \big( V-a+bR\big) \end{cases}

EK1(\(\nu=1\))

EK1(\(\nu=2\))

EK1(\(\nu=3\))

ProbNumDiffEq.jl
黒線：厳密解（reference solution）
\(\mathrm{abstol}=10^{-2}, \mathrm{reltol}=10^{-1}\) として刻み幅制御を適用済

\(V\)

\(R\)

EK0： \( H = E_1\),
EK1： \( H = E_1 - \nabla \bm{f} (E_0 \mathbf{x}_{j\mid j-1}) E_0\)

Python や Julia でパッケージ化はされているが，まだまだ実用レベルではない
- 近似解についての収束解析は議論されている
- 一方で，UQについての理論はほとんど無く，
  「試してみて時々（パラメータによっては）うまく動作する」という域をでていない

感想？？

粒子フィルタとの組み合わせなどによる性能の向上は報告されている

Nordsieck法（古典的解法の一種）との関連が指摘されている

確率的手法

確率的でない（陽には用いない）手法

ODEフィルタ・スムーザー
摂動解法

データ活用手法（逆問題）

ベイズ的再解釈

parareal法を活用した手法

決定論的解法 + 各時間ステップごとに摂動を加える (Conrad et al., 2016)

摂動解法

\( \displaystyle \frac{\mathrm dx}{\mathrm dt} = f(x), \quad x(0) = x_0 \in \mathbb{R}^d\)

決定論的解法の 刻み幅 をランダムに選ぶ (Abdulle, Garegnani, 2020)

\( \Psi_h:\mathbb{R}^d\to \mathbb{R}^d\): Runge–Kutta法などの通常の一段解法として

\widehat{X}_n = \Psi_{h_n} ( \widehat{X}_{n-1}) + \xi_n (h_n), \quad \widehat{X}_0 = x_0

摂動

\widehat{X}_n = \Psi_{h_n} ( \widehat{X}_{n-1}) + \xi_n (h_n), \quad \widehat{X}_0 = x_0

\widehat{X}_n = \Psi_{H_n} ( \widehat{X}_{n-1}) , \quad \widehat{X}_0 = x_0

randomized step size

摂動を加えれば，\(\widehat{X}_n\) は確率変数

摂動の「心」

\widehat{X}_n = \Psi_{h_n} ( \widehat{X}_{n-1}) + \xi_n (h_n), \quad \widehat{X}_0 = x_0

局所誤差を確率変数でモデル化

各 \(t_n\) において分布が生成される

（実際の計算では複数のサンプルによるヒストグラム）

摂動解法の性質（収束解析）

\widehat{X}_n = \Psi_{h_n} ( \widehat{X}_{n-1}) + \xi_n (h_n), \quad \widehat{X}_0 = x_0

仮定

\displaystyle \mathbb{E} \bigg( \max_{n=0,\dots,N} \| \widehat{X}_n - x(t_n) \| \bigg) \leq C h^{\min(p,q)}

定理

(Lie, Stuart, Sullivan 2019; Conrad et al., 2017)

\( \Psi_h\): \(q\) 次精度解法
\( \exist C\), \( \displaystyle \mathbb{E} \big( \| \xi_n(t) \xi_n(t)^\top\|_{\text{F}}\big)\leq C t ^{2p+1}\)
\( f\) についての様々な仮定

仮定

局所誤差 \(\approx \mathrm{O}(h^{q+1})\)

標準偏差 \(\approx \mathrm{O}(h^{p+1/2})\)

\(p \geq q\): ランダムな出力は，決定論的ソルバより悪くはない
\(p < q\): ランダムな出力の大域誤差の期待値は，決定論的ソルバより悪い
Recommend: \(p=q\)

複数サンプルによるヒストグラム

Lorenz system

\displaystyle \begin{cases} \frac{\mathrm dx}{\mathrm dt} = -px + py \\ \frac{\mathrm dy}{\mathrm dt} = -xz+rx-y \\ \frac{\mathrm dz}{\mathrm dt} = xy-bz \end{cases}

\(\Psi_h\): 通常のRunge–Kutta法
\( \xi(h)\): 正規分布
\(q=p=4\)

Lorenz system

各時刻で「ヒストグラム」の形状から
誤差の大きさや計算の信頼度を
定量化できそうな気分になる
ただし，絶対的な根拠があるわけではない

\(x\)

摂動によって，計算の誤差が小さくなることはない（これはたぶん仕方ない）

不満な点

たくさんサンプルを取る必要があるため，計算コストは著しく増大

適切な摂動がよくわからない

\(\Psi_h\) が構造保存解法の場合，摂動によって構造保存性が破壊される

たとえば，正規分布の場合，\(\mathrm{N}(\mathbb{0},\Sigma h^{2q+1})\) とするれば良さそうだが，

適切な \(\Sigma\) はどう選択するのか？

これなら，\(\Psi_h\) が構造保存的な場合，ランダムな出力も構造保存的になる
ただし，正値をとる分布が望ましい

刻み幅をランダムに！

\widehat{X}_n = \Psi_{H_n} ( \widehat{X}_{n-1}) , \quad \widehat{X}_0 = x_0

randomized step size

(Abdulle, Garegnani, 2020)

確率的手法

確率的でない（陽には用いない）手法

ODEフィルタ・スムーザー
摂動解法

データ活用手法（逆問題）

ベイズ的再解釈

parareal法を活用した手法

背景の（超簡略版）問題設定

逆問題：摂動の代わりにデータを使う

\displaystyle \frac{\mathrm d}{\mathrm dt} x(t) = f(x(t)), \quad x(0) = \theta

モデル（微分方程式）：given (モデルの不確実性はないと仮定)
幾つかのシステムパラメータや初期値：unknown

\displaystyle y_1,\dots,y_N, \quad (y_n = Hx(t_n;\theta^*)+\varepsilon_n)

ノイズ付き時系列観測データ：known (観測オペレータ\(H\)も既知)

最尤推定の定式化

\displaystyle \min_{\theta\in\Theta} \sum_{n=1}^N \| Hx(t_n;\theta)-y_n\|^2

\displaystyle \min_{\theta\in\Theta} \sum_{n=1}^N \| H\widehat{x}(t_n;\theta)-y_n\|^2

Ideal

Reality (use approx. sol.)

観測オペレータ

観測 at \(t=t_n\)

※観測オペレータの推定

　 Akita, M., Furihata, 2022

数値解を用いる影響

\displaystyle \min_{\theta\in\Theta} \sum_{n=1}^N \| Hx(t_n;\theta)-y_n\|^2

\displaystyle \min_{\theta\in\Theta} \sum_{n=1}^N \| H\widehat{x}(t_n;\theta)-y_n\|^2

Ideal

Reality (use approx. sol.)

\(\hat{\theta}_{\mathrm{ML}}\)

\(\hat{\theta}_{\mathrm{QML}}\)

バイアス： \( \mathrm{E}_{\theta} [\hat{\theta}_{\mathrm{QML}}] - \theta \neq 0\)

最小二乗誤差（mean squared error）： \( \mathrm{E}_{\theta} \big[\|\hat{\theta}_{\mathrm{ML}}-\theta\|^2\big] \) vs \( \mathrm{E}_{\theta} \big[\|\hat{\theta}_{\mathrm{QML}}-\theta\|^2\big] \)

MLE

Runge–Kutta

中点則

【考察】

数値計算の精度によって，

推定の性能は大きく異りうる

仮定：真のパラメータ \(\theta^*\) は既知

離散化誤差（各時刻の誤差）のモデル

実際には，例えば \(\theta\) の推定と（以下の）離散化誤差の定量化を交互に繰り返す

※ 複数の変数を同時に扱う場合は \( \mathrm{N} (0,\Sigma_n)\)

\( r_n^2 := (y_n - \widehat{x}_n)^2 \sim (\gamma^2 + \sigma_n^2)\chi_1^2\)

残差 \(=\) \((\) データ - 数値解 \()^2\) \(\sim\) \(\chi^2\)-分布（自由度1）

(Matsuda, M., SIAM/ASA UQ, 2021)

観測と数値解から「各時刻での誤差」を推定

以下の目的：

誤差を「確率変数」でモデル化

\( \xi_n = \hat{x}_n - x(t_n) \sim \mathrm{N} (0,\sigma_n^2)\)

【気持ち】 \(\sigma_n\) が誤差のスケールを反映してほしい

観測ノイズがガウシアン \( y_n \sim \mathrm{N}(x(t_n),\gamma^2)\) ならば（かつ，簡単のため観測オペレータ \(H=\mathrm{id.}\) ならば）

データと数値解を比べるモデル →　\(\sigma_1^2,\dots,\sigma_N^2\) を推定

観測データと数値解から \(\sigma_n^2\) を推定

制約付き尤度最大化

\displaystyle \max_{\sigma_1^2,\dots,\sigma_N^2} \underbrace{\log p (r_1^2,\dots,r_N^2 \mid \sigma_1^2,\dots,\sigma_N^2)} _{\displaystyle =\sum_{n=1}^N \log p(r_n^2\mid \sigma_n^2)} \quad \text{s.t.} \quad \sigma_1^2 \leq \cdots \leq \sigma_N^2

凸最適化（最適解は一意）
単調回帰理論（および手法）に帰着
PAVA (pool adjacent violators algorithm) で高速に最適解が求まる
- 計算量：最悪でも \(\mathrm{O}(N^2)\)，通常は \(\mathrm{O}(N)\)
- 微分方程式の一回の時間発展計算に比べて無視できる

FitzHugh–Nagumoモデル

数値実験

\displaystyle \begin{cases} \frac{\mathrm dV}{\mathrm dt} =c \big(V - \frac{V^3}{3} + R \big) \\ \frac{\mathrm dR}{\mathrm dt} = -\frac{1}{c} \big( V-a+bR\big) \end{cases}

\(V\) の誤差について

誤差のだいだいの振る舞いは捉えているようにみえる
当然，「減る振る舞い」は捉えられない（捉えられる方が良いかどうかは問題依存）

正則化によって，単調制約を弱めることも可能

\displaystyle \min_{\sigma_1^2,\dots,\sigma_N^2} - \log p (r_1^2,\dots,r_N^2 \mid \sigma_1^2,\dots,\sigma_N^2) \quad \text{s.t.} \quad \sigma_1^2 \leq \cdots \leq \sigma_N^2

\displaystyle \min_{\sigma_1^2,\dots,\sigma_N^2} - \log p (r_1^2,\dots,r_N^2 \mid \sigma_1^2,\dots,\sigma_N^2) +\lambda \sum_{n=1}^N \max (\eta_i - \eta_{i+1},0)

\displaystyle \eta>0, \quad \eta_n := -\frac{1}{\gamma^2+\sigma_n^2}

正則化

(Matsuda, M., arXiv, 2021)

微分方程式の計算の不確実性定量化は今後重要な研究トピックになると思う
過去10年余り，様々なアイデアが提唱されているが，どれも決定打に欠ける
精度保証分野も含め，数値解析は大きく貢献できるはず！

まとめ

【私見】

【何が難しいのか】

時間発展に伴い誤差が増減する振る舞いを定量化しないといけない
→ さらに，誤差の定量化の信頼性をどのように議論すべきかは非自明
計算コストをかければ，それなりによい定量化手法は作れる
（そもそもあまり計算コストをかけたくないから，ほどほどの計算を行い，
その保証をつけたい，というのが動機の一つ）