Ensemble Kalman Inversion Derivative-Free Optimization Andrew - - PowerPoint PPT Presentation
Ensemble Kalman Inversion Derivative-Free Optimization Andrew Stuart Computing and Mathematical Sciences California Institute of Technology Allen Philanthropies, Mission Control for Earth, NSF, ONR, Schmidt Futures. ISDA 2019, Kobe, Japan
Derivative-Free Optimization
Andrew Stuart
Computing and Mathematical Sciences California Institute of Technology Allen Philanthropies, Mission Control for Earth, NSF, ONR, Schmidt Futures.
ISDA 2019, Kobe, Japan
January 23rd 2019
What Is EKI? EKI For Complex Dynamics EKI With Constraints Conclusions and References
Anderson, J.L. (2001) Lorentzen, R., K. Fjelde, J. Froyen, A. Lage, G. Naevdal and E. Vefring (2001) Naevdal, G., T. Mannseth and E. Vefring (2002) Skjervheim, J. and G. Evensen (2011) Chen, Y. and D. Oliver (2002) Emerick, A. and A. Reynolds (2013) Evensen, G. (2018)
State Space Model
Dynamics Model: vn+1 = Ψ(vn) + ξn, n ∈ Z+ Data Model: yn+1 = Hvn+1 + ηn+1, n ∈ Z+ Noise Structure: ξn ∼ N(0, Σ), ηn ∼ N(0, Γ).
Sequential Optimization Perspective
Predict:
n+1 = Ψ(v (k) n ) + ξ(k) n ,
n ∈ Z+ Model/Data Compromise: J(k)
n (v) = 1
2 C
− 1
2
n+1(v −
v (k)
n+1)2 + 1
2Γ− 1
2 (y (k)
n+1 − Hv)2
Optimize: v (k)
n+1 = argminv J(k) n (v).
v (k)
n+1}.
Standard Inverse Problem Formulation
Find θ from y where G : U → Y, η is noise and y = G(θ) + η, η ∼ N(0, Γ).
State Space Estimation Formulation (Oliver et al [3], Iglesias et al [4])
Dynamics Model: θn+1 = θn, n ∈ Z+ Dynamics Model: wn+1 = G(θn), n ∈ Z+ Data Model: yn+1 = wn+1 + ηn+1, n ∈ Z+ Employ Ensemble Kalman Filter with yn+1 ≡ y.
Standard Inverse Problem Formulation
Find θ from y where G : U → Y, η is noise and y = G(θ) + η, η ∼ N(0, Γ).
State Space Estimation Formulation v = (θ, w), Ψ(v) = (θ, G(θ)) and H = (0, I):
Dynamics Model: vn+1 = Ψ(vn) + ξn, n ∈ Z+ Data Model: yn+1 = Hvn+1 + ηn+1, n ∈ Z+ Employ Ensemble Kalman Filter with yn+1 ≡ y.
EKI: Discrete Time
θ(j)
n+1 = θ(j) n + C θG n (C GG n
+ Γ)−1 y (j)
n
− G(θ(j)
n )
θ(j)(0) = θ(j)
0 .
Let Γ → h−1Γ, θ(j)
n ≈ θ(j)(nh) and h → 0:
EKI: Continuous Time Limit
˙ θ(j) = − 1 J
J
G, G(θ(j)) − y (j)
Γ θ(k),
θ(j)(0) = θ(j)
0 .
Forward Problem
Given (κ, I) ∈ L∞(D; R+) × Rm find (ν, V ) ∈ H1(D) × Rm : −∇ · (κ∇ν) = 0 ∈ D, ν + zℓκ∇ν · n = Vℓ ∈ eℓ, ℓ = 1, . . . , m, ∇ν · n = 0 ∈ ∂D\ ∪m
ℓ=1 eℓ,
∈ eℓ, ℓ = 1, . . . , m.
D ∂D el
Ohm’s Law: V = R(κ) × I.
Inverse Problem
Set κ = exp(u). Given a set of K noisy measurements of voltage V (k) from currents I(k), and Gk(u) = R
y(k) = Gk(u) + η, η ∼ N(0, γ2), k = 1, . . . , K.
1 1.5 2 2.5 3 3.5 4 4.5 5
Figure: True Conductivity.
Parameterization
◮ Continuous level set function. ◮ Lengthscale of level set function. ◮ Smoothness of level set function.
Figure: Five succesive iterations: level set function u.
Figure: Five succesive iterations: thresholded level set function v = S(u).
collaboration with: Cleary, Garbuno-Inigo, Lan, Schneider (2019)
Parameter-Dependent Dynamics
du dt = F(u; θ), u(0) = u0.
Time-Averaged Data
y = GT(θ; u0) = 1 T T ϕ
Central Limit Theorem
GT(θ; u0) = G(θ) + 1 √ T N(0, Σ), y = G(θ) + 1 √ T N(0, Σ).
Governing Dynamics
˙ x = 10 (y − x) ˙ y = r x − y − x z ˙ z = x y − b z ◮ 2-dimensional unknown: θ = [r, b]⊤ ◮ G only available to us approximately via GT. ◮ η only available to us approximately via GT.
EKI: Continuous Time Limit
˙ θ(j) = − 1 J
J
G, G(θ(j)) − y (j)
Γ θ(k),
θ(j)(0) = θ(j) ◮ ϕ(·) is the vector of first and second order moments. ◮ The observation y is simulated from a truth θ∗ ◮ θ∗ = [28, 8/3] ◮ Each y (j) i.i.d ∼ N(y, Γ) ◮ G only available to us approximately via GT.
5 10 15 20 25 30 5 10 15 20 25 30 10
1
10
2
10
3
10
4
10
5
10
6
10
7
Objective Function
25 26 27 28 29 30 r 250 500 750 1000 1250 1500
1 2|y
G( )|2 2.0 2.5 3.0 3.5 4.0 b 250 500 750 1000 1250 1500
1 2|y
G( )|2
◮ Marginals of the objective function Φ(u) = 1
2y − G(θ)2.
◮ In blue, noisy evaluations of the misfit Φ(θ) ◮ In orange: the GP mean from emulation of log Φ. ◮ In shaded orange: is the GP 95% uncertainty bands. ◮ k(θ, θ′) = σ2 exp(−|ℓ, θ − θ′|2) + δ2.
Bayesian Inversion
2.60 2.65 2.70 b 27.8 27.9 28.0 28.1 r 0.000 0.015 0.030 0.045 0.060 0.075 0.090 0.105 0.120
(a) MCMC using noisy Φ(θ)
2.60 2.65 2.70 b 27.8 27.9 28.0 28.1 r 0.000 0.003 0.006 0.009 0.012 0.015 0.018 0.021
(b) MCMC using GP
Figure: MCMC-RWM Sample Density Plots
MCMC-RWM: Comparisons using ∼ 5 × 104 samples
◮ Noisy Φ: 1% acceptance rate; not converged. ◮ GP Φ: 30% acceptance rate; converged.
Forward Model
Moisture Conservation: ∂q ∂t + v · ∇q = −q − qref (T; θ) τq(q, T; θ) Energy Conservation: ∂T ∂t + v · ∇T = T − Tref (q, T; θ) τT(q, T; θ) + RAD + ... ◮ Coupled with mass and momentum conservation equations. ◮ Model computes precipitation rates; initially Pq = PT. ◮ Adjusts reference profiles and/or timescales until Pq = PT. ◮ Enthalpy constraints imposed. ◮ O(105) unknowns ◮ Employs 2 unknown parameters:
◮ θRH: reference relative humidity ◮ θτ: default relaxation timescale
EKI: Discrete Time
θ(j)
n+1 = θ(j) n + C θG n (C GG n
+ Γ)−1 y (j)
n
− G(θ(j)
n )
θ(j)(0) = θ(j)
0 .
◮ ϕ(θ) = [RH(q, T; θ), P(q, T; θ), F(q, T; θ)]⊤
◮ RH: relative humidity ◮ P: daily precipitation rate ◮ ·: denotes longitudinal average ◮ F: longutudinal probability of extreme precipitation (90th %ile)
◮ ϕ(θ) is evaluated over all latitudes at fixed altitude of 5 km ◮ ϕ(θ) ∈ R96
Objective Function
0.600 0.625 0.650 0.675 0.700 0.725
RH
100 200 300 400 500
1 2|y
G( )|2 1 2 3 4 [hrs] 50 100 150 200 250 300
1 2|y
G( )|2
Figure: GP trained from EKI data
In orange, the GP mean from emulation of Φ.
MCMC
0.690 0.695 0.700 0.705 0.710
RH
1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 [hrs] 0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 0.016
(a) GP from uniform samples
0.690 0.695 0.700 0.705 0.710
RH
1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 [hrs] 0.000 0.003 0.006 0.009 0.012 0.015 0.018 0.021 0.024
(b) GP from EKI
Figure: MCMC Sample density plot
RWM with ∼ 106 iterations
collaboration with Albers, Blancquart, Levine, Esmaeilzadeh-Seylabi [6]
Wang D, Chen Y and Cai X 2009 Janji´ c T, McLaughlin D, Cohn S E and Verlaan M 2014
State Space Model
Dynamics Model: vn+1 = Ψ(vn) + ξn, n ∈ Z+ Data Model: yn+1 = Hvn+1 + ηn+1, n ∈ Z+ Noise Structure: ξn ∼ N(0, Σ), ηn ∼ N(0, Γ).
Sequential Optimization Perspective
Predict:
n+1 = Ψ(v (k) n ) + ξ(k) n ,
n ∈ Z+ Model/Data Compromise: J(k)
n (v) = 1
2 C
− 1
2
n+1(v −
v (k)
n+1)2 + 1
2Γ− 1
2 (y (k)
n+1 − Hv)2
Optimize: v (k)
n+1 = argminv∈A J(k) n (v)
Constraint Set: A = {v : Fv = f , Gv g}.
Governing Dynamics
∂ ∂z
∂z
∂t2 = 0, d(H, t) = d0(t), ∂d(0, t)/∂z = 0, d(z, 0) = 0, ∂d(z, 0)/∂t = 0.
Observation Operator
G(θ) = ∂2d(0, t)/∂t2.
Parameterization
c(z) = c0 0 ≤ z ≤ z0 c0
n z0 ≤ z ≤ z1 αc0
n z1 ≤ z ≤ H . θ = (c0, k, z0, n, z1, α).
Constraints
0 ≤ c0 ≤ 1000, 0 ≤ k ≤ 100, 0 ≤ z0 ≤ z1, 0 ≤ n ≤ 1, z0 ≤ z1 ≤ H, 1 ≤ α ≤ 10.
Evolution of Ensemble
Percentage Of Violations
5 10 15 20 Iteration cs0 < 0 k < 0 z0 < 0 n < 0 z0 > z1 , < 1 cs0 > 1000 k > 100 n > 1 z1 > H , > 10 10 20 30 40 50 60
Velocity
10 2 10 4
cs (m/s)
z=H
uy 7 uj=0 7 uj=40
◮ Kalman state estimation may be adapted to parameter estimation. ◮ Continuous time formulation gives insight and new algorithms. ◮ Fit dynamical systems to time-averaged data. ◮ Use Gaussian Process to remove finite time effects. ◮ Applications in climate modeling. ◮ Constraints can be incorporated into EKI. ◮ Applications in seismology. ◮ Tikhonov regulared EKI ([7], 12 noon, Neil Chada, 23/01/19).
[1] R. Kalman A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1960), 35–45. [2] G. Evensen. Data Assimilation: The Ensemble Kalman Filter. Springer, 2006. [3] D.S Oliver, A.C. Reynolds and N Liu Inverse Theory For Petroleum Reservoir Characterization And History Matching. CUP 2008 [4] M. A. Iglesias, K. Law, A. M. Stuart Ensemble Kalman method for inverse problems, Inverse Problems, 29 (4), 2018. [5] N. K. Chada, M. Iglesias, L. Roininen, A. M. Stuart Geometric and Hierarchical Ensemble Kalman Inversion, Inverse Problems, 34 (5), 2018. [6] D. Albers, P.-A. Blancquart, M.D. Levine, E. Esmaeilzadeh Seylabi and A. M. Stuart Ensemble Kalman inversion with constraints, arXiv:1901.05668 [7] N. Chada, A.M. Stuart, X.T. Tong Tikhonov regularizarion within ensemble Kalman inversion, arXiv (to appear).