Ensemble Kalman Methods With Constraints

TL;DR

This paper develops a general framework for incorporating equality and inequality constraints into Ensemble Kalman methods, enhancing their applicability in state and parameter estimation problems with prior information.

Contribution

It introduces a novel methodology and theoretical justification for constrained Ensemble Kalman methods, supported by numerical experiments.

Findings

Effective enforcement of constraints in Ensemble Kalman methods

Theoretical validation of the constrained methodology

Numerical experiments demonstrating improved estimation accuracy

Abstract

Ensemble Kalman methods constitute an increasingly important tool in both state and parameter estimation problems. Their popularity stems from the derivative-free nature of the methodology which may be readily applied when computer code is available for the underlying state-space dynamics (for state estimation) or for the parameter-to-observable map (for parameter estimation). There are many applications in which it is desirable to enforce prior information in the form of equality or inequality constraints on the state or parameter. This paper establishes a general framework for doing so, describing a widely applicable methodology, a theory which justifies the methodology, and a set of numerical experiments exemplifying it.

Figures

Equations

$$\text{Dynamics Model:}{v}_{j+1}$$

$$\text{Data Model:}{y}_{j+1}$$

$$\text{Probabilistic Structure:}{v}_0$$

$${\hat{v}}_{j+1}^{(n)}$$

$${\hat{m}}_{j+1}$$

$${\hat{C}}_{j+1}$$

$$v_{j+1}^{(n)}=\underset{v}{\mathrm{argmin}}{I}_{\mathrm{filter},\mathrm{j},\mathrm{n}}(v)$$

$$I_{\mathrm{filter},j,n}(v):=$$

$$\displaystyle\mathopen{}\mathclose{{}\left\{\begin{array}[]{ll}\frac{1}{2}v^{T}\Bigr{(}H^{T}{\Gamma}^{-1}H+\widehat{C}_{j+1}^{-1}\Bigr{)}v-\Bigr{(}\widehat{C}_{j+1}^{-1^{T}}\widehat{v}_{j+1}^{(n)}+H^{T}{\Gamma^{-1}}^{T}y_{j+1}^{(n)}\Bigr{)}^{T}v&\mbox{if }v-\widehat{v}_{j+1}^{(n)}\in\mathcal{R}(\widehat{C}_{j+1}).\\ \infty&\mbox{otherwise.}\end{array}}\right.$$

$$S_{j+1}$$

$$K_{j+1}$$

$$y_{j+1}^{(n)}$$

$$v_{j+1}^{(n)}$$

$$v_{j+1}^{(n)}=\widehat{v}_{j+1}^{(n)}+\underset{(a,v^{\prime})\in\mathcal{A}}{\mathrm{argmin}}\,\Bigr{(}\frac{1}{2}|y^{\prime}-Hv^{\prime}|^{2}_{\Gamma}+\frac{1}{2}\langle a,v^{\prime}\rangle\Bigr{)}$$

$$v_{\epsilon}=\underset{v\in\mathcal{H}_{1}}{\mathrm{argmin}}\,\Bigl{(}\frac{1}{2}|y_{j+1}^{(n)}-Hv|^{2}_{\Gamma}+\frac{1}{2}|\widehat{v}_{j+1}^{(n)}-v|^{2}_{\widehat{C}_{\epsilon}}\Bigr{)}$$

$$\mathcal{L}(v^{\prime },a,\lambda )=\frac{1}{2}|y^{\prime }-H{v}^{\prime }{|}_{\Gamma }^2+\frac{1}{2}⟨a,v^{\prime }⟩+⟨\lambda ,\hat{C}a-v^{\prime }⟩$$

$$-H^T{\Gamma }^{-1}(y^{\prime }-H{v}^{\prime })+\frac{1}{2}a-\lambda$$

$$\frac{1}{2}{v}^{\prime }+\hat{C}\lambda$$

$$v^{\prime }-\hat{C}a$$

$$-H^T{\Gamma }^{-1}(y^{\prime }-H\hat{C}a)+a=0.$$

$$a=(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}{H}^T{\Gamma }^{-1}{y}^{\prime }.$$

$$v_{j+1}^{(n)}$$

$$={\hat{v}}_{j+1}^{(n)}+\hat{C}a$$

$$={\hat{v}}_{j+1}^{(n)}+\hat{C}(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}{H}^T{\Gamma }^{-1}{y}^{\prime }$$

$$={\hat{v}}_{j+1}^{(n)}+\hat{C}(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}{H}^T{\Gamma }^{-1}(y_{j+1}^{(n)}-H{\hat{v}}_{j+1}^{(n)}).$$

$$K=\hat{C}(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}{H}^T{\Gamma }^{-1}$$

$$v_{j+1}^{(n)}=(I-KH){\hat{v}}_{j+1}^{(n)}+K{y}_{j+1}^{(n)}.$$

$$(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}{H}^T{\Gamma }^{-1}=H^T(H\hat{C}{H}^T+\Gamma {)}^{-1}.$$

$$(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}(H^T{\Gamma }^{-1}H\hat{C}{H}^T+H^T)=H^T$$

$$(H^T{\Gamma }^{-1}H\hat{C}+I{)}^{-1}(H^T{\Gamma }^{-1})(H\hat{C}{H}^T+\Gamma )=H^T.$$

$$(H^T{\Gamma }^{-1}H+{\hat{C}}_{ϵ}^{-1})v_{ϵ}=H^T{\Gamma }^{-1}{y}_{j+1}^{(n)}+{\hat{C}}_{ϵ}^{-1}{\hat{v}}_{j+1}^{(n)}.$$

$$v_{ϵ}=({\hat{C}}_{ϵ}-{\hat{C}}_{ϵ}{H}^T(H{\hat{C}}_{ϵ}{H}^T+\Gamma {)}^{-1}H{\hat{C}}_{ϵ})(H^T{\Gamma }^{-1}{y}_{j+1}^{(n)}+{\hat{C}}_{ϵ}^{-1}{\hat{v}}_{j+1}^{(n)}).$$

$$({\hat{C}}_{ϵ}-{\hat{C}}_{ϵ}{H}^T(H{\hat{C}}_{ϵ}{H}^T+\Gamma {)}^{-1}H{\hat{C}}_{ϵ}){\hat{C}}_{ϵ}^{-1}=(I-{\hat{C}}_{ϵ}{H}^T(H{\hat{C}}_{ϵ}{H}^T+\Gamma {)}^{-1}H)$$

$$({\hat{C}}_{ϵ}-$$

$$={\hat{C}}_{ϵ}{H}^T(I-(H{\hat{C}}_{ϵ}{H}^T+\Gamma {)}^{-1}H{\hat{C}}_{ϵ}{H}^T){\Gamma }^{-1}$$

$$={\hat{C}}_{ϵ}{H}^T(H{\hat{C}}_{ϵ}{H}^T+\Gamma {)}^{-1}\Gamma {\Gamma }^{-1}$$

$$={\hat{C}}_{ϵ}{H}^T(H{\hat{C}}_{ϵ}{H}^T+\Gamma {)}^{-1}$$

Full text

††: Inverse Problems

, , , ,

Keywords: ensemble Kalman methods, equality and inequality constraints, derivative-free optimization, convex optimization

1 Introduction

1.1 Overview

Kalman filter based methods have been enormously successful in both state and parameter estimation problems. However, a major disadvantage of such methods is that they do not naturally take constraints into account. The ability to constrain a system often has a number of advantages that can play an important role in state and parameter estimation: they can be used to enforce physicality of modeled systems (non-negativity of physical quantities, for example); relatedly they can be used to ensure that computational models are employed only within state and parameter regimes where the model is well-posed; and finally the application of constraints may provide robustness to outlier data. Resulting improvements in algorithmic efficiency and performance, by means of enforcing constraints, has been demonstrated in the recent literature in a diverse set of fields, including process control [1], biomechanics [2], cell energy metabolism [3], medical imaging [4], engine health estimation [5], weather forecasting [6], chemical engineering [7], and hydrology [8]. Within the Kalman filtering literature the need to incorporate constraints is widely recognized and has been addressed in a systematic fashion by viewing Kalman filtering from the perspective of optimization. Indeed this optimization perspective leads naturally to many extensions, and to the incorporation of constraints in particular. Including constraints in Kalman filtering, via optimization, lends itself to an elegant mathematical framework, to a practical computational framework, and has potential in numerous applications. Surveys of the work may be found in the papers of Aravkin, Burke and co-workers [9, 10] and our work in this paper may be viewed as generalizing their perspective to the ensemble setting.

In the probabilistic view of filtering methods, constraints may be introduced by moving beyond the Gaussian assumptions that underpin Kalman methods and imposing constraints through the prior distributions on states and/or parameters. This, however, can create significant computational burden as the resulting distributions cannot be represented in closed form, through a finite number of parameters, in the way that Gaussian distributions can be. Here we circumvent this issue by taking the viewpoint that ensemble Kalman methods constitute a form of derivative-free optimization methodology, eschewing the probabilistic interpretation. The ensemble is used to calculate surrogates for derivatives. With this optimization perspective, constraints may be included in a natural way. Standard ensemble Kalman methods employ a quadratic optimization problem encapsulating relative strengths of belief in the predictions of the model and the data; these optimization problems have explicit analytic solutions. To impose constraints the optimization problem is solved only within the constraint set; when the constraints form a non-empty closed convex set, this constrained optimization problem has a unique solution.

In this introductory section, we give a literature review of existing work in this setting, we describe the contributions in this paper, and we outline notation used throughout.

1.2 Literature Review

Overviews of state estimation using Kalman based methods may be found in [11, 12, 13, 14]. The focus of this article is on ensemble based Kalman methods, introduced by Evensen in [15] and further developed in [16, 11]. The extension of the ensemble Kalman methodology to parameter estimation and inverse problems is overviewed in [17], especially for oil reservoir applications, and in an application-neutral formulation in [18]. Equipping Kalman-based methods with constraints can be desirable for a variety of inter-linked reasons described in the previous subsection: to enforce known physical boundaries in order to improve estimation accuracy; to operationalize filtering of a model which is ill-posed in subsets of its state or parameter space; and to provide robustness to noisy data and outlier events.

In extending the Kalman filter to non-Gaussian settings, a number of methods may be considered. Particle filters provide the natural methodology if propagation of probability distributions is required for state [19] or parameter [20] estimation. In the optimization setting, there are three primary methodologies: the extended Kalman filter, the unscented Kalman filter and the ensemble Kalman filter. The extended Kalman filter is based on linearization of the nonlinear system and therefore needs the computation of derivatives for propagation of the state covariance; this makes them unattractive in high dimensional problems. Unscented and ensemble Kalman filters, on the other hand, can be considered as particle-based methods which are derivative-free. In the unscented Kalman filter, the particles (sigma points) are chosen deterministically and are propagated through the nonlinear system to approximate the covariance, which is then corrected using the Kalman gain to compute the new sigma points. In the ensemble Kalman filter, the particles (ensemble members) are chosen randomly from the initial ensemble and are propagated through the dynamical system and corrected using the Kalman gain without needing to maintain the covariance.

In [21], and more recently in [22], overviews of different ways to impose constraints in linear and nonlinear state estimation are presented. To ensure that the estimates satisfy the constraints, moving horizon based estimators that solve a constrained optimization problem have been proposed [23, 24]. The paper [25] proposed a recursive nonlinear dynamic data reconciliation (RNDDR) approach based on extended Kalman filtering to ensure that state and parameter estimates satisfy the imposed bounds and constraints. The updated state estimates in this method are obtained by solving an optimization problem instead of using the Kalman gain. The resulting covariance calculations are, however, still similar to the Kalman filter: that is, unconstrained propagation and correction involving the Kalman gain, which can affect the accuracy of the estimates. To eliminate this deficiency, [26] proposed a Kullback-Leibler based method to update states and error covariances by solving a convex optimization problem involving conic constraints.

On the other hand, the paper [27] combined the concept of the unscented transformation [28] with the RNDDR formulation. In the prediction step, they propose step sizes to scale sigma points asymmetrically to better approximate the covariance information in the presence of lower and upper bounds. Then, for the update of each sigma point, they solve a constrained optimization problem. One disadvantage of this procedure is that the chosen step sizes for scaling the sigma points can only ensure the bound constraints. The paper [1] also tested various algorithms based on constrained optimization, projection [29] and truncation [5] to enforce bound constraints on unscented Kalman filtering. The paper [30] developed a class of estimators named constrained unscented recursive estimators to address the limitations of the unscented RNDDR method using optimization-based projection algorithms for obtaining sigma points in the presence of convex, non-convex and bound constraints.

As mentioned earlier, since the corrected covariance is used to compute the sigma points, unscented formulations always require enforcing constraints in both propagation and correction/update steps. In contrast, ensemble-based methods only require constraints to be enforced in the update step. In this context, the paper [8] tested projection and accept/reject methods to constrain ensemble members in a post-processing step, after application of the unconstrained ensemble Kalman filter. In the former, they project the updated ensemble members to the feasible space if they violate the constraints and in the latter they enforce the updated ensemble members to obey the constraints by resampling the dynamic and/or data model errors. On the other hand, [31, 32] proposed updating the state estimates in ensemble Kalman filtering by solving a constrained optimization problem while truncating the Gaussian distribution of the initial ensemble. The paper [6] demonstrated how to enforce a physics-based conservation law on an ensemble Kalman filtering based state estimation problem by formulating the filter update as a set of quadratic programming problems arising from a linear data acquisition model subject to linear constraints. Here we develop this body of work on constraining ensemble Kalman techniques, providing a unifying framework with an underpinning theoretical basis.

1.3 Our Contribution

The preceding literature review demonstrates that the imposition of constraints on state and parameter estimation procedures is highly desirable. It also indicates that ensemble Kalman methods offer the most natural context in which to attempt to do this, as extended Kalman methods do not scale well to high dimensional state or parameter space, whilst the unscented filter does not lend itself as naturally to the incorporation of constraints.

In this paper we build on the application-specific papers [8, 6] which demonstrate how to impose a number of particular constraints on ensemble based parameter and state estimation problems respectively. We formulate a very general methodology which is application-neutral and widely applicable, thereby making the ideas in [8, 6] accessible to a wide community of researchers working in inverse problems and state estimation. We also describe a straightforward mathematical analysis which demonstrates that the resulting algorithms are well-defined since they involve the solution of quadratic minimization problems subject to convex constraints at each step of the algorithm; these optimization problems have a unique solution. And finally we showcase the methodology on two applications, one from biomedicine and one from seismology. All of the algorithms discussed are clearly stated in pseudo-code.

Section 2 outlines the ensemble Kalman (EnKF) methodology for state estimation, with and without constraints. In section 3 the same program is carried out for ensemble Kalman inversion (EKI). Section 4 describes the numerical experiments which illustrate the foregoing ideas.

1.4 Notation

Throughout the paper we use $\mathbb{N}$ to denote the positive integers $\{1,2,3,\cdots\}$ and $\mathbb{Z}^{+}$ to denote the non-negative integers $\mathbb{N}\cup\{0\}=\{0,1,2,3,\cdots\}.$ The matrix $I_{M}$ denotes the identity on $\mathbb{R}^{M}.$ We use $|\cdot|$ to denote the Euclidean norm, and the corresponding inner-product is denoted $\langle\cdot,\cdot\rangle.$ A symmetric, square matrix $A$ is positive definite (resp. positive semi-definite) if the quadratic form $\langle u,Au\rangle$ is positive (resp. non-negative) for all $u\neq 0$. By $|\cdot|_{B}$ we denote the weighted norm defined by $|v|_{B}^{2}=v^{*}B^{-1}v$ for any positive-definite $B$. The corresponding weighted Euclidean inner-product is given by $\langle\cdot,\cdot\rangle_{B}:=\langle\cdot,B^{-1}\cdot\rangle.$ We use $\otimes$ to denote the outer product between two vectors: $(a\otimes b)c=\langle b,c\rangle a.$

2 Ensemble Kalman State Estimation

2.1 Filtering Problem

Consider the discrete-time dynamical system with noisy state transitions and noisy observations in the form:

[TABLE]

We assume that $\mathcal{H}_{1},\mathcal{H}_{2}$ are finite dimensional Hilbert spaces. Then $v_{j}\in\mathcal{H}_{1}$, and $\Psi:\mathcal{H}_{1}\mapsto\mathcal{H}_{1}$ is the state-transition operator. The operator $H:\mathcal{H}_{1}\mapsto\mathcal{H}_{2}$ is the linear observation operator and $y_{j}\in\mathcal{H}_{2}.$ The covariance operators $C_{0},\Sigma$ are assumed to be invertible. The objective of filtering is to estimate the state $v_{j}$ of the dynamical systems at time $j$, given the data $\{y_{\ell}\}_{\ell=1}^{j}.$

Remark 2.1.

• •

We may extend the methodology in this paper to the setting where $\mathcal{H}_{1},\mathcal{H}_{2}$ are separable infinite dimensional Hilbert spaces. The covariance operators $C_{0},\Sigma$ are assumed trace-class on $\mathcal{H}_{1}$, and $\Gamma$ on $\mathcal{H}_{2}$ to ensure that the initial condition $v_{0}$ and the noises $\xi_{j}$ and $\eta_{j}$ live in $\mathcal{H}_{1},\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ (respectively) with probability one. The update formulae we derive require operator composition and inversion, together with minimization of quadratic functionals on $\mathcal{H}_{1}$ subject to convex constraints. Provided all of these operations can be carried out, then the methods derived here are well-defined in the general Hilbert space setting. This fact is important because it means that the methods derived have a robustness to mesh refinement and similar procedures arising when the problem of interest is specified via a partial differential equation, or other infinite dimensional problem.

• •

We restrict attention to linear observation operators $H$ because this leads to solvable quadratic optimization problems within the context of Kalman-based methods. In principle, a non-linear observation operator could be used, but the optimization problems defining the algorithms arising in this paper might not have a unique solution in this setting.

2.2 Ensemble Kalman Filter

The ensemble Kalman filter is a particle-based sequential optimization approach to the state estimation problem. The particles are denoted by $\{{v}_{j}^{(n)}\}^{N}_{n=1}$ and represent a collection of $N$ candidate state estimates at time $j$. The method proceeds as follows. The state of all the particles at time $j+1$ are predicted using the dynamics model to give $\{\widehat{v}_{j+1}^{(n)}\}^{N}_{n=1}$. The resulting empirical covariance of the particles is then used to define the objective function $I_{{\rm filter},j,n}(v)$, which encapsulates the model-data compromise. This is minimized in order to obtain the updates $\{{v}_{j+1}^{(n)}\}^{N}_{n=1}.$ To understand the origin of this optimization perspective on ensemble Kalman methods we argue as follows. In equation (4.10) of [33], it is shown that the data incorporation step of the Kalman filter may be written as a quadratic optimization problem for the state. In equation (4.15) of [33], the ensemble Kalman filter is written by using this quadratic minimization principle with an empirically (from the ensemble) computed covariance.

The prediction step is

[TABLE]

Here we have $\xi_{j}^{(n)}\sim N(0,\Sigma)$ i.i.d.. Because the empirical covariance contains only $N-1$ independent pieces of information, (1c) is sometimes scaled by $N-1$ and not $N$; making this change would lead to no changes in the statements and proofs of all the theorems, and would only affect the definition of covariance within the algorithms.

Let $\mathcal{R}(\widehat{C}_{j+1})$ denote the range of $\widehat{C}_{j+1}$. The update step is then

[TABLE]

where

[TABLE]

It can be useful to rewrite the objective function for the optimization problem in an equivalent and more standard form for input to software:

[TABLE]

The $y_{j+1}^{(n)}$ are either identical to the data $y_{j+1}$, or found by perturbing it randomly.

Note that $\widehat{C}_{j+1}$ is an operator of rank at most $N-1$, and thus can only be invertible when $N-1$ is larger than the dimension of $\mathcal{H}_{1}$. For moderate- and high-dimensional systems, it is often impractical to satisfy this condition. However, the minimizing solution can be found by regularizing $\widehat{C}_{j+1}$ by addition of $\epsilon I$ for $\epsilon>0$, deriving the update equations and then letting $\epsilon\to 0.$ We give the resulting formulae, and then justify them immediately afterwards, in the following subsubsection. Alternatively it is possible to directly seek a solution in $\mathcal{R}(\widehat{C}_{j+1})$, which is a subspace of dimension $N-1$; this is done in the subsequent subsubsection.

2.2.1 Formulation In The Original Variables

The well-known Kalman update formulae arising from solution of the minimization problem (2), (5) are as follows:

[TABLE]

Here $\eta_{j}^{(n)}\sim N(0,\Gamma)$ i.i.d. and the constant $s$ takes value [math] or $1$. When $s=1$ the $y_{j+1}^{(n)}$ are referred to as perturbed observations. The choice $s=1$ is made to ensure the correct statistics of the updates in the linear Gaussian setting when a probabilistic viewpoint is taken, and more generally to introduce diversity into the ensemble procedure when an optimization viewpoint is taken. Derivation of the formulae may be found in [33]. In brief the formulae arise from completing the square in the objective function $I_{{\rm filter},j,n}(\cdot)$ and then applying the Sherman–Morrison formula to rewrite the updates in the data space rather than state space; the latter is advantageous in many applications where $\mathcal{H}_{2}$ has dimension much smaller than $\mathcal{H}_{1}.$

We summarize with the following pseudo-code:

An equivalent formulation of the minimization problem is now given by means of a penalized Lagrangian approach to incorporate the property that the solution of the optimization problem lies in the range of the empirical covariance. The perspective is particularly useful when further constraints are imposed on the solution of the optimization problem.

Theorem 2.2**.**

Suppose that the dimensions of $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are finite. Let $j$ be in $\mathbb{Z}^{+}$ and $1\leq n\leq N$. Define $y^{\prime}=y_{j+1}^{(n)}-H\widehat{v}_{j+1}^{(n)}$. Then the update formulae (6), which follow from the minimization problem (2), (5), may be given alternatively as

[TABLE]

where $\mathcal{A}=\{(a,v^{\prime})\in\mathcal{H}_{1}\times\mathcal{H}_{1}:\widehat{C}_{j+1}a=v^{\prime}\}$ and the argmin is projected from the pair $(a,v^{\prime})$ onto the $v^{\prime}$ coordinate only. Moreover $v_{j+1}^{(n)}=\underset{\epsilon\rightarrow 0}{\mathrm{lim}}\,v_{\epsilon}$ with

[TABLE]

and $\widehat{C}_{\epsilon}=\widehat{C}_{j+1}+\epsilon I$.

Proof.

For notational convenience denote $\widehat{C}=\widehat{C}_{j+1}$. The objective function $I_{{\rm filter},j,n}(v)$ appearing in (5) is infinite if and only if $v-\widehat{v}_{j+1}^{(n)}$ is in the range of $\widehat{C}$. Thus, since the range of $\widehat{C}$ is non-empty, we may confine minimization to the set of $(a,v^{\prime})\in\mathcal{A}.$ Note that $\widehat{C}$ is in general not invertible, as it has rank $N-1$, which may be less than the dimension of $\mathcal{H}_{1}$. The set $\mathcal{A}$ thus comprises all $v^{\prime}$ in the range of $\widehat{C}$ (a convex set) and, for each such $v^{\prime}$ the set of $a$ solving $\widehat{C}a=v^{\prime}$; such an $a$ is unique upto translations in the null-space of $\widehat{C}.$ Thus $\mathcal{A}$ is a convex set. Notice that, for such pairs $(a,v^{\prime})\in\mathcal{A}$, $\langle a,v^{\prime}\rangle=|v^{\prime}|^{2}_{\widehat{C}}$ with $v^{\prime}$ lying in the range of the operator $\widehat{C}$. Although the element $a$ is uniquely defined only up to translations in the nullspace of $\widehat{C}$, such translations do not change the value of the inner product $\langle a,v^{\prime}\rangle$. The restriction of $\widehat{C}$ over the constraint set is positive definite which means that the quadratic objective function, now depending only on $v^{\prime}$, is strongly convex. Therefore the problem has a unique solution and its Lagrangian is written as:

[TABLE]

To express optimality conditions compute the derivatives and set them to zero:

[TABLE]

The last two equations imply that $\widehat{C}(2\lambda+a)=0.$ Thus we set $\lambda=-\frac{1}{2}a$ and drop the second equation, replacing the first by

[TABLE]

Solving the resulting equation for $a$ gives

[TABLE]

From this formula it follows that

[TABLE]

If we define

[TABLE]

then we see that

[TABLE]

This is precisely the form of the ensemble Kalman update, and to complete the proof of the first part of the theorem it remains to show that this defintion of $K$ agrees with the formulae given in (6); this amounts to verifying the identity

[TABLE]

To verify this we start from the matrix identity

[TABLE]

noting that it may be factored to write

[TABLE]

Inverting $(H\widehat{C}H^{T}+\Gamma)$ on the right gives the desired identity (12).

We now study the alternative representation of the minimization problem (2), (5), by (8). We first note that $H^{T}\Gamma^{-1}H+\widehat{C}_{\epsilon}^{-1}$ is strictly positive definite and hence the related quadratic function is strongly convex. As a consequence we have existence and uniqueness of the solution, and the optimality condition becomes,

[TABLE]

Then if we apply Woodbury matrix identity we obtain

[TABLE]

Note that the matrix multiplying $\widehat{v}_{j+1}^{(n)}$ is

[TABLE]

and that the matrix multiplying $y_{j+1}^{(n)}$ is

[TABLE]

so that

[TABLE]

Finally, as $A\mapsto A^{-1}$ is continuous over the set of invertible matrices, letting $\epsilon\rightarrow 0$ gives:

[TABLE]

which concludes the proof. ∎

2.2.2 Formulation In Range Of The Covariance

The minimization problem for each individual particle has a solution which, when suitably shifted, lies in the range of the empirical covariance. This allows us to seek the solution of the minimization problem as a linear combination of a given set of vectors, and to minimize over the scalars which define this linear combination. This reformulation of the optimization problem is widely employed in a variety of applications, such as weather forecasting, where the number of ensemble members $N$ is much smaller than the dimension of the data space; this is because the inversion of $S$ to form the Kalman gain $K$ takes place in the data space.

In order to implement the minimization in the $N$ dimensional subspace we note that $I_{{\rm filter},j,n}(v)$ is infinite unless

[TABLE]

for some $a\in{\mathbb{R}^{n}}.$ From the structure of $\widehat{C}_{j+1}$ given in (1c) it follows that

[TABLE]

Here each unknown parameter $b_{m}\in\mathbb{R}$ and $b:=\{b_{m}\}_{m=1}^{N},$ is the unknown vector to be determined. This form for $v$ follows from the fact that

[TABLE]

which in turn implies that

[TABLE]

Note that the unknown vector $b$ depends on $n$ as we need to solve the constrained minimization problem for each of the particles, indexed by $n=1,\dots,N$; we have suppressed the dependence of $b$ on $n$ for notational simplicity.

The expression (13) for $v$ in terms of the $e^{(m)}$ can be substituted into (5) to obtain a functional $J_{{\rm filter},j,n}(b)$ to be minimized over $b\in{\mathbb{R}^{N}},$ because $v$ is an affine function of $b.$ Equation (13) may be written in compact form as

[TABLE]

where $B$ is the linear mapping from ${\mathbb{R}^{N}}$ into $\mathcal{H}_{1}$ defined by

[TABLE]

We now identify $J_{{\rm filter},j,n}(b)$. We note that (15) is solved by taking

[TABLE]

Although $a$ is not unique, the non-uniqueness stems only from translations in the nullspace of $\widehat{C}_{j+1}$. Translations do affect the values taken by the $b_{m}$, but do not affect the vector $v$ given by (16) because they result in changes to $b$ which are in the null-space of $B$. Furthermore, for any such solution, independently of which $a$ is chosen,

[TABLE]

Using this and (16) in the definition of $I_{{\rm filter},j,n}(v)$ we obtain

[TABLE]

and hence, from (5),

[TABLE]

Once $b$ is determined it may be substituted back into (16) to obtain the solution to the minimization problem.

The preceding considerations also yield the following result, concerning the unconstrained Kalman minimization problem; its proof is a corollary of the more general Theorem 2.4 from the next subsection, which includes constraints in the minimization problem.

Corollary 2.3**.**

Suppose that the dimensions of $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are finite. Given the prediction (1a), the unconstrained Kalman update formulae may be found by minimizing $J_{{\rm filter},j,n}(b)$ from (17) with respect to $b$ and substituting into (16).

We summarize the ensemble Kalman state estimation algorithm, using minimization over the vector $b$, in the following pseudo-code:

2.3 Constrained Ensemble Kalman Filter

In this subsection we introduce linear equality and inequality constraints on the state variable into the ensemble Kalman filter. We make prediction according to (1), and then incorporate data by solving the minimization problem (5) subject to the additional constraints

[TABLE]

Here $F$ and $G$ are linear mappings which, respectively, take the state $v$ into the number of equality and inequality constraints; the notation $\preceq$ denotes inequality componentwise.

2.3.1 Formulation In The Original Variables

The preceding considerations lead to the following algorithm for ensemble Kalman filtering subject to constraints. The existence of a solution to the constrained minimization follows from Theorem 2.4 below.

2.3.2 Formulation In Range Of The Covariance

The linear constraints (18) can be rewritten in terms of the vector $b$, by means of (16), as follows:

[TABLE]

We may thus predict and then optimize the objective function $J_{{\rm filter},j,n}(b)$, given by (17), subject to the constraints (19). Implementation of this leads to following algorithm for ensemble Kalman filtering subject to constraints:

Justification for the use of this algorithm, working in the constrained space parameterized by $b$, is a consequence of the following:

Theorem 2.4**.**

Suppose that the dimensions of $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are finite. The problem of finding $v^{(n)}_{j+1}$ as the minimizer of $I_{{\rm filter},j,n}(v)$ subject to the constraints (18) is equivalent to finding $b$ to minimize $J_{{\rm filter},j,n}(b)$ subject to the constraints (19) and then using (16) to find $v_{j+1}^{(n)}$ from $b$. Furthermore, both of these constrained minimization problems have a unique solution provided that the constraint sets are non-empty.

Proof.

For notational convenience set $\widehat{v}=\widehat{v}^{(n)}_{j+1}$, $y=y^{(n)}_{j+1}$, $y^{\prime}=y-H\widehat{v}$, $\widehat{C}=\widehat{C}_{j+1}$ and $\widehat{C}_{\epsilon}=\widehat{C}_{j+1}+\epsilon I$.

Denote

[TABLE]

and

[TABLE]

The part of the statement of Theorem 2.4 concerning existence of a minimizer is a consequence of the Lemma 2.5 stated and proved below. The second part, concerning the equivalence of minimization over $b$ and over $v$ (or $v^{\prime}$) was shown in equations (13)-(17). This concludes the proof. ∎

Notice that in the following lemma, the variables $v_{\epsilon}$, $\widehat{v}$, and $v^{*}$ are analogous to $v_{j+1}^{(n)}$, $\widehat{v}_{j+1}^{(n)}$, and the minimizer over $v^{\prime}$ from (7), respectively.

Lemma 2.5**.**

Suppose that the constraint sets of (20) and (21) are non empty, then $v^{*}$ exists and is unique and for all $\epsilon>0$, $v_{\epsilon}$ exists and is unique. Furthermore $\underset{\epsilon\rightarrow 0}{lim\,v_{\epsilon}}={\widehat{v}+v^{*}}$.

Proof.

The proof is broken into two parts. In the first we prove existence and uniqueness of a solution by using the idea that the constraint relating $a$ and $v^{\prime}$ renders the problem convex. In the second part of the proof we study the $\epsilon\to 0$ limit of the regularized solution, extracting convergent subsequences through compactness, and demonstrating that they converge to the desired limit. To prove existence and uniqueness of the solution of (20), notice that it can be reformulated as

[TABLE]

and that the restriction of $\widehat{C}$ over its range is strictly positive definite. Hence $J$ is a strongly convex function being minimized over a non empty closed convex set. From standard theory $v^{*}$ exists and is unique. Then as $\widehat{C}_{\epsilon}$ is strictly positive definite, the same type of arguments provide existence and uniqueness of $v_{\epsilon}$.

Now we prove the second part of the lemma. We note that ${\widehat{v}+v^{*}}$ matches the constraints of (21). It follows that for all $\epsilon>0$, $J_{\epsilon}(v_{\epsilon})\leq J_{\epsilon}({\widehat{v}+v^{*}})$. Then let us prove that $J_{\epsilon}({\widehat{v}+v^{*}})\underset{\epsilon\rightarrow 0}{\rightarrow}J({\widehat{v}+v^{*}})$. First denote by $\lambda_{1}\leq\dots\leq\lambda_{N-1}$ the strictly positive eigenvalues of $\widehat{C}$ (recall that $\widehat{C}$ is symmetric positive semidefinite and that $\mathrm{rank}(\widehat{C})=N-1$ almost surely). Hence $\widehat{C}_{\epsilon}^{-1}=\sum_{k=1}^{N-1}\frac{1}{\lambda_{k}+\epsilon}a_{k}a_{k}^{T}+\sum_{k=N}^{\mathrm{dim}(\mathcal{H}_{1})}\frac{1}{\epsilon}a_{k}a_{k}^{T}$ where the $a_{k}$’s are the eigenvectors of $\widehat{C}$ (the first and second sums respectively gather the vectors of the range and of the nullspace of $\widehat{C}$) . As $v^{*}$ lies in the range of $\widehat{C}$, it holds that $|v^{*}+\widehat{v}-\widehat{v}|_{\widehat{C}_{\epsilon}}^{2}=\,|v^{*}|_{\widehat{C}_{\epsilon}}^{2}=\sum_{k=1}^{N-1}\frac{1}{\lambda_{k}+\epsilon}(a_{k}^{T}v^{*})^{2}$. Now as the $a_{k}$’s do not depend on $\epsilon$, by letting $\epsilon$ tending to zero, this quantity will tend to

[TABLE]

Therefore it holds that $J_{\epsilon}({\widehat{v}+v^{*}})\underset{\epsilon\rightarrow 0}{\rightarrow}J({\widehat{v}+v^{*}})$. From this we deduce that there exists $\delta>0$ such that for all $0<\epsilon<\delta$, $J_{\epsilon}(v_{\epsilon})\leq J({\widehat{v}+v^{*}})+1$.

Then set $w_{\epsilon}=v_{\epsilon}-\widehat{v}=w_{\epsilon}^{0}+w_{\epsilon}^{1}$ where $w_{\epsilon}^{0}$ lies in the nullspace of $\widehat{C}$ and $w_{\epsilon}^{1}$ in its range (recall that for a symmetric matrix nullspace and range are orthogonal) and see that $J_{\epsilon}(v_{\epsilon})=\frac{1}{2}|y^{\prime}-Hw_{\epsilon}|_{\Gamma}^{2}+\frac{1}{2}|w_{\epsilon}|_{\widehat{C}_{\epsilon}}^{2}$. It holds that $\frac{1}{2}|w_{\epsilon}|_{\widehat{C}_{\epsilon}}^{2}\leq J_{\epsilon}(v_{\epsilon})\leq J({\widehat{v}+v^{*}})+1$ for $\epsilon$ sufficiently small. Furthermore $|w_{\epsilon}|_{\widehat{C}_{\epsilon}}^{2}=|w_{\epsilon}^{0}|_{\widehat{C}_{\epsilon}}^{2}+|w_{\epsilon}^{1}|_{\widehat{C}_{\epsilon}}^{2}=\frac{1}{\epsilon}|w_{\epsilon}^{0}|^{2}+|w_{\epsilon}^{1}|_{\widehat{C}_{\epsilon}}^{2}$, and since this quantity is bounded from above we deduce that $w_{\epsilon}^{0}\underset{\epsilon\rightarrow 0}{\rightarrow}0$ and that $w_{\epsilon}^{1}$ is bounded. Let $(\epsilon_{m})_{m\in\mathbb{N}}$ be a sequence of positive real numbers such that $\epsilon_{m}\underset{m\rightarrow\infty}{\rightarrow}0$, and from the preceding extract a converging subsequence (denoted $(\epsilon_{m})_{m\in\mathbb{N}}$ for simplicity) such that $(w_{\epsilon_{m}}^{1})_{m\in\mathbb{N}}$ converges to a limit denoted $w^{*}$. As $w_{\epsilon_{m}}^{1}$ lies in $\mathcal{R}(\widehat{C})$, we can use the eigenvalue decomposition of $\widehat{C}$ to show that $|w_{\epsilon_{m}}^{1}|_{\widehat{C}_{\epsilon_{m}}}^{2}\underset{m\rightarrow\infty}{\rightarrow}|w^{*}|_{\widehat{C}}^{2}$. This limiting identity, and the fact that $w_{\epsilon}^{0}$ has limit [math], may be used to establish the first equality within the following chain of equalities and inequalities:

[TABLE]

Now note that $w^{*}$ matches all the constraints of (20). Indeed $w_{\epsilon_{m}}^{1}$ lies in the range of $\widehat{C}$ which is a closed space, also $v_{\epsilon_{m}}-\widehat{v}=w_{\epsilon_{m}}^{0}+w_{\epsilon_{m}}^{1}\underset{m\rightarrow\infty}{\rightarrow}w^{*}$. It is clear that $v_{\epsilon_{m}}-\widehat{v}$ matches the equality and inequality constraints of (20) for all $m$ and hence passing to the limit we have that $w^{*}$ satisfies the equalities and inequalities.

From the uniqueness of the minimizer of (20) we have that $w^{*}$ is equal to $v^{*}$. In particular this means that $v^{*}$ is the unique cluster point of the original sequence $(w_{\epsilon_{m}}^{1})_{m\in\mathbb{N}}$. Since the original sequence was arbitrarily chosen, we conclude that $\underset{\epsilon\rightarrow 0}{{\rm lim}\,v_{\epsilon}}={\widehat{v}+v^{*}}$. ∎

Remark 2.6.

Notice that the proof remains true if we take general convex inequalities. We simply need the constrained sets to be closed and convex; however we have restricted to linear equality and inequality constraints for simplicity and because these arise most often in practice.

3 Ensemble Kalman Inversion

3.1 Inverse Problem

In this section we show how a generic inverse problem may be formulated as a partially observed dynamical system. This enables the machinery from the preceding section 2 to be used to solve inverse problems.

We are interested in the inverse problem of finding $u\in\mathcal{H}_{1}$ from $y\in\mathcal{H}_{2}$ where

[TABLE]

Time does not appear (explicitly) in this equation (although $G$ may involve solution of a time-dependent differential equation, for example). In order to use the ideas from the previous section, we introduce a new variable $w=G(u)$ and rewrite the equation as

[TABLE]

The key point about writing the equation this way is that the data $y$ is now linearly related to the variable $v=(u,w)^{T}$ and now we may apply the ideas of the previous section to the model by introducing the following dynamical system, taking $y_{j+1}=y$ as the given data:

[TABLE]

If we introduce the new variables

[TABLE]

and write $v_{j}=(u_{j},w_{j})^{T}$, we may write the dynamical system in the form

[TABLE]

which is exactly in the same form as in the previous section. We note that

[TABLE]

3.2 Ensemble Kalman Inversion

The prediction step and the Kalman gain are defined as in (5), and the solution of the optimization problem is given by (6). We now simplify these formulae using the specific structure on $\Psi$, $v$, $H$ arising in the inverse problem and given in (22); this results in block form vectors and matrices. First we note that

[TABLE]

Here

[TABLE]

and

[TABLE]

The covariance $C^{ww}_{j+1}$ denotes the empirical covariance of the ensemble in data space, $C^{uu}_{j+1}$ denotes the empirical covariance of the ensemble in space of the unknown $u$, and $C^{uw}_{j+1}$ denotes the empirical cross-covariance from data space to the space of the unknown.

Noting that $S_{j+1}=(C_{j+1}^{ww}+\Gamma)^{-1}$ we obtain

[TABLE]

Combining equation (24) with the update equation within (6) it follows that

[TABLE]

and

[TABLE]

and hence that

[TABLE]

Thus we have derived the EKI update formula:

[TABLE]

We note also that

[TABLE]

However $w_{j+1}^{(n)}$ is not needed to update the state and so plays no role in this unconstrained EKI algorithm. (It may be used, however, to impose constraints on observation space, as discussed in the next subsection.)

In summary we have derived the following algorithm for solution of the unconstrained inverse problem:

3.3 Ensemble Kalman Inversion With Constraints

3.3.1 Formulation In The Original Variables

We now consider imposing constraints on the optimization step arising in ensemble Kalman inversion. As in the unconstrained case we do this by formulating the problem as a special case of the partially observed dynamical system, subject to constraints, from the previous section.

To this end we formulate the constraints in the space of the unknown and the data as follows:

[TABLE]

The algorithm proceeds by predicting according to equation (1), and then optimizing (5), all using the specific structure (22), and with the optimization subject to the constraints (27), written in the notation of the general Kalman updating formulae in (29), detailed below; in particular the rewrite (29) of the constraints expresses everything in terms of the variable $v$. We may summarize the constraints as follows, to allow direct application of the ideas of the previous section. To this end define

[TABLE]

Then the constraints (27) may be written as

[TABLE]

See Algorithm 6 for the resulting pseudo-code.

3.3.2 Formulation In Range Of The Covariance

We describe an alternative way to approach the derivation of the EKI update formulae. We apply Theorem 2.4 with the specific structure (22), (29) arising from the dynamical system used in EKI. To this end we define

[TABLE]

where $b$ is the vector of $N$ scalar weights $b_{m}$ and

[TABLE]

Once this quadratic form has been minimized with respect to $b$ then the update formula (13) gives

[TABLE]

Note that the vector $\{b_{m}\}$ depends on the particle label $n$; as in the previous section, we have suppressed this dependence for notational convenience. We may now impose linear equality and inequality constraints on both $u$ and $w=G(u)$ (i.e. in parameter and data spaces) and minimize (30) subject to these constraints. To be more specific if we impose the constraints (29) expressed in the variable $b$:

[TABLE]

Here $F,G,f$ and $g$ are given by (28), $B$ is defined by (31) and

[TABLE]

See Algorithm 7 for the resulting pseudo-code.

Remark 3.1.

As in the previous section, the result holds true for general convex inequality constraints; the linear case is considered for simplicity of exposition, and because it is most frequently arising in practice.

Remark 3.2.

The EKI algorithm, with or without constraints, has the following invariant subspace property: define $\mathcal{A}=span(u^{(n)}_{0})_{n\in\{1,\cdots,N\}}$, then for all $j$ in $\{0,\dots,J\}$ and for all $n$ in $\{1,\cdots,N\}$, then the $u^{(n)}_{j}$ defined by the three algorithms in this section all lie in $\mathcal{A}$. This is a direct consequence of writing the update formulae in terms of $b$ and noting (32).

We can now state a result analogous to Theorem 2.4, and with proof that is a straightforward corollary of that result, using the specific structure (22):

Theorem 3.3**.**

Suppose that the dimensions of $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ are finite. Suppose also that the specific structure (22) is applied. The problem of finding $u_{j+1}^{(n)}$ from the minimizer of $I_{{\rm filter},j,n}(v)$, defined in (5) and subject to the constraint (27), is equivalent to finding $b$ that minimizes (30), subject to (33), and then using (32) to find $u_{j+1}^{(n)}$ from $b$. Furthermore, both of these constrained minimization problems have a unique solution provided that the constraint sets are non-empty.

4 Numerical Results

This section contains numerical results which demonstrate the benefits of imposing constraints on ensemble Kalman methods. Subsection 4.1 concerns an application of state estimation (using EnKF) in biomedicine, using real patient data, whilst subsection 4.2 concerns an application of inversion (using EKI) in seismology and employs simulated data. When comparing results from the two experiments, recall that iterations of EKI correspond to an algorithmic dynamics intended to converge to a single distribution (over ensemble members) on the parameters for which we invert, whereas iterations of EnKF correspond to the incorporation of new data at every physical measurement time, and thus the distribution (over ensemble members) is not necessarily expected to converge as the iteration progresses. In both applications, minimizations were performed in Matlab using the default interior-point method (fmincon) through the general quadratic programming function, quadprog.

4.1 State Estimation

Here we present an application of the constrained EnKF to the tracking and forecasting of human blood glucose levels. We use self-monitoring data collected by an individual with Type 2 Diabetes. We use the "P1" data set described by Albers et al. in [34]; this dataset includes measurements of blood glucose and consumed nutrition, and is publicly available on physionet.org. For more information on the data, and on an unconstrained data assimilation approach using the unscented Kalman filter, see [34]. We model the glucose-insulin system with the ultradian model proposed by [35]. The primary state variables are the glucose concentration, $G$, the plasma insulin concentration, $I_{p}$, and the interstitial insulin concentration, $I_{i}$; these three state variables are augmented with a three stage delay $(h_{1},h_{2},h_{3})$ which encodes a non-linear delayed hepatic glucose response to plasma insulin levels. The resulting ordinary differential equations have the form:

[TABLE]

Here $m_{G}(t)$ represents a known rate of ingested carbohydrates, $f_{1}(G)$ represents the rate of glucose-dependent insulin production, $f_{2}(G)$ represents insulin-independent glucose utilization, $f_{3}(I_{i})G$ represents insulin-dependent glucose utilization and $f_{4}(h_{3})$ represents delayed insulin-dependent hepatic glucose production; the functional forms of these parameterized processes can be found in the appendix, along with a description of model parameters.

In the EnKF setting, we write $u=[I_{p},I_{i},G,h_{1},h_{2},h_{3}]$, and use (34) to define $F$ such that

[TABLE]

where $\theta$ contains model parameters. We then extend the state vector in order to perform joint parameter estimation: $v=[u,R_{g}]^{T}$.

For the purposes of this paper, the function $m_{G}(t)$ may be viewed as known; it is determined from data describing meals consumed by the patient. Since insulin ($I_{p}$ and $I_{i}$) and delay variables ($h_{1}$, $h_{2}$, and $h_{3}$) are not measured, whilst glucose is measured, we define the measurement operator to be $H=[0,0,1,0,0,0,0]$. The discrete time forward model is obtained by integrating the deterministic model in (34) between consecutive measurement time-points and applying an identity map to $R_{g}$. Because these time-points may not be equally spaced, and because the time-dependent forcings (meals) will differ in different time-intervals, this leads to a map of the form

[TABLE]

This is a slight departure from the methodology outlined in section 2, where $\Psi$ does not depend on $j$ (autonomous dynamics) but is a straightforward extension which the reader can easily provide.

We present EnKF results from a single patient’s data when run with and without constraints (Algorithms 1 and 3 respectively). We performed joint state-parameter estimation, augmenting the state with parameter $R_{g}$ (see Appendix for details of where this parameter appears) and adding identity-map dynamics for parameter $R_{g}.$ The following constraints were imposed:

[TABLE]

Figure 1 compares the overall distribution of updated state means over time when running EnKF with and without these state constraints. While individual particles in this experiment often violated the constraints, the overall updated means did not. Nevertheless, enforcement of lower-bound constraints shifts up the state distribution slightly. Note that upper bound constraints were never violated in this experiment.

Figure 2 shows a two-dimensional state projection of updated particles at a given time step before and after applying the constrained optimization. Note that particles may additionally violate constraints in unplotted dimensions—this explains why one particle whose unconstrained update appears to live within the constraints is in fact differently updated under the constrained optimization. Time step 126 was selected for illustrative purposes, and was the measurement event in which particles most often violated the constraints.

Figure 3 depicts the overall frequency of constraint violations. We observe that the the measured state (blood glucose) never violated a constraint, nor did the inferred parameter $R_{g}$. However, other model states did often violate constraints, and up to $30\%$ $(4/13)$ of particles simultaneously violated the constraints at a single time-step.

By adding constraints, we ensure that all the simulations which constitute the ensemble method are biologically plausible.

4.2 Inverse Problem

Here we present application of the constrained EKI in seismology. We study near-surface site characterization in which we invert for the shear wave velocity profile of the geomaterials in the earth shallow crust, using downhole array data. For forward modeling, we consider a semi-discrete form of the following wave equation in a horizontally stratified heterogeneous soil layer:

[TABLE]

Here $d(z,t)$ is the displacement field of the wave response as a function of spatial variable $z\in(0,H)$ and time variable $t\in(0,T]$. The function $c_{s}(z)$ is the shear wave velocity function. We impose the following boundary and initial conditions:

[TABLE]

where $d_{0}(t)$ is the prescribed displacement at depth $z=H$. Generally, the shear wave velocity changes as a piecewise constant function with depth. If the layering information, i.e., the total number of layers and their thickness, is not available or is poorly characterized, it is desired to use a generic function for site characterization, such as this:

[TABLE]

See, for example, [36]. In the constrained EKI setting, $u=(c_{s0},k,z_{0},n,z_{1},\alpha)$ and

[TABLE]

For the numerical example studied here, $G^{u}$ and $g^{u}$ are determined by enforcing the constraints $0\leq c_{s0}\leq 1000$, $0\leq k\leq 100$, $0\leq z_{0}\leq z_{1}$, $0\leq n\leq 1$, $z_{0}\leq z_{1}\leq H$, and $1\leq\alpha\leq 10$. We generate the initial ensemble by drawing samples from uniform distributions and discard members that violate the enforced constraints. In order to avoid very large velocities at $z=z_{1}$, we also discard members with $c_{s}(z_{1})>5000$m/s. If we perform parameter learning using the unconstrained EKI, the experiment fails at $j=1$ because of incapability of the dynamic model to propagate unphysical values of the shear wave velocity $c_{s}$.

All results shown use Algorithm 7. Figure 4 shows the ensemble distribution of $u$ at $j=2$ before and after enforcing constraints whilst Figure 5 shows the evolution of the updated ensemble. Note that parameter $k$ saturates with an ensemble close to the upper bound of $100$ imposed through constraints on this parameter; however experiments in which we imposed different upper bounds on this parameter lead to different estimates for $k$, with little change to the estimated velocity profile and we conclude that this parameter suffers from identifiability issues. (Note that Figure 4 displays the updated ensemble distribution at a single step in the sequence of ensemble updates, comparing the effect of imposing constraints with neglecting them; in contrast Figure 1 shows the distribution over all measurement time-points of the ensemble means. The figures thus illustrate different phenomena).

Moreover, Figure 6a shows the map of violation for different constraints enforced on parameters whilst Figure 6b shows the estimated generic $c_{s}$ profile after 40 iterations compared to the true profile and the initial estimate. Figure 6a shows the key role employed by the enforcing of constraints. In this case the addition of constraints ensures that all the simulations which constitute the ensemble method are physically meaningful, and also that the forward model remains well-posed.

5 Conclusions

Constraints arise naturally in many state and parameter estimation problems. We have shown how convex constraints may be incorporated into ensemble Kalman based state or parameter estimation algorithms with relatively few changes to existing code: the standard algorithm is applied and for any ensemble member which violates a constraint, a quadratic optimization problem subject to convex constraints is solved instead. We have written the resulting algorithms in easily digested pseudo-code, we have developed an underpinning theory and we have given illustrative numerical examples.

Two primary directions suggest themselves in this area. The first is the use of these methods in applications. As indicated in the introduction, our general formulation is inspired by the two papers [8, 6] from the geosciences and we have demonstrated applicability to problems from biomedicine and seismology; but many other potential application domains are ripe for application of ensemble Kalman methodology, because of its black-box and derivative-free formulation, and the ability to impose constraints in a straightforward fashion will help to extend this methodology. The second is the theoretical analysis of these methods: can the inclusion of constraints be used to deduce improved accuracy of state or parameter estimates; or can the inclusion of constraints be used to demonstrate improved performance as measured, for example, by proportion of model runs which are physically (or biologically etc.) plausible? Furthermore, although the imposition of constraints is reasonable, it is not clear that it may not lead to pathologies in algorithmic performance and ruling out, or understanding, the occurrence of such pathological behaviour may be important.

This work was funded by NIH-NLM grant RO1 LM012734. AMS was also funded by AFOSR Grant FA9550-17-1-0185 and by ONR grant N00014-17-1-2079.

Appendix

We give the details of the ultradian model of glucose-insulin dynamics used as the forward model in subsection 4.1. An example of the induced dynamics is given in Figure 7.

[TABLE]

where, for $N$ meals at times $\{t_{j}\}_{j=1}^{N}$ with carbohydrate composition $\{m_{j}\}_{j=1}^{N}$

[TABLE]

and

[TABLE]