The Reduced-Order Hybrid Monte Carlo Sampling Smoother

0.42 MB·

by  Ahmed Attia

Computer Science Technical Report
/
CSTR-1 2016
January 5, 2016
6
1
0
Ahmed Attia, Razvan S¸tefa˘nescu, and Adrian
2

n
a
J Sandu

2

]
A
N
“The Reduced-Order Hybrid Monte Carlo
.
s
c
[

1 Sampling Smoother”
v
9
2
1
0 ComputationalScienceLaboratory
0
1. ComputerScienceDepartment
0
6 VirginiaPolytechnicInstituteandStateUniversity
1
:
v Blacksburg,VA24060
i
X
Phone: [phone]
r
a
Fax: [phone]
Email: [email protected],[email protected]
Web: http://csl.cs.vt.edu
InnovativeComputationalSolutions
The Reduced-Order Hybrid Monte Carlo Sampling Smoother
A.Attia∗,R.S¸tefa˘nescu,andA.Sandu
ComputationalScienceLaboratory,
DepartmentofComputerScience,
VirginiaPolytechnicInstituteandStateUniversity,
Blacksburg,Virginia,24060,USA
∗Correspondingauthor
Emailaddress:[email protected],[email protected], [email protected](A.Attia∗,R.S¸tefa˘nescu,andA.Sandu)
PreprintsubmittedtoIJNMF January5,2016
The Reduced-Order Hybrid Monte Carlo Sampling Smoother
A.Attia∗,R.S¸tefa˘nescu,andA.Sandu
ComputationalScienceLaboratory,
DepartmentofComputerScience,
VirginiaPolytechnicInstituteandStateUniversity,
Blacksburg,Virginia,24060,USA
Abstract
Hybrid Monte-Carlo (HMC) sampling smoother is a fully non-Gaussian four-dimensional data assimilation al-
gorithm that works by directly sampling the posterior distribution formulated in the Bayesian framework. The
smootherinitsoriginalformulationiscomputationallyexpensiveduetotheintrinsicrequirementofrunningthe
forwardandadjointmodelsrepeatedly. HerewepresentcomputationallyefficientversionsoftheHMCsampling
smoother based on reduced-order approximations of the underlying model dynamics. The schemes developed
herein are tested numerically using the shallow-water equations model on Cartesian coordinates. The results re-
vealthatthereduced-orderversionsofthesmootherarecapableofaccuratelycapturingtheposteriorprobability
density,whilebeingsignificantlyfasterthantheoriginalfullorderformulation.
Keywords: DataAssimilation,HamiltonianMonte-Carlo,Smoothing,Reduced-OrderModeling,Proper
OrthogonalDecomposition.
Contents
1 Introduction 3
2 DataAssimilation 5
2.1 4D-Vardataassimilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 SmoothingbysamplingandtheHMCsamplingsmoother . . . . . . . . . . . . . . . . . . . . . . 8
3 Four-DimensionalVariationalDataAssimilationwithReduced-OrderModels 10
3.1 Reducedordermodeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.2 Reducedorder4D-Vardataassimilation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4 Reduced-OrderHMCSamplingSmoothers 12
∗Correspondingauthor
Emailaddress:[email protected],[email protected], [email protected](A.Attia∗,R.S¸tefa˘nescu,andA.Sandu)
PreprintsubmittedtoIJNMF January5,2016
4.1 Samplinginthereduced-orderspace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.2 Samplinginthefullspaceusingapproximategradients . . . . . . . . . . . . . . . . . . . . . . . 13
5 PropertiesoftheDistributionsSampledwithReduced-OrderModels 14
5.1 Projectionoftheposteriordistributionforlinearmodelandobservationoperators . . . . . . . . . 16
5.2 Approximatingthelikelihoodfunctionusingreducedordermodels . . . . . . . . . . . . . . . . . 19
6 NumericalResults 22
6.1 TheSWEmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
6.2 Smoothingexperimentalsettings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
6.3 Numericalresults . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.4 Computationalcosts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
7 ConclusionsandFutureWork 26
1. Introduction
Manylarge-scalepredictionproblemssuchasatmosphericforecastingareformulatedasinitialvalueproblems.
Theuncertaintyoftheassociatedmodelinitialconditionscanbedecreasedbycombiningimperfectforecastspro-
ducedbypropagatingthemodeldynamicswithrealmeasurementscollectedatdiscretetimesoveranassimilation
window. Themodelstateandtheobservationsarebothuncertainandcanbemodeledasrandomvariables. The
probabilitydistributiondescribingtheknowledgeabouttheinitialstatesystem,beforeincorporatingobservations
information, is known as the prior distribution. The likelihood function accounts for the discrepancies between
themeasurementsandmodel-predictedobservations. DataassimilationappliesBayes’theoremtoobtainaposte-
riorprobabilitydistribution,namedtheanalysis,thatcharacterizestheknowledgeaboutthesystemstategiventhe
observedmeasurements. Inpractice,duetothestatespacehigh-dimensionalityofmanyrealisticmodels,itisim-
possibletoexactlydescribetheposteriordistributionandseveralassumptionsandapproximationsarenecessary.
Widely accepted assumptions are that the background and the observation errors are characterized by Gaussian
distributions,withnocorrelationsbetweenobservationerrorsatdifferenttimeinstances.
Twofamiliesofmethodologiesaregenerallyfollowedinordertogenerateaccurateestimatesofthetruesys-
temstate.Ensemble-basedstatisticalmethodsseektoapproximatetheposteriorprobabilitydensityfunction(PDF)
based on an ensemble of model states, while the variational approaches estimate the true state of the system by
searching for the state that maximizes the posterior PDF. Among the variational techniques, the 4D-Var method
3
achievesthisgoalbysearchingforalocalminimumofanobjectivefunctioncorrespondingtothenegativeloga-
rithmoftheposteriordistribution. 4D-Varfindsthemaximumposterior(MAP)estimateofthetruestateanddoes
notdirectlyestimatetheuncertaintyassociatedwiththeanalysisstate. Forscenariosincludingnonlinearobserva-
tionoperatorsandstatemodels,theanalysisdistributionisnotGaussian,andthe4D-Varalgorithmmaybetrapped
inalocalminimumofthecostfunctionleadingtoincorrectconclusions.
Accuratesolutionofthenon-Gaussiandataassimilationproblemsrequiresaccountingforallregionsofhigh
probabilityintheposteriordistribution.TherecentlydevelopedhybridMonte-Carlo(HMC)samplingsmoother[4,
5]isafourdimensionaldataassimilationschemedesignedtosolvethenon-Gaussiansmoothingproblembysam-
plingfromtheposteriordistribution.ItreliesonanacceleratedMarkovchainMonte-Carlo(MCMC)methodology
whereaHamiltoniansystemisusedtoformulateproposaldensities. TwoissuesareimportantwhenusingHMC
sampling strategy. First, it requires the formulation of the posterior negative-log function gradient. Secondly,
theinvolvedHamiltoniansystemispropagatedusingasymplecticintegratorwhoseparametersmustbecarefully
tunedinordertoachievegoodperformance.
Whileproducingconsistentdescriptionoftheupdatedsystemuncertainty(e.g., theanalysiserrorcovariance
matrix),theoriginalformulationoftheHMCsamplingsmootheriscomputationallyexpensivewhencomparedto
the4D-Varapproach. Thisisduetothelargenumberofgradientevaluationsrequired,whichtranslatesintomany
forwardandadjointmodelsruns. Inthecaseoflargescaleproblemsthecomputationalcostbecomesprohibitive.
In this present study we propose a practical solution by approximating the gradient using information obtained
fromlower-dimensionalsubspacesviamodelreduction.
Reduced order modeling refers to the development of low-dimensional systems that represent the important
characteristicsofahigh-dimensionalorinfinitedimensionaldynamicalsystem. Typicallythisisachievedbypro-
jectingmodeldynamicsontoalowerdimensionalspaces. Constructionoflowrelevantmanifoldscanbeachieved
usingthereducedbasismethod[8,26,43,47,18,36],dynamicmodedecomposition[46,52,62,10]andProper
Orthogonal Decomposition (POD) [32, 37, 30, 38]. The latest is the most prevalent basis selection method for
nonlinear problems. Data analysis using POD and method of snapshots [54, 55, 56] is conducted to extract ba-
sisfunctions, fromexperimentaldataordetailedsimulationsofhigh-dimensionalsystems, forsubsequentusein
Galerkin projections that yield low dimensional dynamical models. By coupling POD with empirical interpola-
tionmethod(EIM)[8], discretevariantDEIM[14, 16, 15, 20]orbestpointsinterpolationmethod[41], onecan
obtainfastapproximationsofreducedordernonlineartermsindependentofthedimensionofhigh-fidelityspace.
Othersuchapproachesincludemissingpointestimation[3]andGauss-Newtonwithapproximatedtensors[12,13]
relyinguponthegappyPODtechnique[23].
The present manuscript develops practical versions of the HMC sampling smoother using approximate dy-
4
namicalinformationobtainedviaPOD/DEIMreducedordermodels[57,59,19]. Tworeducedordervariantsare
proposeddifferentiatedbythechoiceofthesamplingspaceusedtogeneratetheproposals. Inthefirstcaseweare
sampling in a reduced order space while the second approach samples directly from the high fidelity space. For
bothscenariosthenegativelogarithmoftheposteriordistributionreassemblesoneoftheflavorsofreducedorder
4D-Var objective functions [24, 11, 63], where the prior term is either estimated in the reduced or full space re-
spectively. Whilereducedordermodelsareemployedtoestimatetheposteriordistributionnegativelogarithmand
it’sgradient,theassociatedHamiltoniansystemgeneratesproposalsinareducedorhigh-fidelityspacedepending
onthenatureofthesamplingspace.
Thechoiceofreducedordermanifoldsiscrucialfortheaccuracyofthereducedsamplers. Thesmoothersmay
sufferfromthefactthatthebasiselementsarecomputedfromareferencetrajectorycontainingfeatureswhichare
quite different from those of the current proposal trajectory. To overcome the problem of unmodelled dynamics
in the POD-basis we propose to update the basis from time to time according to the current state similarly as in
thecaseofadaptivereducedorderoptimization[1,45,35]. Inspiredbytherecentadvancesinreducedorderdata
assimilation field where it was shown that accurate reduced order Karush-Kuhn-Tucker conditions with respect
totheirfullordercounterpartshighlyincreasetheaccuracyoftheobtainedanalysis[61],wechosetoupdatethe
reducedordermanifoldsbasedonhigh-fidelityforwardandadjointtrajectoriescorrespondingtocurrentproposal
as well as the associated gradient of the negative logarithm of the high-fidelity posterior distribution. The basis
is refreshed once after several HMC smoothers iterations. The numerical results using swallow-water equations
modeloncartesiancoordinatesrevealthatthereduced-orderversionsofthesmootherareaccuratelycapturingthe
posteriorprobabilitydensity,whilebeingsignificantlyfasterthantheoriginalfullorderformulation.
Thepaperisorganizedasfollows. Section2reviewsthefour-dimensionaldataassimilationproblemandthe
originalHMCsmoother. Section3reviewsreduced-ordermodeling,andintroducesthereducedorder4D-Vardata
assimilationframework. Section4presentsarigorousdescriptionoftheproposedversionsoftheHMCsmoother.
In section 5 several theoretical properties of the distributions sampled with reduced order models are derived.
NumericalresultsarepresentedinSection6whileconclusionsaredrawninSection7.
2. DataAssimilation
One can solve a data assimilation (DA) problem by describing the posterior distribution in the Bayesian for-
malism given the prior and the likelihood function. In the four-dimensional DA (4DDA) context, the main goal
is to describe the posterior distribution of system state at the initial time of a specific assimilation window. The
assimilation window is defined based on the availability of observations at specific discrete time instances. The
knowledgeaboutthesystemattheinitialtimet definesthedistribution(theprior)Pb(x )ofthestatex ∈ RNvar
0 0 0
5
before absorbing knowledge presented by the observations captured over the assimilation window. Based on a
set of observations {yk = y[tk] ∈ Rm}k=0,1,...,Nobs, at the discrete time points {tk}k=0,1,...,Nobs in the interval [t0,tF],
the sampling distribution (likelihood function) is defined as P(y0,y1,...,yNobs|x0), and the posterior distribution
describing the updated information about the model state at the initial time, given the measurements, takes the
generalform
P(x0|y0,y1,...,yNobs)= Pb(x0P)P(y(0y,0y,1y,1.,....,.y,NyoNbso)bs|x0)). (1)
FullyandaccuratelydescribingthisgeneralposteriorprobabilitydensityfunctioninDAliteratureisanintractable
problem, and usually simplifying assumptions are generally made. As we mentioned in Section 1, a frequent
supposition is that the background and observation errors are characterized by Gaussian distributions. If the ob-
servationsareassumedtobeindependentfromthemodelstates, andtheassociatederrorcharacteristicsofthese
observationsarenotcorrelatedintime,theposteriordistributiontakestheform
Ä ä
Pa(x0)=P(x0|y0,y1,...,yNobs)∝exp −J(x0) , (2a)
Nobs
J :RNvar →R, J(x )= 1(cid:107)x −xb(cid:107)2 + 1(cid:88)(cid:107)y −H (x )(cid:107)2 , (2b)
0 2 0 0 B−01 2 k=0 k k k R−k1
where xb = xb[t ] is a background/forecast state, the matrices B and R , k=0,1,...,Nobs, are the covariance
0 0 0 k
matrices associated with the background and measurement errors respectively. The observation operator H :
k
RNvar → Rm, at time instance t , maps a given state x ∈ RNvar to the observation space. The dimension of the
k k
observationspaceisusuallymuchsmallerthanthesizeofthestatespace,thatism(cid:28)Nvar. Theassociatednorms
overthestateandobservationspacesaredefinedas:
(cid:107)a −b(cid:107)2 =(a−b)TC(a−b), (3)
C
wherea, bbelongtoeitherRNvarorRmandCisamatrixofRNvar×NvarorRm×mdimensions.
The model state x = x[t ], is obtained from the initial state x , by propagating the model dynamics to time
k k 0
pointt ,i.e.
k
x =M (x )=M (x ), (4)
k t0→tk 0 0,k 0
where,M :RNvar →RNvar, k=1,...,Nobsrepresentsthediscretizedmathematicalmodelreflectingthestate
t0→tk
dynamics. ThefunctionJ givenin(2)isquadratic,andconsequentlythedistribution(2)isGaussiandistribution,
onlyifthestatesx arelinearlyrelatedtothesysteminitialconditionx ,andtheobservationsy arelinearlyrelated
k 0 k
to the model states x . In virtually all practical settings, the time dependent models are complicated resulting in
k
6
a nonlinear bond between states at different time instances. More difficulty is added since for example, in the
fieldofatmosphericsciences,highly-nonlinearobservationoperatorsareoftenconstructedtorelatenewtypesof
measurementstostatevariables.
2.1. 4D-Vardataassimilation
The strongly-constrained 4D-Var DA scheme is an optimization algorithm that searches for the maximum
aposteriori estimate (MAP) by seeking a local minimizer of the cost function (2b), constrained by the model
equation(4). Precisely,theconstrainedoptimizationproblemisdefinedas
Nobs
1 1(cid:88)
minJ(x )= (cid:107)x −xb(cid:107)2 + (cid:107)y −H (x )(cid:107)2 ,
x0 0 2 0 0 B−01 2 k=0 k k k R−k1 (5)
x =M (x ) k=1,...,Nobs.
k t0→tk 0
Gradient-basedschemesaregenerallyemployedtofindthecorrespondinginitialconditionswhichinthecontrol
theorylanguageareknownascontrolvariables.UsingtheLagrangemultipliertechniquetheconstrainedoptimiza-
tionproblem(5)isredesignedasanunconstrainedminimizationproblemwiththeassociatedfirstordernecessary
optimalityconditions:
Forwardmodel: x =M (x ), k=1,...,Nobs, (6a)
k t0→tk 0
Adjointmodel: λNobs =HTNobsR−N1obs(cid:0)yNobs−HNobs(xNobs)(cid:1), (6b)
λk =MTt0→tk+1λk+1+HTkR−k1(cid:0)yk−H(xk)(cid:1), k=Nobs−1,..,0, (6c)
Costfunctiongradient: ∇ J(x )=−B−1(cid:0)xb−x (cid:1)−λ =0. (6d)
x0 0 0 0 0 0
The Jacobians of the model and observation operators are denoted by M , k = 0,1,...,Nobs−1 and
t0→tk+1
H , k = 0,1,...,Nobs while the adjoint solution λ ∈ RNvar, k = 0,1,...,Nobs, provides an efficient way to
k k
compute the gradient (6d). The nonlinear optimization procedure is computationally expensive and either low-
resolutionmodels(incremental4D-Var[58]),oralternativelyreduced-ordermodelsareused[61]toalleviatethis
drawback.Itiswell-knownthat4D-Vardoesnotinherentlyprovideameasureoftheuncertaintyabouttheupdated
state (e.g., analysis error covariance matrix) and usually hybrid methods are considered to account for this type
of information. This approach results in some inconsistency between the analysis state and the analysis error
covariancematrixespeciallywhentheyareobtainedusingdifferentalgorithms.
7
2.2. SmoothingbysamplingandtheHMCsamplingsmoother
Monte-Carlosmoothingreferstotheprocessofrepresenting/approximatingtheposteriordistribution(2)using
anensembleofmodelstatessampledfromthatposterior. EnsembleKalmansmoother(EnKS)[22]isanextension
of the well-known ensemble Kalman Filter [31] to the case where observation are assimilated simultaneously.
EnKSproducesaminimum-varianceunbiasedestimate(MVUE)ofthesystemstatebyestimatingtheexpectation
oftheposteriorE [x ]usingthemeanofanensembleofstates. ThestrictGaussianityandlinearityassumptions
Pa 0
imposedbyEnKS,usuallyresultinpoorperformanceofthesmoother.
Pure sampling of the posterior distribution (1) using a Markov Chain Monte-Carlo (MCMC) [39, 40] tech-
nique is known in theory to provide more accurate estimates without strictly imposing linearity or Gaussianity
constraints. MCMC is a family of Monte-Carlo schemes tailored to sample a given distribution (up to a propor-
tionality constant), by constructing a Markov chain whose stationary distribution is set as the target distribution.
Bydesign, anMCMCsampleris guaranteedtoconvergetoitsstationarity. However, thechoiceoftheproposal
density, the convergence rate, the acceptance rate, and the correlation level among sampled points are the main
buildingblocksofMCMCresponsibleforitsperformanceandefficiency. PracticalapplicationofMCMCrequires
developingacceleratedchainsthosecanattainstationarityfast,andthenexplorethestatespaceefficientlyinvery
fewsteps. OneoftheMCMCsamplersmainlydesignedforcomplicatedPDFsandlargedimensionalspacesisthe
Hamiltonian/Hybrid Monte-Carlo sampler (HMC). HMC was firstly presented in [21] as an accelerated MCMC
samplingalgorithm. Thesamplermainlyusesinformationaboutthegeometryoftheposteriortoguideitsstepsin
ordertoavoidrandomwalkbehaviourandvisitmorefrequentlyregionswithhighprobabilitywiththecapability
ofjumpingbetweenseparatedmodesofthetargetPDF.
HMCsampling. AsallotherMCMCsamplers,HMCsamplesfromaPDF
π(x)∝exp(−JN(x)); x∈RNvar, (7)
whereexp(−JN(x))istheshapefunctionofthedistribution,andJN : RNvar → RisthePDFnegative-log. The
power of MCMC, and consequently HMC, is that only the shape function, or alternatively the negative-log, is
needed, while thescalingfactor isnot strictlyrequiredas inthe caseofstandard applicationof Bayes’theorem.
HMCworksbyviewingxasapositionvariableinanextendedphasespaceconsistingofpoints(p,x) ∈ R2Nvar,
wherep∈RNvarisanauxiliarymomentumvariable. TheHamiltoniandynamicsismodeledbythesetofordinary
differentialequations(ODEs):
dx
=∇ H,
dt p
(8)
dp
=−∇ H,
dt x
8
where H = H(p,x)istheconstanttotalenergyfunctionofthesystem,knownasthetheHamiltonianfunctionor
simplytheHamiltonian. AstandardformulationoftheHamiltonianinthecontextofHMCis:
1
H(p,x)= pTM−1p+ JN(x) , (9)
(cid:32)(cid:32)(cid:32) (cid:32)(cid:32)(cid:32)
2(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:124) (cid:123)(cid:122) (cid:125)
(cid:124) (cid:123)(cid:122) (cid:125) potentialenergy
kineticenergy
whereM∈RNvar×Nvar isapositivedefinitematrixknownasthemassmatrix. Thisparticularformulationleadsto
acanonicaldistributionofthejointstate(p, x)proportionalto:
Å ã
1
exp(−H(p, x))=exp − pTM−1p π(x). (10)
2
The exact flow Φ : R2Nvar → R2Nvar; Φ (cid:0)p[0],x[0](cid:1) = (cid:0)p[T],x[T](cid:1) describes the time evolution of
T T
the Hamiltonian system (8) and is practically approximated using a numerical integrator that is symplectic as
φ :R2Nvar →R2Nvar; φ (cid:0)p[0],x[0](cid:1)≈(cid:0)p[T],x[T](cid:1). Theuseofasymplecticnumericalintegratortoapproxi-
T T
matetheexactHamiltonianflowresultsinchangesoftotalenergy. Traditionalwisdomrecommendssplittingthe
pseudo-timestepT intomsmallerstepsofsizehinordertosimulatetheHamiltoniantrajectorybetweenpoints
(cid:0) (cid:1) (cid:0) (cid:1)
p[0],x[0] , p[T],x[T] more accurately. The symplectic integrator of choice is position (or velocity) Verlet.
New higher order symplectic integrators were proposed recently and tested in the context of filtering for data
assimilation[50,6]. OnestepofsizehofthepositionVerlet[51,50]integratorisdescribesasfollows
h
x[h/2]=x[0]+ M−1p[0], (11a)
2
p[h]=p[0]−h∇ J(x[h/2]), (11b)
x
h
x[h]=x[h/2]+ M−1p[h]. (11c)
2
While,themassmatrixisauser-definedparameter,itcanbedesignedtoenhancetheperformancethesampler[6].
Thestepparametersofthesymplecticintegratorm, hcanbeempiricallychosenbymonitoringtheacceptancerate
inapreprocessingstep.Specifically,theparametersoftheHamiltoniantrajectorycanbeempiricallyadjustedsuch
astoachieveaspecificrejectionrate. Generallyspeaking,thestepsizeshouldbechosentoachievearejectionrate
between25%and30%,andthenumberofstepsshouldgenerallybelarge[40].
AdaptiveversionsofHMChavebeenalsoproposedwiththecapabilityofadjustingit’sstepparameters. No-
U-Turnsampler(NUTS)[29]isaversionofHMCcapableofautomaticallytuningitsparameterstoprohibitthe
samplerfromretracingitsstepsalongtheconstructedHamiltoniantrajectory. AnotherHMCsamplerthattunesits
parametersautomaticallyusingthird-orderderivativeinformationistheRiemannmanifoldHMC(RMHMC)[25].
9