arXiv:1307.5076v1 [] 18 Jul 2013

Computational Science
Laboratory Technical Report
CSL-TR-3-2013
July 22, 2013
Alexandru Cioaca,
Adrian Sandu, and Eric de Sturler
Efficient methods for computing
observation impact
in 4D-Var data assimilation
Computational Science Laboratory
Computer Science Department
Virginia Polytechnic Institute and State University
Blacksburg, VA 24060
Phone: (540)-231-2193
Fax: (540)-231-6075
Email: sandu@cs.vt.edu
Web: http://csl.cs.vt.edu

Innovative Computational Solutions

Low-rank Approximations
for Computing Observation Impact
in 4D-Var Data Assimilation
Alexandru Cioaca, Adrian Sandu

Abstract
We present an efficient computational framework to quantify the impact of
individual observations in four dimensional variational data assimilation. The
proposed methodology uses first and second order adjoint sensitivity analysis,
together with matrix-free algorithms to obtain low-rank approximations of observation impact matrix. We illustrate the application of this methodology to
important applications such as data pruning and the identification of faulty
sensors for a two dimensional shallow water test system.

Preprint submitted to Elsevier

July 22, 2013

Contents
1 Introduction

1

2 4D-Var Data Assimilation
2.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Computational aspects . . . . . . . . . . . . . . . . . . . . . . . .

2
2
3

3 Observation Impact
3.1 4D-Var sensitivity to observations .
3.2 Observation impact matrix . . . .
3.2.1 Main features . . . . . . . .
3.2.2 Computational ingredients .

3
3
5
5
6

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

4 Efficient Computation of Observation Sensitivity and Impact
4.1 An iterative (serial) approach for matrix-free low-rank approximations with SVD . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 An ensemble-based (parallel) approach for matrix-free low-rank
approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
9

5 Applications
10
5.1 Test problem: shallow water equations . . . . . . . . . . . . . . . 10
5.2 Data assimilation setup . . . . . . . . . . . . . . . . . . . . . . . 11
5.3 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . 12
5.3.1 Validating the low-rank approximations of the observation
impact matrix . . . . . . . . . . . . . . . . . . . . . . . . 12
5.3.2 Observation impact . . . . . . . . . . . . . . . . . . . . . 13
5.3.3 Pruning observations based on sensitivity values . . . . . 17
5.3.4 Identifying faulty data using observation sensitivity . . . . 18
6 Conclusions and future work

19

2

1. Introduction
Data assimilation is a dynamic-data driven application that integrates information from physical observations with numerical model predictions. This
paper describes a systematic approach to quantify the contribution of each observation data point in improving the model state estimates. The focus is on
the four-dimensional variational (4D-Var) data assimilation approach, which fits
a model trajectory against time-distributed observations, to obtain maximum
likelihood estimates of the initial state and model parameters.
Computing the “observation impact” (OI) means quantifying the individual
(or group) contributions of “assimilated” data points to reduce the uncertainty
in state or model parameters. OI provides a numerical measure to distinguish
between crucial and redundant information, and benefits applications such as
big data pruning, error detection/correction, intelligent sensor placement, and
other decision-making processes.
There are several approaches available in the scientific literature to compute
OI. In early studies OI was associated with the energy of perturbations along
the dominant directions of maximum error growth [1, 2]. Other approaches
study observation impact using metrics from information theory [3, 4], statistical design [5], and robust control [6]. More recently, observation impact has
been assessed through adjoint-sensitivity analysis, as a counterpart measure to
observation sensitivity [7–9]. The current state-of-the-art uses the second-order
sensitivity equations of 4D-Var system[10–12], a powerful approach followed by
the present research study. The efficient implementation of this methodology
is hindered by several computational limitations related to the calculation of
second-order derivative approximations.
This work develops a systematic and efficient approach to compute OI in a
sensitivity analysis framework. Observation impact is formalized as a sensitivity matrix which maps changes in observation space to associated changes in
solution space. This matrix is usually not available in explicit (full) form. Our
computational approach makes smart use of tangent linear and adjoint models
[13–15] to propagate first and second order sensitivities through the data assimilation system, and obtain products of the OI matrix with user defined vectors.
Matrix-free linear algebra tools are necessary for this. Eigenvalue and singular
value decompositions are used to obtain low-rank approximations that capture
the main features of observation impact, while significantly reducing the computational burden. As of recently, low-rank approximations of matrices have
become very popular in image processing [16], information retrieval [17], and
machine learning [18] to extract correlations and remove noise from data. Two
alternative ways to compute a low-rank approximation for observation impact
are given, one of serial nature and one highly parallel. Several applications of OI
are exemplified using the two dimensional shallow water equations test problem.
The remaining part of the paper is structured as follows. Section 2 reviews
the formulation of 4D-Var data assimilation. Section 3 reviews the 4D-Var
sensitivity equations necessary for the observation impact matrix (following the
original derivation [11]). We discuss the structural and functional properties of
1

the impact matrix and the required numerical tools for evaluating its action.
Section 4 introduces low-rank approximations of the observation impact matrix.
Numerical experiments to exemplify and validate the observation impact matrix
are presented in Section 5. Conclusions and future directions of research are
outlined in Section 6.
2. 4D-Var Data Assimilation
2.1. Formulation
Data assimilation is the process by which model predictions are constrained
with real measurements [19, 20]. It combines three sources of information: an
a priori estimate of the initial state of the system (“background”), knowledge
of the physical laws governing the behavior of the system (captured by a numerical model), and measurements of the real system (“observations”). Data
assimilation solves the inverse problem for improving estimates of model states,
initial and boundary conditions, or various system parameters.
Four-dimensional variational (4D-Var) assimilation is formulated as a PDEconstrained nonlinear optimization. The improved initial state xa0 (“analysis”)
is obtained by minimizing the following cost function:
T
1
b
x0 − xb0 · B−1
(1a)
J (x0 ) =
0 · (x0 − x0 )
2
N
1X
T
(Hk (xk ) − yk ) · R−1
+
k · (Hk (xk ) − yk ) ,
2
k=1

xa0

=

arg min J (x0 )
x0

subject to xk = Mt0 →tk (x0 ) .

(1b)

Here M denotes the numerical model used to evolve the initial state vector x0
in time. Hk denotes the observation operator at assimilation time tk , and maps
the model state xk ≈ x(tk ) to the observation space. Rk is the observation
error covariance matrix. A practical assumption is that observations are independent, hence uncorrelated, which reduces Rk to a diagonal matrix containing
observation variances. The covariance matrices B0 and Rk are predefined, and
the quality of their approximation influences the resulting analysis.
The first term of the 4D-Var cost functional (1a) quantifies the mismatch
between the initial solution (x0 ) and the background state (xb0 ) at the initial
time (t0 ). This mismatch is computed in a least-squares sense, scaled by the
inverse background error covariance matrix B0 . The second term measures the
mismatch between the model trajectory (initialized from x0 ) and observations
yk taken at times tk , k = 0, . . . , N scaled by the inverse observation error
covariances Rk . Through B0 and Rk the inverse problem takes into account
the uncertainty in the data and model predictions. 4D-Var can be viewed as a
method of Bayesian inference, which computes the maximum likelihood initial
solution conditioned by the observations.
When assimilating observations only at the initial time t0 the method is
known as three-dimensional variational (3D-Var), as the additional time dimension is not present.
2

2.2. Computational aspects
The minimization problem (1b) is computed numerically using gradientbased iterative nonlinear solvers such as quasi-Newton, nonlinear conjugate gradients, and truncated Newton methods.
The iterative solver starts from an initial guess and advances to a minimizer,
along directions of descent computed using the gradients of J . The descent
directions depend on all three available sources of information; intermediate
iterations can be viewed as the result of partially assimilating information given
by observations.
The iterative nature of the solution approach has several important consequences. The resulting analysis is not the global minimizer, but rather is a local
one. When the computational cost per iteration is high, the available computational resources constrain the number of iterations that can be performed. Thus
in practice analyses are almost always suboptimal.
We will call the numerical model M (1b) the forward model (fwd) or forecast model. Applying gradient-based methods to minimize (1b) requires the
derivatives of J with respect to the initial model states. This can be achieved
using adjoint modeling [13, 14], a methodology that has been successfully implemented in optimization, sensitivity analysis and uncertainty quantification
[15, 21, 22]. The adjoint models (adj) can be based either on the linearization
of the original differential equations (“continuous” adjoint approach), or the linearization of the numerical method (“discrete” adjoint approach). A convenient
approach to generate discrete adjoint models is automatic differentiation [23],
which takes as input the source code of the fwd model, performs line by line
differentiation, and returns the source code of the adjoint models.
3. Observation Impact
3.1. 4D-Var sensitivity to observations
We now establish a relation between the analysis xa0 and the observations
yk . This is done using sensitivity analysis, and expresses how slight changes
in the observational data translate into changes in the resulting analysis. This
section follows the 4D-Var sensitivity approach of Daescu [11].
Consider the problem of finding a vector x = (x1 , x2 , ..., xn )T ∈ Rn that
minimizes the twice continuously differentiable cost function:
xa = arg min J (x, u) .
x

(2)

The function depends on the state vector x and parameter vector u ∈ Rm .
For any value ū of the parameter, the solution obeys the first order optimality
condition
∇x J (xa , ū) = 0 .
(3)
Assume that the Hessian of the function is positive definite at the optimum,
∇2x,x J (xa , ū) > 0. The implicit function theorem applied to (3) guarantees

3

there exists a vicinity of ū where the optimal solution is a smooth function of
the parameters, xa = xa (u), and has the sensitivity
−1
∇u xa (u) = −∇2u,x J (u, xa ) · ∇2x,x J (u, xa )
.
We apply this procedure to the 4D-Var cost function (1a) with the parameters
u being the observational data yk . The first-order necessary condition reads:
∇x0 J (xa0 )

=

B−1
0

a

x −x

b



+

N
X

MT0,k HTk R−1
k (Hk (xk ) − yk ) = 0 ,

(4)

k=1

where M0,k = M0t0 →tk (x0 ) is the tangent linear propagator associated with the
numerical model M, and Hk = Hk0 (xk ) is the linearized observation operator
at time tk . Differentiating (4) with respect to the observations yk
∇2yk ,x0 J (xa0 ) = −Rk Hk M0,k
leads to the following expression for the sensitivity of the analysis to observations:
 a T
∂x0
a
= R−1
(5)
∇yk x0 =
k Hk M0,k A0 ,
∂yk
−1
A0 = ∇2x0 ,x0 J (xa0 )
.
(6)
Consider a verification cost functional Ψ : Rn → R that measures the discrepancy between the analysis and a verification solution xv0 :
Ψ(xa0 ) =

1 a
(x − xv0 )T C (xa0 − xv0 ) .
2 0

(7)

Here the verification solution xv0 is also defined at t0 . The matrix C is a weighting matrix corresponding to a particular norm or restricting the verification to
a subdomain of the solution space. Using chain-rule differentiation and (5), the
sensitivity to observations (the gradient of Ψ with respect to observations yk )
is:
∇yk Ψ(xa0 ) = ∇yk xa0 · ∇xa0 Ψ(xa0 )
= R−1
k

Hk M0,k A0 C (xa0

(8)
−

xv0 ) .

The sensitivity (8) can be computed via the following steps:
• The sensitivity of the verification function with respect to the analysis is
∇xa0 Ψ(xa0 ) = C (xa0 − xv0 ) .
• A “supersensitivity” is obtained through solving a linear system with the
matrix the 4D-Var Hessian eqn:inverse-hessian
µ = A0 · ∇xa0 Ψ(xa0 ) .

(9)

This step dominates the computational cost of the sensitivity calculation.
4

• Finally, the verification sensitivity to observations yk is the supersensitivity vector propagated to time tk through the tangent linear model,
mapped through the observation selection operator, and scaled with the
inverse error covariance:
∇yk Ψ(xa0 ) = R−1
k Hk M0,k · µ .

(10)

The sensitivities of Ψ to the background xb0 and to the error covariance matrices
B and R can be derived in a similar fashion [11].
3.2. Observation impact matrix
3.2.1. Main features
Equation (5) defines a matrix whose elements are the sensitivities of each
component of the analysis vector xa0 to each component of the observation vector
yk assimilated at time tk . The observation impact matrix collects the sensitivities to all observations as follows:

  −1

R1 H1 M0,1
∇y1 xa0

 

H2 M0,2 
 a T  ∇y2 xa0   R−1

  2

∂x
 A0 .
=
=
(11)
T=
.. 
..




∂y
.
 .  


 

∇yN xa0
R−1
H
M
N
0,N
N
Each of the n columns of T represents the sensitivities of one particular
model state to all observations. For real applications n ∼ 107 − 1010 . Each row
of T contains the sensitivities of each state to one particular observation. In
typical data assimilation applications the number of observations (rows) is twothree orders of magnitude smaller than the number of model states (columns).
We now seek to understand the structure of the observation impact matrix
(11), which is the transpose of the sensitivity matrix ∂xa /∂y. The impact matrix is the product of submatrices of type R−1
k Hk M0,k , with the inverse of
the 4D-Var Hessian (6). The Hessian is symmetric and positive-definite when
evaluated at the minimum, but can lose positive-definiteness when evaluated at
an inexact analysis, such as when the minimization of 4D-Var was incomplete.
The inverse of the 4D-Var Hessian at the minimum approximates the analysis
error covariance [24]. This aposteriori error covariance characterizes the 4DVar process and quantifies the uncertainty reduction due to assimilating the
observations. Each column of this matrix represents the error covariance corresponding to a certain model state. The tangent linear model M0,k propagates
this error covariance to the observational time tk . The propagated perturbations
are mapped to observation space through Hk , and then scaled with R−1
k . The
inverse observation error covariance can be interpreted as a measure of trust assigned to each data point. Large values of the result are associated to points in
space and time where better measurements can benefit the assimilation process.

5

This metric can be assessed not only at times when we have observations available, but also at intermediate times, and can prove to be useful for deploying
adaptive observations.
Small changes in observations ∆y lead to a change in the analysis ∆xa0 which,
to first order, can be calculated using the observation impact matrix (11) as
∆xa0 = TT · ∆y .

(12)

We call the change in analysis ∆xa0 the observation impact. It is computed via
the following steps:
• Each observation change ∆yk is scaled and pulled back to the initial time
via the adjoint model:
∆yk = MT0,k HTk R−1
k · ∆yk .
• The aggregated contribution of all observation changes is
∆y =

N
X

∆yk ,

k=0

and in practice can be computed via a single adjoint run.
• The observation impact is obtained by solving a linear system whose matrix is the 4D-Var cost function Hessian
∆xa0 = A0 · ∆y .
3.2.2. Computational ingredients
For real models the matrices that appear in (11) are very large. The calculations of observation impact rely on matrix-vector products. The linear operator
M0,k is applied by running the tangent linear model from t0 to tk . Application
of the linearized observation mapping operator Hk and scaling by R−1
k can be
performed directly as they are relatively inexpensive.
The inverse Hessian has been successfully used in practical applications to
compute the supersensitivity. This was done by solving iteratively a linear
system for the supersensitivity, the system matrix being the 4D-Var Hessian.
Iterative solvers of Krylov-type require only Hessian-vector products. These
products can be computed by running the second-order adjoint model. Furthermore, the linear solvers can be preconditioned for faster convergence and several
matrix-free methods are readily available [25]. When the second-order adjoint
is not available, several approximations can be used as follows.
1. The finite difference of gradients
∇2x0 ,x0 J (xa0 ) · u ≈

∇x0 J (xa0 +  · u)T − ∇x0 J (xa0 )T
.


requires two first-order adjoint runs, which can be performed in parallel.
The accuracy of this approximation is typically low.
6

2. The Gauss-Newton approximation of the Hessian
∇2x0 ,x0 J (xa0 ) · u ≈ B−1
0 ·u+

N
X

MT0,k HTk R−1
k Hk M0,k · u

k=1

discards second order terms which contain the observation residual, so
this approximation is independent of the observation data. It requires one
tangent-linear model run followed by one first-order adjoint run.
3. Limited memory quasi-Newton approximations are based on the sequence
of solutions and gradients generated during the numerical optimization
procedure. An example of such an approximation is L-BFGS [26].
4. Efficient Computation of Observation Sensitivity and Impact
The computation of observation sensitivity or impact is a non-trivial task
since for practical problems of interest it is not feasible to build the full matrices.
We seek to develop computationally inexpensive approximations of the impact
matrix that capture the most important features, and whose accuracy can be
improved with an increased computational effort. Our approach is based on
the Singular Value Decomposition (SVD), a powerful tool to generate low-rank
approximations of large matrices. We present two algorithms, one iterative
(inherently serial), and one ensembled-based (inherently parallel).
Consider the SVD of (11) and the corresponding low-rank approximations
T = U S VT ,

T(p) = U(p) S(p) V(p) ,

where U and V are orthogonal, S is diagonal, U(p) and V(p) are the right and
the left singular vectors associated with the largest p singular values, and S(p)
has these dominant singular values on the diagonal formed.
T(p) has the smallest ”reconstruction error” in both Frobenius kT − T(p) kF
and L2 norms kT − T(p) k2 among all the rank p approximations of T. The
accuracy increases as more dominant singular modes are added. The cut-off
threshold is particular to the problem under study and can be determined from
the singular value spectrum decay.
4.1. An iterative (serial) approach for matrix-free low-rank approximations with
SVD
An ideal iterative SVD algorithm for our problem uses one matrix-vector
product per iteration, and reveals one new singular pair with each iteration.
The singular vectors are discovered in decreasing order of the magnitude of their
singular values. Thus, running the algorithm for p iterations would generate the
leading p singular values and their associated singular vectors.
There are no classic algorithms to compute iteratively the SVD of matrix
available only in operator form. We change our problem to computing the

7

leading eigenvectors of the product between the observation impact matrix and
its transpose, T∗ T = US2 U∗ , where
T∗ T
M0,k

= A∗0
=

N
X

∗

M0,k M0,k A0 =

k=1
R−1
k Hk

N
X

T∗k Tk ,

(13)

k=1

M0,k ,

Tk = M0,k A0 .

This problem is solved using Krylov-based approaches (Lanczos [27], Arnoldi
[28]), e.g., by the Jacobi-Davidson algorithm [29], available in the JDQZ software
library [30].
Our algorithm computes a low-rank approximation of A0 , the inverse 4D-Var
Hessian. An iterative procedure can be used to compute the smallest eigenvalues and the corresponding eigenvectors of the Hessian matrix, which are the
dominant eigenpairs of A0 . The low-rank approximation of the inverse Hessian
reads:


−1
∗
A0 = (V D V∗ ) = V D−1 V∗ ≈ V(p) D−1
V
.
(14)
(p)
(p)
We replace the inverse Hessian in (13) with its low-rank approximation:



 ∗
−1
∗
∗
M
M
V
D
V
T∗k Tk ≈
V(p) D−1
V
0,k
(p)
0,k
(p)
(p)
(p)
(p)

Wk

=

−1
∗
∗
V(p) D−1
(p) Wk Wk D(p) V(p) ,

=

M0,k V(p) .

(15)

The columns of Wk are the Hessian eigenvectors V(p) propagated forward
and scaled by the tangent linear model M0,k . The p tangent linear models can
be performed in parallel. This allows us to approximate (13) as:
T∗ T ≈
W

=

−1
∗
V(p) D−1
(p) W D(p) V(p) ,
N
X

Wk∗ Wk .

k=1

With the eigendecomposition:
−1
∗
D−1
(p) W D(p) = Vred Dred Vred ,

(16)

the approximation of (13) becomes

∗
T∗ T ≈ V(p) Vred Dred V(p) Vred .

(17)

This represents a truncated singular vector decomposition of T, with Dred the
matrix of dominant singular values, and Vp Vred the matrix of left singular
vectors.
Algorithm 1 summarizes the main computational steps. The computational
cost is dominated by the first step, where the expensive second-order adjoint
model is run repeatedly to generate the Hessian eigenpairs. The iterative approach is suited for applications that benefit from an iterative improvement of
the low-rank approximation. The methodology can be applied to any data assimilation system for which first and second order adjoint models are available.
8

Algorithm 1 Iterative algorithm for low-rank approximations
Solve iteratively the eigenvalue problem for the 4D-Var Hessian (14)
Map newly generated eigenvectors through the tangent linear model (15)
3: Compute the truncated SVD of the resulting matrix (16)
4: Project the left singular vectors onto the eigenvector base of the 4D-Var
Hessian and build the low-rank approximation of T (17)
1:
2:

4.2. An ensemble-based (parallel) approach for matrix-free low-rank approximations
This approach uses a “randomized SVD” algorithm [31] to compute the
Moore-Penrose pseudoinverse [32] of the Hessian. The 4D-Var Hessian matrix
A−1
0 is available only in operator form, i.e., only matrix vector products can be
evaluated by running the second order adjoint. The randomized algorithm is as
follows:
1: Draw p random vectors and form a matrix Ω.
−1
2: Compute the product Y = A0 Ω using Hessian-vector multiplications, i.e.,
running the second order order adjoint model for each column.
3: Construct the QR decomposition Y = QR.
Each of the above steps can be performed in parallel.
The columns of Q form an orthonormal basis for the range of Y. Randomized
SVD uses a series of algebraic manipulations, starting from the observation that
Q is also the orthonormal factor in the QR decomposition of A−1
0 :
A−1
0 = QB,

B = Q∗ A−1
0 ,

B∗ = A−1
0 Q.

(18)

Next, compute an SVD of B:
∗
B = UB ΣB VB

(19)

and replace (19) in (18) to obtain the SVD of A−1
0 :
A−1
0

∗
∗
= Q UB ΣB VB
= U A ΣB V B
.

(20)

The left singular vectors of A represent the projections of the left singular
vectors of B onto the columns of Q. The singular values and right singular
vectors of A are the same as those of B. The pseudoinverse of the 4D-Var
Hessian A−1
0 reads:
+
∗
A+
(21)
0 ≈ V B ΣB U A .
The observation impact matrix is approximated using the tangent linear
model M0,k and the pseudoinverse A+
0:
T≈

N
X

M0,k A+
0 .

k=1

The computational flow is summarized in Algorithm 2.
9

(22)

Algorithm 2 Sampling algorithm for low-rank approximations
Build the matrix B, through parallel second-adjoint runs (18)
Compute a full SVD of B (19)
−1
3: Project the left singular vectors of B in Q and form the SVD of A0 (20)
+
4: Compute the Hessian pseudoinverse A0 (21)
5: Build the impact matrix T through parallel tangent linear runs (22)

1:
2:

Computing the rows of B (18) is done as matrix-vector products through
second-order adjoint runs that can be performed in parallel. The last step,
which propagates the components of the pseudoinverse through the linearization
of the model, is achieved by multiple tangent linear model runs in parallel. The
tangent model results are checkpointed at each of the observation times, so that
only one run across the entire time horizon is necessary for each input vector.
5. Applications
We illustrate several applications of the observation impact matrix in 4DVar using the two-dimensional shallow water equations. We first describe the
system and its numerical discretization, then present the 4D-Var implementation and the experimental setting for data assimilation. Observation sensitivity
is computed both in full and using low-rank approximation to assess how well
the latter captures the essential features. We apply this analysis to three applications, namely, detecting the change in impact from perfect data to noisy
data, pruning the least important observations, and detecting faulty sensors.
5.1. Test problem: shallow water equations
The two-dimensional shallow-water equations (2D swe) [33] approximate
the movement of a thin layer of fluid inside a basin:
∂
∂
∂
h+
(uh) +
(vh) = 0
∂t
∂x
∂y


∂
∂
∂
1
(uh) +
(uvh) = 0
u2 h + gh2 +
∂t
∂x
2
∂y


∂
∂
∂
1 2
2
(vh) +
(uvh) +
v h + gh
= 0.
∂t
∂x
∂y
2

(23)

Here h(t, x, y) is the fluid layer thickness, and u(t, x, y) and v(t, x, y) are the
components of the velocity field of the fluid. The gravitational acceleration is
denoted by g.
We consider a spatial domain Ω = [−3, 3]2 (spatial units), and an integration window is t0 = 0 ≤ t ≤ tf = 0.1 (time units). Boundary conditions are
considered periodic. The space discretization is realized using a finite volume
scheme, and the time integration uses a fourth Runge-Kutta scheme, following

10

Table 1: Normalized CPU times of different sensitivity models. The forward
model takes one time unit to run.
fwd
tlm
foa
soa

1
2.5
3.7
12.8

fwd + tlm
fwd + foa
fwd + tlm + foa + soa

3.5
4.7
20

the method Lax-Wendroff [34]. The model uses a square -q × q uniform spatial discretization grid, which brings the number of model (state) variables to
n = 3 q2 .
We use the automatic differentiation tool TAMC [35, 36] to build various
sensitivity models, as follows. The tangent-linear model (tlm) propagates perturbations forward in time. The first-order adjoint model (foa) propagates
perturbations backwards in time, and efficiently computes the gradient of a
scalar cost functional defined on the model states. The second-order adjoint
model (soa) computes the product between the Hessian of the cost function
and a user-defined vector [22].
The overhead introduced by the sensitivity models is considerable. Table 1
presents the CPU times of tlm, foa, and soa models (normalized with respect
to that of one fwd run). One soa integration is about 3.5 times more expensive
than a single first-order adjoint run, while the foa takes 3.7 times longer than
the forward run. These relative costs depend on the discretization methodology
and the implementation of sensitivities. Our previous research showed how to
build efficient adjoint models by reusing computations performed in the forward
run [22]. For example, for the shallow water model, the alternative continuous
adjoints we built required a fraction of the forward model CPU time to run.
5.2. Data assimilation setup
The 4D-Var system is set up for a simple version of the “circular dam”
problem [37]. The reference initial height field h is a Gaussian bell of a width
equal to 1 length units centered on the grid, and the reference initial velocity
vector components are constant u = v = 0. The physical interpretation is that
the front of water falls to the ground (h decreases) under the effect of gravity
and creates concentric ripples which propagate towards the boundaries. Figures
1 represent snapshots of the reference trajectory at initial and final time.
The computational grid is square and regular with q = 40 grid points in
each direction, for a total of 4800 model variables (states). The simulation time
interval is set to 0.01 seconds, using N = 100 timesteps of size 0.0001 (time
units).
The h component of the a priori estimate (background) xb is generated by
adding a correlated perturbation to the h reference solution at initial time. The
background error covariance B0 corresponds to a standard deviation of 5% of
the reference field values. For the u and v components we use white noise to

11

40

130

40

130

35

125

35

125

30

120

30

120

25

115

25

115

20

110

20

110

15

105

15

105

10

100

10

100

5

95

5

95

0
0

90

0
0

10

20

30

40

(a) Initial time

10

20

30

40

90

(b) Final time

Figure 1: The height field h at the beginning and at the end of the reference
trajectory.
prescribe perturbations. The spatial error correlation uses a Gaussian decay
model, with a correlation distance of five grid points.
Synthetic observations are generated at the final time, by adding random
noise to the reference trajectory. Since the observation errors are assumed uncorrelated, the observation error covariance matrix R is diagonal. The standard
deviation for observation noise is 1% of the largest absolute value of the observations for each variable. The observation operator H is linear and selects observed
variables at specified grid points. For the following experiments, we consider
observations of all variables at each grid point.
The minimization of the 4D-Var cost function is performed with the L-BFGSB solver [26] using a fixed number of 100 iterations.
5.3. Experimental results
5.3.1. Validating the low-rank approximations of the observation impact matrix
In the first experiment, we study the effect of two possible sources of errors:
data noise and SVD truncation errors. These issues are inherent to performing
the data assimilation and the observation impact analysis.
We apply our computational methodology to the data assimilation scenario
introduced in Section 5.2. To assess how the reanalysis is affected by the presence of noise, we first assimilate perfect observations (i.e., reference model values), then assimilate the same observations with added small noise. Figure 2
reveals that the convergence of the numerical optimization procedure is similar
in the two scenarios, and leads to a similar decrease in the root mean square
(RMS) error for each variable.
We now apply the systematic approach presented in Section 4 to each one
of the two assimilation scenarios, in order to compute the sensitivity to observations (8) when the verification solution xv0 is the background xb0 and C is the
identity matrix:
∇yk Ψ(xa0 ) = T (xa0 − xb0 ) .
12

(24)

10
RMSE

RMSE
1

PERFECT DATA
NOISY DATA

1.3

PERFECT DATA
NOISY DATA

PERFECT DATA
NOISY DATA

1.3

10
RMSE

2

10

1.2

10

10

1.2

10
0

50
BFGS iterations

(a) h

100

0

50
BFGS iterations

(b) u

100

0

50
BFGS iterations

100

(c) v

Figure 2: Root mean square error decrease for each variable versus the number
of L-BFGS iterations.
This procedure provides the sensitivity of the 4D-Var increment norm xa0 − xb0
to data. Large sensitivity values will be associated to observations that had a
more important contribution to the 4D-Var correction.
Figure 3 plots the observation sensitivity (24) for each one of h, u and v, for
both perfect and noisy scenarios. The standard deviation for observation noise
is equal to 1% of the largest absolute value of the observations for each variable.
We notice the observation sensitivity is locally correlated. For the h observations
∇yk Ψ(xa0 ) is almost symmetric, while for u and v is aligned along the East-West
and South-North directions. Observation sensitivities for perfect and for noisy
data exhibit similar spatial features and numerical values. Reasonably small
noise in data does not significantly affect the sensitivity fields.
The results in Figure 3 are obtained with the full sensitivity matrix, built
offline. We compare these results with a low-rank approximation of the sensitivity matrix computed using Algorithm 1. The low rank approximation of the
sensitivity is

(∇yk Ψ(xa0 ))(p) = T(p) · xa0 − xb0 .
(25)
Figure 4 displays the singular value spectrum of the observation impact matrix
and the reconstruction error of the low-rank truncated values for observation
sensitivity. It can be noticed that the decay in singular value spectrum resembles
that of the truncation error. Based on this, a good choice for the rank of the
approximation is p ≈ 1600, which corresponds to one third of the full workload.
In Figure 5 we plot the low-rank approximation (25) of sensitivity to h
observations and the corresponding truncation error. A visual comparison with
the full-rank observation sensitivity plotted in Figure 3(a) reveals that the lowrank approximation captures well the main features of the sensitivity field.
5.3.2. Observation impact
The computational methodology developed in Section 4 is now applied to
compute the forward observation impact. Specifically, we map changes in the
observations (innovation vector ∆y = y − H(x)) to changes in the analysis
(∆xa0 ) using the relation (12).
13

40

40

5
4

30

30

4

20

2

y

y

3
20

2
1

10

10

0

0
0
0

20
x

0
0

40

(a) perfect h observations

20
x

40

(b) noisy h observations

40

40
1

0.5

0

20

−0.5

10

1

30

0.5
y

y

30

20

0

10

−0.5

−1
0
0

20
x

−1
0
0

40

(c) perfect u observations

40

(d) noisy u observations

40

40

30

1

30

1

20

0

20

0

y

y

20
x

10

10

−1

−1
0
0

20
x

0
0

40

(e) perfect v observations

20
x

40

(f) noisy v observations

Figure 3: Sensitivity field (24) for perfect and for noisy observations of h, u,
and v.

14

5

10
2

10

RMS

Value

0

10
0

10

−2

10

−5

10

−10

0

1600
3200
Mode

10

4800

(a) Singular value spectrum

0

1600
3200
Rank

4800

(b) Truncation error

Figure 4: Singular value decay for the observation impact matrix T (11) and
the corresponding truncation error norm.

40

40

5
4

30

0.4

30

0.2

20

0

10

−0.2

20

y

y

3
2
1

10

0
0
0

20
x

−0.4
0
0

40

(a) Low-rank estimate

20
x

40

(b) Truncation error field

Figure 5: Low-rank approximation of observation sensitivity (25) for h data and
the associated truncation error field for 1600 modes.

15

40

40

0.5

30

20

y

y

30
0

10
0
0

0.5

20

0

10
20
x

40

−0.5

0
0

−0.5
20
x

40

(a) Full-rank impact for center obser- (b) Low-rank impact for center observation
vation
40
0.05

30

0.05

30

20

0

20

10
0
0

y

y

40

−0.05
20
x

10
−0.05
0
0

40

0

20
x

40

(c) Full-rank impact for corner obser- (d) Low-rank impact for corner observation
vation

Figure 6: The impact of two h observations from assimilation time t1 00, placed
in the center and in the corner of the grid. Both full rank and reduced rank
solutions are shown. 500 modes are used for the reduced rank approximation.
We compute the observation impact in full-rank and low-rank corresponding
to two observations of h, one in the center (grid coordinates (20, 20)) and one in
the corner (at (5, 5)) at the locations represented with white markers. This is
achieved by multiplying TT with a vector containing just the innovation brought
by the observation whose impact we want to evaluate (all other vector entries
have value zero). The resulting impact fields are plot in Figure 6. The spatial
features of the observation impact have a local radial correlation in both cases.
This means the information carried by the observation is spread in its proximity
by the 4D-Var process; this is also true across the periodic boundaries of the
shallow water system. Moreover, the low-rank approximations are able to pick
up the important features of the full-rank calculations, and provide impacts of
a similar magnitude.
Having already computed the SVD of the observation impact matrix, we also
look at the directions in the data space ∆y that have the largest impact on the
analysis, and the directions in the 4D-Var correction space ∆x that benefit most
from the assimilation. These directions are given by the dominant left and right
singular vectors of TT , respectively. Figure 7 plots the first dominant left (and

16

8

40
1000

30

0

20

2

30
y

y

x 10

40

20

0

−1000
10

10

−2

−2000
0
0

20
x

0
0

40

20
x

40

(a) First dominant left singular vector (b) First 500 dominant left singular
(model space)
vectors (model space)
6

40

x 10

40

3
30

10

30

20

0

20

2

10
0
0

y

y

1
−1

10

−10
20
x

0
−2

0
0

40

20
x

40

(c) First dominant right singular vector (d) First 500 dominant right singular
(observation space)
vectors (observation space)

Figure 7: Principal directions of growth in observation and solution space, as
defined by the dominant left and right singular vectors of the observation impact
matrix.
right) singular vectors corresponding to the h variable, and also the composition
of the most important 500 directions for the observation and solution space using
the formula
vdomin =

500
X

s2i ∗ vi ,

(26)

i=1

where si and vi are the singular pair corresponding to the i-th mode.
5.3.3. Pruning observations based on sensitivity values
In this experiment we illustrate how sensitivity analysis can be used to select
the most useful data points.
For this, we compute the 4D-Var reanalysis and then apply the observation
impact methodology to compute the sensitivity of the cost function Ψ to each
one of the 4800 observations (24). We then split our observation set in three
subsets of 1600 observations, corresponding to each one of h, u and v. Within
each subset, we order the observations by their sensitivity; the top 800 observations form the high sensitivity set, and the ones in the bottom 800 form the
17

RMSE [log scale]

RMSE [log scale]

100

FULL
HIGH
LOW

10

1

1.3

10

FULL
HIGH
LOW

1.2

10

10
100
BFGS iterations [log scale]

1

10
100
BFGS iterations [log scale]

(a) h RMS error decrease versus the (b) u RMS error decrease versus the
number of L-BFGS iterations.
number of L-BFGS iterations.

30

1.3

10

y

RMSE [log scale]

40

1.2

10

FULL
HIGH
LOW

20
10
0
0

1
10
100
BFGS iterations [log scale]

10

20
x

30

40

(c) v RMS error decrease versus the (d) Location of high (red) and low imnumber of L-BFGS iterations.
pact observations

Figure 8: Location of high and low sensitivity observations, and their use in
data assimilation.

low sensitivity set. This procedure partitions our initial 4800 observations in
two halves with respect to observation sensitivity, taken variable-wise. The 4DVar data assimilation process is repeated using either one of the high or low
data sets.
Figure 8(a) shows the distribution of high versus low observations on our
computational grid. We run 4D-Var to assimilate only the high, followed by
assimilating only the low data points, and compare the decrease in the RMS
true error in h in Figure 8(b). The convergence results show that assimilating
observations of larger sensitivity yields slightly better results. When the number
of optimizer iterations increases, the two scenarios exhibit similar performance,
comparable to using all the observations.
5.3.4. Identifying faulty data using observation sensitivity
Finally, we consider the problem of detecting isolated faulty sensors and
try to solve it using our new computational capabilities. We prescribe certain
sensors to exhibit the wrong behavior and expect the observation sensitivity
analysis to reflect this.

18

We perform the 4D-Var scenario above again, using perfect observations at
each grid point, with the exception of two points where the observation values
are increased by a factor of 10 (thus simulating two faulty instruments). The two
points are located at coordinates (20, 20), the grid center where fluid height field
has large values initially, and at coordinates (10, 10), which is initially outside
the water height profile.
The 4D-Var data assimilation uses the faulty data together with the correct
values. Figure 9(a) plots the correction in the variable h, i.e., the quantity
(xa0 − xb0 ). The correction field is smooth and there is no indication that there
are very large errors in some data points. Figure 9(b) plots the supersensitivity
vector µ obtained by mapping the increment through the inverse 4D-Var Hessian
(9). The supersensitivity field exhibits a clear structure consisting of large
values in concentric circles around the two locations with faulty observations.
Figure 9(c) shows the observation sensitivity y, obtained by propagating the
supersensitivity through the tangent linear model (10). The observation impact
methodology clearly identifies two points of abnormally large sensitivity. These
are exactly the locations of the faulty sensors, and because their values are very
different than the surrounding field, they have an abnormally large impact on
the analysis. This is precisely identified by the observation impact metrics.
6. Conclusions and future work
This research is motivated by dynamic data driven applications which use
measurements of the real system state to constrain computer model predictions.
We develop efficient approaches to compute the observation impact in 4D-Var
data assimilation, i.e., the contribution of each data point to the final analysis. Quantification of observation impact is important for applications such
as data pruning, tuning the parameters of data assimilation, identification of
faulty sensors, and optimal configuration of sensor networks. Immediate applications include numerical weather predicition, climate and air-quality forecast,
hydrology, renewable energy systems and biohazard proliferation.
While 4D-Var is the state of the art data assimilation methodology, and is
widely used in many operational centers, it is one of the most computationally
challenging approaches to data assimilation. The computation of observation
impact adds to the cost of performing analyses. We review the derivation of
the sensitivity equations that lead to the observation impact matrix, and make
in-depth comments about its structure, functionality and computation requirements. One contribution of this work is to develop highly efficient implementation strategies for the 4D-Var sensitivity equations via a smart use of adjoint
models. Two matrix free algorithms are proposed, one serial and one parallel, to
compute approximations of observation impact matrix. These algorithms compute SVD-based low-rank approximations, which capture the most important
features of observation impact matrix. The accuracy of the generated approximations scales with the computational cost.
A second contribution of this work is to illustrate the use of observation impact in solving two practical applications: pruning observations of small impact,
19

40

40
10
30

5

20

0

10

−5

y

y

30

10

20

0

10
−10

0
0

−10
20
x

0
0

40

(a) 4D-Var increment

20
x

40

(b) Supersensitivity field

40
35

0

30
−10
25
−20

20
15

−30

10
−40
5
0
0

−50
10

20

30

40

(c) Sensitivity to observations

Figure 9: Observation sensitivity field when the assimilated data is corrupted at
two locations with coordinates (10,10) and (20,20). The location of the faulty
sensors is unknown to the data assimilation system, but is retrieved via the
observation impact methodology.

20

and detecting faulty sensors. Numerical experiments that validate the proposed
computational framework are carried out with a two dimensional shallow water
equations system.
Several future research directions emerge from this study. On the computational side, the impact approximation algorithms can be further developed to
achieve superior performance. A rigorous analysis of approximation errors will
be developed to guide the choice of the number of iterations and the truncation
level. On the application side, our will be integrated in real large scale applications to provide a measure of importance of different measurements in real time
assimilation. In hindsight, the impact can be used to prune the data subsets and
to detect erroneous data points. In foresight, it can be used to design efficient
strategies of sensor placement for targeted and adaptive observations.
Acknowledgements
This work was supported by the National Science Foundation through the
awards NSF DMS-0915047, NSF CCF-0635194, NSF CCF-0916493 and NSF
OCI-0904397, and by AFOSR DDDAS program through the awards FA9550–
12–1–0293–DEF and AFOSR 12-2640-06.
References
[1] T. Palmer, T. Gelaro, J. Barkmeijer, and R. Buizza. Singular vectors,
metrics, and adaptive observations. Journal of the Atmospheric Sciences,
55(4):633–653, 1998.
[2] W. Liao, A. Sandu, and G. Carmichael. Total energy singular vector analysis for atmospheric chemical transport models. Monthly Weather Review,
134(9):2443–2465, 2006.
[3] D. Zupanski, A.Y. Hou, S.Q. Zhang, M. Zupanski, C.D. Kummerow, and
S.H. Cheung. Applications of information theory in ensemble data assimilation. Quarterly Journal of the Royal Meteorological Society, 133(627):1533–
1545, 2007.
[4] K. Singh, A. Sandu, and M. Jardak. Information theoretic metrics to
characterize observations in variational data assimilation. Proceedings of
the International Conference on Computational Science, 9(1):1047–1055,
2012.
[5] L.M. Berliner, Z.Q. Lu, and C. Snyder. Statistical design for adaptive
weather observations. Journal of the Atmospheric Sciences, 56(15):2536–
2552, 1999.
[6] W. Kang and L. Xu. Optimal placement of mobile sensors for data assimilations. Tellus A, 64(0), 2012.

21

[7] R. Gelaro, Y. Zhu, and R.M. Errico. Examination of various-order adjointbased approximations of observation impact. Meteorologische Zeitschrift,
16(6):685–692, 2007.
[8] R. Gelaro and Y. Zhu. Examination of observation impacts derived from
observing system experiments (OSEs) and adjoint models. Tellus A,
61(2):179–193, 2009.
[9] Y. Tremolet. Computation of observation sensitivity and observation impact in incremental variational data assimilation. Tellus A, 60(5):964–978,
2008.
[10] F.X. Le Dimet and H.E. Ngodock. Sensitivity analysis in variational data
assimilation. Journal of the Meteorological Society of Japan, 75(1):245–255,
1997.
[11] D.N. Daescu. On the sensitivity equations of four-dimensional variational
(4D-Var) data assimilation. Monthly Weather Review, 136(8):3050–3065,
2008.
[12] D.N. Daescu and R. Todling. Adjoint sensitivity of the model forecast to
data assimilation system error covariance parameters. Quarterly Journal
of the Royal Meteorological Society, 136(653):2000–2012, 2010.
[13] D.G. Cacuci. Sensitivity theory for nonlinear systems. I: Nonlinear functional analysis approach. Journal of Mathematical Physics, 22:2794, 1981.
[14] Z. Wang, I.M. Navon, F.X. Le Dimet, and X. Zou. The second order adjoint
analysis: theory and applications. Meteorology and Atmospheric Physics,
50(3):3–20, 1992.
[15] A. Sandu and L. Zhang. Discrete second order adjoints in atmospheric chemical transport modeling. Journal of Computational Physics,
227(12):5949–5983, 2008.
[16] J. Ye. Generalized low rank approximations of matrices. Machine Learning,
61(1-3):167–191, 2005.
[17] M.W. Berry and R.D. Fierro. Low-rank orthogonal decompositions for
information retrieval applications. Numerical Linear Algebra with Applications, 3(4):301–327, 1996.
[18] I.S. Dhillon and D.S. Modha. Concept decompositions for large sparse text
data using clustering. Machine Learning, 42(1-2):143–175, 2001.
[19] R. Daley. Atmospheric data analysis. Cambridge University Press, Cambridge, 1991.
[20] E. Kalnay. Atmospheric modeling, data assimilation and predictability.
Cambridge University Press, Cambridge, 2002.

22

[21] A. Sandu, D.N. Daescu, G.R. Carmichael, and T. Chai. Adjoint sensitivity
analysis of regional air quality models. Journal of Computational Physics,
204(1):222–252, 2005.
[22] A. Cioaca, A. Sandu, and M. Alexe. Second-order adjoints for solving PDEconstrained optimization problems. Optimization Methods and Software,
27(4-5):625–653, 2011.
[23] A. Griewank. On automatic differentiation. Mathematical Programming:
recent developments and applications, 6:83–107, 1989-.
[24] I.Y. Gejadze, F.X. Le Dimet, and V. Shutyaev. On analysis error covariances in variational data assimilation. SIAM Journal on Scientific Computing, 30(4):1847–1874, 2008.
[25] A. Cioaca, A. Sandu, and E. de Sturler. Efficient methods for computing
observation impact in 4d-var data assimilation. Technical report, 2012.
[26] C. Zhu, R.H. Byrd, P. Lu, and J. Nocedal. Algorithm 778: L-BFGS-B:
Fortran subroutines for large-scale bound-constrained optimization. ACM
Transactions on Mathematical Software (TOMS), 23(4):550–560, 1997.
[27] C. Lanczos. An iteration method for the solution of the eigenvalue problem
of linear differential and integral operators. United States Government
Press Office, 1950.
[28] W.E. Arnoldi. The principle of minimized iterations in the solution of the
matrix eigenvalue problem. Quarterly of Applied Mathematics, 9(1):17–29,
1951.
[29] G.L.G. Sleijpen and H.A. Van der Vorst. A Jacobi–Davidson iteration
method for linear eigenvalue problems. SIAM Review, 42(2):267–293, 2000.
[30] D.R. Fokkema and M.B. van Gijzen. Short manual for the JDQZ-package,
1999.
[31] E. Liberty, F. Woolfe, P.G. Martinsson, V. Rokhlin, and M. Tygert. Randomized algorithms for the low-rank approximation of matrices. Proceedings of the National Academy of Sciences, 104(51):20167–20172, 2007.
[32] T.O. Lewis and T.G. Newman. Pseudoinverss of postive semidefeinite matrices. SIAM Journal on Applied Mathematics, 16(4):701–703, 1968.
[33] D.D. Houghton and A. Kasahara. Nonlinear shallow fluid flow over an
isolated ridge. Communications on Pure and Applied Mathematics, 21(1):1–
23, 1968.
[34] R. Liska and B. Wendroff. Composite schemes for conservation laws. SIAM
Journal of Numerical Analysis, 35(6):2250–2271, 1998.

23

[35] R. Giering. Tangent linear and adjoint model compiler. User Manual,
TAMC Version, 4, 1997.
[36] R. Giering and T. Kaminski. Recipes for adjoint code construction. ACM
Transactions on Mathematical Software, 24(4):437–474, 1998.
[37] K. Anastasiou and C.T. Chan. Solution of the 2d shallow water equations using the finite volume method on unstructured triangular meshes.
International Journal for Numerical Methods in Fluids, 24(11):1225–1245,
1997.

24

