Gappy POD and GNAT
David Amsallem
Stanford University
February 05, 2014
David Amsallem Gappy POD and GNAT
Outline
• Nonlinear partial differential equations
• An issue with the model reduction of nonlinear equations
• The gappy proper orthogonal decomposition
• The discrete empirical interpolation method (DEIM)
• The Gauss-Newton with approximated tensors method (GNAT)
• Two applications
David Amsallem Gappy POD and GNAT
Nonlinear PDE
• Parametrized partial differential equation (PDE)
L(W, x, t; µ) = 0
• Associated boundary conditions
B(W, xBC , t; µ) = 0
• Initial condition
W0 (x) = WIC (x, µ)
• W = W(x, t) ∈ Rq : state variable
• x ∈ Ω ⊂ Rd , d ≤ 3: space variable
• t ≥ 0: time variable
• µ ∈ D ⊂ Rp : parameter vector
David Amsallem Gappy POD and GNAT
Discretization of nonlinear PDE
• The PDE is then discretized in space by one of the following methods
• Finite Differences approximation
• Finite Element method
• Finite Volumes method
• Discontinuous Galerkin method
• Spectral method....
• This leads to a system of Nw = q × Nspace ordinary differential equations
(ODEs)
dw
= f (w, t; µ)
dt
in terms of the discretized state variable
w = w(t; µ) ∈ RNw
with initial condition w(0; µ) = w(µ)
• This is the high-dimensional model (HDM)
David Amsallem Gappy POD and GNAT
Model reduction of nonlinear equations
• High-dimensional model (HDM)
dw
(t) = f (w(t), t)
dt
• Reduced-order modeling assumption using a reduced basis V
w(t) ≈ Vq(t)
• q(t): reduced (generalized) coordinates
• Inserting in the HDM equation
dq
V (t) ≈ f (Vq(t), t)
dt
• Nw equations in terms of k unknowns q
• Galerkin projection
dq
(t) = VT f (Vq(t), t)
dt
David Amsallem Gappy POD and GNAT
Issue with the model reduction of nonlinear equations
• Galerkin projection
dq
(t) = VT f (Vq(t), t)
dt
• k equations in terms of k unknowns
• To evaluate fk (Vq(t), t):
1 Compute w(t) = Vq(t)
2 Evaluate f (Vq(t), t)
3 Left-multiply by VT : VT f (Vq(t), t)
• The computational cost associated with these three steps scales linearly
with the dimension Nw of the HDM
• Hence no significant speedup can be expected when solving the
projection-based ROM
David Amsallem Gappy POD and GNAT
The Gappy POD
• First applied to face recognition (Emerson and Sirovich, ”Karhunen-Loeve
Procedure for Gappy Data” 1996)
• Procedure
1 Build a database of m faces (snapshots)
2 Construct a POD basis V for the database
3 For a new face f , record a few pixels f1 , · · · , fn
4 Using the POD basis V, approximately reconstruct the new face f
David Amsallem Gappy POD and GNAT
The Gappy POD
• First applied to face recognition (Emerson and Sirovich, ”Karhunen-Loeve
Procedure for Gappy Data” 1996)
David Amsallem Gappy POD and GNAT
The Gappy POD
• Other applications
• Flow sensing and estimation (Willcox, 2004)
• Flow reconstruction
• Nonlinear model reduction
David Amsallem Gappy POD and GNAT
Nonlinear function approximation by gappy POD
• Approximation of the nonlinear function f in
dq
= VT f (Vq(t), t)
dt
• The evaluation of all the entries in the vector f (·, t) is expensive (scales
with Nw )
• Only a small subset of these entries will be evaluated (gappy approach)
• The other entries will be reconstructed either by interpolation or a
least-squares strategy using a pre-computed specific reduced-order
basis Vf
• The solution space is still reduced by any preferred model reduction
method (by POD for instance)
David Amsallem Gappy POD and GNAT
Nonlinear function approximation by gappy POD
A complete model reduction method should then provide algorithms for
• Selecting the evaluation indices I = {i1 , · · · , iNi }
• Selecting a reduced-order bases Vf for the nonlinear function
• Reconstructing the complete approximated nonlinear function vector
f̂ (·, t)
David Amsallem Gappy POD and GNAT
Construction of a POD basis for f
• Construction of a POD basis Vf of dimension kf
1 Collection of snapshots for the nonlinear function from a transient simulation
F = [f (w(t1 ), t1 ), · · · , f (w(tmf ), tmf )] ∈ RNw ×mf
2 Singular value decomposition
F = Uf Σf ZTf
3 Basis truncation (kf mf )
Vf = [uf ,1 , · · · , uf ,kf ]
David Amsallem Gappy POD and GNAT
Reconstruction of an approximated nonlinear function
• Assume ki indices have been chosen
I = {i1 , · · · , iki }
• The choice of indices will be specified later
• Consider the Nw -by-ki matrix
h i
P = ei1 , · · · , eiki
• At each time t, for a given value of the state w(t) = Vq(t), only the
entries in the function f corresponding to those indices will be evaluated
fi1 (w(t), t)
..
PT f (w(t), t) =
.
fiki (w(t), t)
• This is cheap if ki Nw
• Usually only a subset of the entries in w(t) will be required to construct
that vector (case of sparse Jacobian)
David Amsallem Gappy POD and GNAT
Discrete Empirical Interpolation Method
• Case where ki = kf : interpolation
• Idea: fˆij (w, t) = fij (w, t), ∀w ∈ RNw , ∀j = 1, · · · , ki
• This means that
PT f̂ (w(t), t) = PT f (w(t), t)
• Remember that f̂ (·, t) belongs to the span of the vectors in Vf , that is
f̂ (Vq(t), t) = Vf fr (q(t), t)
• Then
PT Vf fr (q(t), t) = PT f (Vq(t), t)
• Assuming PT Vf is nonsingular
fr (q(t), t) = (PT Vf )−1 PT f (Vq(t), t)
• In terms of f̂ (·, t):
f̂ (·, t) = Vf (PT Vf )−1 PT f (·, t) = ΠVf ,P f (·, t)
• This results in an oblique projection of the full nonlinear vector
David Amsallem Gappy POD and GNAT
Oblique projection of the full nonlinear vector
f̂ (·, t) = Vf (PT Vf )−1 PT f (·, t) = ΠVf ,P f (·, t)
• ΠV,W = V(WT V)−1 WT : oblique projector onto V orthogonally to W
David Amsallem Gappy POD and GNAT
Reduced-order dynamical system
• Case where ki > kf : least-squares reconstruction
• Idea: fˆij (w, t) ≈ fij (w, t), ∀w ∈ RNw , ∀j = 1, · · · , Ni in the least squares
sense
• Idea: minimize
fr (q(t)) = argmin kPT Vf yr − PT f (Vq(t), t)k2
yr
T ki ×kf
• Note that M = P Vf ∈ R is a skinny matrix
• One can compute its singular value decomposition
M = UΣZT
• The left inverse of M is then defined as
M† = ZΣ† UT
where Σ† = diag( σ11 , · · · , σ1r , 0, · · · , 0) if Σ = diag(σ1 , · · · , σr , 0, · · · , 0) with
σ1 ≥ · · · σr > 0
• Then
f̂ (q(t)) = Vf (PT Vf )† PT f (Vq(t), t)
David Amsallem Gappy POD and GNAT
Greedy function sampling
• This selection takes place after the vectors [vf,1 , · · · , vf,kf ] have already
been computed by POD
• Greedy algorithm (Chaturantabut et al. 2010):
1: [s, i1 ] = max{|vf,1 |}
2: Vf = [vf,1 ], P = [ei1 ]
3: for l = 2 : kf do
4: Solve PT Vf c = PT vf,l for c
5: r = vf,l − Vf c
6: [s, il ] = max{|r|}
7: Vf = [Vf , vf,l ], P = [P, eil ]
8: end for
David Amsallem Gappy POD and GNAT
Model reduction at the fully discrete level
• Semi-discrete level: dw
dt (t) = f (w(t), t)
• Subspace assumption w(t) ≈ Vq(t)
dq
V (t) ≈ f (Vq(t), t)
dt
• Fully discrete level (implicit, backward Euler scheme)
q(n+1) − q(n)
V (n)
≈ f Vq(n+1) , t(n+1)
∆t
• Fully discrete residual
(n+1) (n+1) q(n+1) − q(n)
(n+1) (n+1)
rD (q )=V − f Vq , t
∆t(n)
• Model reduction by least-squares (Petrov-Galerkin)
(n+1)
q(n+1) = argmin krD (y)k2
y
• rD (q(n+1) ) is nonlinear ⇒ use the gappy POD idea
David Amsallem Gappy POD and GNAT
Gappy POD at the fully discrete level
• Gappy POD procedure for the fully discrete residual rD
• Algorithm
1 Build a reduced basis Vr ∈ RNw ×kr for rD
2 Construct a sample mesh I (indices i1 , · · · , iki ) by a greedy procedure
3 Consider the gappy approximation
“ ”†
(n+1)
rD (q(n+1) ) ≈ Vr rkr (q(n+1) ) ≈ Vr PT Vr PT r(n+1) (Vq(n+1) )
4 Solve
q(n+1) = argmin kVr rkr (y)k2
y
= argmin krkr (y)k2
y (1)
‚“ ‚
‚ T ”† T (n+1) ‚
= argmin ‚
‚ P Vr P r (Vy) ‚
‚
y 2
David Amsallem Gappy POD and GNAT
Gauss-Newton for nonlinear least squares problem
• Nonlinear least squares problem miny kr(y)k2
• Equivalent function to be minimized
f (y) = 0.5kr(y)k22 = r(y)T r(y)
• Gradient
∇f (y) = J(y)T r(y)
∂r
where J(y) = ∂y (y)
• Iterative solution using Newton’s method y(k+1) = y(k) + ∆y(k+1) with
∇2 f (y(k) )∆y(k+1) = −∇f (y(k) )
• What is ∇2 f (y)?
N
X ∂ 2 ri
∇2 f (y) = J(y)T J(y) + (y)ri (y)
i=1
∂y2
• Gauss-Newton method
∇2 f (y) ≈ J(y)T J(y)
David Amsallem Gappy POD and GNAT
Gauss-Newton for nonlinear least squares problem
• Gauss-Newton method y(k+1) = y(k) + ∆y(k+1) with
J(y(k) )T J(y(k) )∆y(k+1) = −J(y(k) )T r(y(k) )
• This is the normal equation for
∆y(k+1) = argmin
J(y(k) )z + r(y(k) )
z 2
• QR decomposition of the Jacobian
J(y(k) ) = Q(k) R(k)
• Equivalent solution using the QR decomposition (assume R(k) is full
rank)
−1 T
∆y(k+1) = −J(y(k) )† r(y(k) ) = − R(k) Q(k) r(y(k) )
David Amsallem Gappy POD and GNAT
Gauss-Newton with Approximated Tensors
• GNAT = Gauss-Newton + Gappy POD
• Minimization problem
†
min
PT Vr PT r(n+1) (Vy)
y 2
• Jacobian †
JD (y) = PT Vr PT J(n+1) (Vy)
• Define a small dimensional operator (constructed offline)
†
A = PT Vr
• Least-squares problem at iteration k
∆y(k) = argmin
APT J(n+1) (Vy(k) )Vz + APT r(n+1) (Vy(k) )
z 2
(k) (k) T (n+1)
• GNAT solution using QR decomposition Q R = AP J (Vy(k) )V
−1 T
∆y(k) = − R(k) Q(k) APT r(n+1) (Vy(k) )
David Amsallem Gappy POD and GNAT
Gauss-Newton with Approximated Tensors
• Further developments
• Concept of reduced mesh
• Concept of output mesh
• Error bounds
• GNAT using Local reduced bases
• More details in Carlberg et al., JCP 2013
David Amsallem Gappy POD and GNAT
Application 1: compressible Navier-Stokes equations
• Flow past the Ahmed body (automotive industry benchmark)
• 3D compressible Navier-Stokes equations
• Nw = 1.73 × 107
• Re = 4.48 × 106 , M∞ = 0.175 (216km/h)
• More details in Carlberg et al., JCP 2013
David Amsallem Gappy POD and GNAT
Application 1: compressible Navier-Stokes equations
• Model reduction (POD+GNAT): k = 283, kf = 1, 514, ki = 2, 268
Method CPU Time Number of CPUs Relative Error
Full-Order Model 13.28 h 512 –
ROM (GNAT) 3.88 h 4 0.68%
David Amsallem Gappy POD and GNAT
Application 2: design-optimization of a nozzle
• Full model: Nw = 2, 048, p = 5 shape parameters
• Model reduction (POD+DEIM): k = 8, kf = 20, ki = 20
p1
p2
p3
p4
p5
0
x
L
min kM (w(µ)) − Mtarget k2
µ∈R5
s.t. f (w(µ)), µ) = 0
David Amsallem Gappy POD and GNAT
Application 2: design-optimization of a nozzle
Method Offline CPU Time Online CPU Time Total CPU Time
Full-Order Model – 78.8 s 78.8 s
ROM (GNAT) 5.08 s 4.87 s 9.96 s
1.2 Target
0.75 I nitial Guess
O ptim ization (Herm ite RB F )
0.7 1.1
0.65
1
0.6
M (x)
A(x)
0.55 0.9
0.5 0.8
0.45
0.7
0.4
Target
0.35 I nitial Guess 0.6
O ptim ization (Herm ite RB F )
0.3
0.5 1 1.5 0.5 1 1.5
x x
David Amsallem Gappy POD and GNAT
References
• R. Everson, L. Sirovich. Karhunen–Loeve procedure for gappy data.
Journal of the Optical Society of America A 1995; 12(8):1657–1664.
• M. Barrault, Y. Maday, N.C. Nguyen, A.T. Patera. An empirical
interpolation method: application to efficient reduced basis discretization
of partial differential equations. Comptes Rendus de lAcademie des
Sciences Paris 2004; 339:667–672.
• K. Willcox. Unsteady flow sensing and estimation via the gappy proper
orthogonal decomposition. Computers and Fluids 2006; 35:208–226.
• S. Chaturantabut, D.C. Sorensen. Nonlinear Model Reduction via
Discrete Empirical Interpolation. SIAM Journal on Scientific Computing
2010; 32:2737–2764.
David Amsallem Gappy POD and GNAT
References (continued)
• K. Carlberg, C. Bou-Mosleh, C. Farhat. Efficient non-linear model
reduction via a least-squares Petrov–Galerkin projection and
compressive tensor approximations. International Journal for Numerical
Methods in Engineering 2011; 86(2):155–181.
• K. Carlberg, C. Farhat, J. Cortial, D. Amsallem. The GNAT method for
nonlinear model reduction: effective implementation and application to
computational fluid dynamics and turbulent flows. Journal of
Computational Physics 2013; 242(C): 623–647.
• D. Amsallem, M. Zahr, Y. Choi, C. Farhat. Design Optimization Using
Hyper-Reduced-Order Models, submitted for publication – preprint
http://www.stanford.edu/∼amsallem/HROM Opt.pdf
David Amsallem Gappy POD and GNAT
