The motivation for writing this article was Research Debt by Chris Olah and Shan Carter. With it I am trying to reduce my research debt by writing an "extended introduction" to two manuscripts I have worked on:
- [1] "A strongly convergent numerical scheme from Ensemble Kalman inversion" by Dirk Blömker, Claudia Schillings, and myself.
- [2] "Well Posedness and Convergence Analysis of the Ensemble Kalman Inversion" by Dirk Blömker, Claudia Schillings, Simon Weißmann, and myself.
Warning 1: In order to convey the general ideas (i.e. the "gist"), we will at some points sacrifice mathematical exactness for clarity. For the rigorous details you can always consult the papers [1,2].
Warning 2: This is still work in progress.
Bayesian inverse problems
We want to solve the following inverse problem: Given an observation \(y\) given by
vation
\(\in\mathbb R^K\) unknown parameter \(\in \mathcal X\)
operator
\(\mathcal X\to\mathbb R^K\) observational noise \(\mathcal N(0,\Gamma)\)
we want to infer the unknown parameter \(u.\)
\[ \text{Inverse Problem:}\qquad y \text{ (given)} \stackrel{?}{\mapsto} u, \]The problem here is that \(y\) is polluted by (unknown) noise and that \(\mathcal G\) might be strongly non-injective (i.e. not one-to-one). Thus we cannot do this: \(u = \mathcal G^{-1}(y - \eta).\)
Formally, we can write \( y - \mathcal G(u) = \eta \sim \mathcal N(0,\Gamma) \) and hence the likelihood functional is given by
\[\begin{aligned}\text{Likelihood:}\qquad\mathbb P(y|u) &= C \cdot \exp\left( -\tfrac{1}{2}\left\|\Gamma^{-\frac{1}{2}}(y - \mathcal G(u))\right\|^2\right) \\ & =: C\cdot \exp\left( -\Phi(u;y)\right) . \end{aligned}\]The function \(\Phi(u;y)\) is also called the negative loglikelihood.
If we employ a Bayesian approach, we define a prior distribution \(\mu_0\) (we could also write this as \(\mathbb P(u)\)) on parameter space and by Bayes' theorem\(\mathbb P(u|y) \propto \mathbb P(y|u)\cdot \mathbb P(u)\) we get the posterior distribution \(\mu\) (alternatively, \(\mathbb P(u|y)\)) by
\[\text{Posterior: }\qquad\mu(\mathrm d u) \propto \exp\left( -\Phi(u;y)\right) \cdot \mu_0(\mathrm d u) \]In other words, while we might think of the inverse problem as
which is not solvable due to the noise term and possible non-invertibility of \(\mathcal G,\) the Bayesian approach suggests to "settle" for
\[\text{Bayesian solution:}\qquad (y, \mu_0) \mapsto \mu. \]Are we done?
We know that the problem is solvable. Practically, for most real-world applications, the analytical form of the posterior is still too complicated to make use of it. One way of actually working with the posterior is to employ Monte Carlo methods. Here, we will think about the Ensemble Kalman filter instead.
The Kalman filter
Easiest form of the Kalman method:
\[\begin{aligned} \text{unknown parameter:}\qquad u &\sim \mathcal N(m, C)\\ \text{observational noise:}\qquad \eta&\sim \mathcal N(0, \Gamma)\\ \text{observation:} \qquad y &= A \cdot u + \eta \end{aligned}\]Then with "Kalman gain" \(K = CA^T\left(\Gamma + ACA^T\right)^{-1}\), we obtain
Example: Infer form of linear functions from noisy measurements
We consider a linear function \(f: \xi\mapsto m\cdot \xi + t.\) where the parameters \(m,t\) are unknown. We choose Gaussian prior with variance \(\sigma^2\) on both parameters. Now we get some information on the parameters by observing \(f\) at \(\xi = -1,0,1\), perturbed by some observational noise. In the language from above, we can state this as follows:
Infer \(m, t\) in \[f: \xi\mapsto m\cdot \xi + t\] from data points \(f(x_i) + \eta_i\) at all three positions \(x_i = -1,0,1.\)
\[\begin{aligned}\text{parameter: }u &= \begin{pmatrix}m \\ t\end{pmatrix} \sim \mathcal N\left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma^2 & 0 \\ 0 & \sigma^2 \end{pmatrix} \right) \\ \text{observational noise: }\eta &\sim \mathcal N \left(\begin{pmatrix}0\\0\\0\end{pmatrix},\begin{pmatrix} \gamma^2 & 0 & 0 \\ 0 &\gamma^2 & 0\\ 0 & 0 & \gamma^2\end{pmatrix}\right)\\ \text{evaluation matrix: }A & = \begin{pmatrix} -1 & 1 \\ 0 & 1 \\ 1 & 1 \end{pmatrix}\\ \text{noisy evaluation: }y &= A\cdot u + \eta\end{aligned} \]Instructions
You can play with this example by using the controls below:
- Choose "true" parameter values The notion of true parameters is not strictly Bayesian. \(m,t\).
- Choose a observational noise level \(\gamma\).
- Generate a noisy observation by clicking "Sample points".
- Let the Kalman filter compute the posterior distrbution of the parameters \(m,t\) given the data by clicking "Make inference".
- Sample parameter values \(\tilde m, \tilde t\) from the posterior distribution by clicking "Sample Posterior" and plot the corresponding function \(\xi \mapsto \tilde m\cdot \xi + \tilde t.\) Check whether this is a plausible fit for the data (keep in mind your set noise level).
Algorithm only uses this
This is a great way to do Bayesian inference: It is relatively straightforward, can be applied to systems that change in time, for example for position estimation in autonomous vehicles or for parameter estimation of dynamical systems. The downside of this approach is that it is only exact for linear forward operators as above and it does not scale well to high dimensions. Usually, the Kalman filter is applied to a dynamical system which is time-adaptively "filtered". This means that the calculation of the Kalman gain which involves inversion of a high-dimensional matrix, has to be done in every time step.
Caveat: Only exact for linear forward operator and costly calculation of covariance update (unfeasible in realistic settings).
The Ensemble Kalman filter (EnKF)
The EnKF does not work with a probability distribution , but it manages a collection of "particles" comprising an "ensemble" which is a surrogate for the distribution.
Kalman filter: . EnKF: .
The (normal) Kalman filter does this:
Calculate posterior mean and covariance from prior mean and covariance (using observation).
Kalman filter: prior mean, covariance \(\mapsto\) posterior mean, covariance
In contrast, the Ensemble Kalman filter does this:
Evolve ensemble particles. Mean and covariance are implicitly given by (empirical) mean and covariance of this ensemble.
EnKF: prior ensemble particles \(\mapsto\) posterior ensemble particles
Let's apply the EnKF method to our simple problem of inferring an unknown parameter \(u \in \mathbb R^n\) where we only know a noisy observation \(y = Au+z\). Here, \(z \sim N(0,\Gamma)\) is an (unknown) noise term with (known) covariance.
- ChooseYou could do this by sampling from your prior distribution. an initial
ensemble with \(J\) members of "prior" parameters, or particles, \(\{{u^{(i)}}\}_{i=1}^J\) with \({u^{(i)}} \in \mathbb R^n\). - Artificially perturbThere is good reason for that: First, this "relativizes" the strength of the data as perceived by the EnKF. It makes sense not to enforce 100 identical instances of the same data if this data is not exact. In contrast, 100 "wiggly" data points with "wiggle width" similar to the original noise's standard deviation gives a more realistic image. Second, this will be necessary in order for the posterior empirical covariance to be "correct" as we will see later. the data \(y\) for each particle: \[{y^{(i)}} = y + \epsilon_i\] where \(\epsilon_i \sim N(0, \Gamma)\) are sampled.according to the same measure as the original observation error \(z\).
- The ensemble of particles has an initial empirical mean \[\bar u = \frac{1}{J}\sum_{i=1}^J {u^{(i)}}\] and initial empirical covariance In case you were wondering: We use the biased version with \(\frac{1}{J}\) instead of \(\frac{1}{J-1}\) for simplicity. \[C^\text{emp} = \frac{1}{J}\sum_{i=1}^J ({u^{(i)}}-\bar u) \cdot ({u^{(i)}}-\bar u)^T\]
- Now we transform each particle (we write \({u^{(i)}}^*\) for the new particle position): \[{u^{(i)}}^* = {u^{(i)}} + C^\text{emp}A^T\left(AC^\text{emp}A^T + \Gamma\right)^{-1}({y^{(i)}} - A{u^{(i)}})\]
- This new ensemble implicitly constitutes a posterior measure with empirical mean \[ \bar u^* = \frac{1}{J}\sum_{i=1}^J {u^{(i)}}^* \] and covariance \[C^{*\text{,emp}} = \frac{1}{J}\sum_{i=1}^J ({u^{(i)}}^*-\bar u^*) \cdot ({u^{(i)}}^*-\bar u^*)^T.\]
- Initial ensemble \(\{{u^{(i)}}\}_{i=1}^J\)
- Perturb observation: \[{y^{(i)}} = y + \epsilon_i\] where \(\epsilon_i \sim N(0, \Gamma)\)
- Empirical mean \[\bar u = \frac{1}{J}\sum_{i=1}^J {u^{(i)}}\] Empirical covariance \[C^\text{emp} = \frac{1}{J}\sum_{i=1}^J ({u^{(i)}}-\bar u) \cdot ({u^{(i)}}-\bar u)^T\]
- Map to posterior ensemble \[{u^{(i)}}^* = {u^{(i)}} + C^\text{emp}A^T\left(AC^\text{emp}A^T + \Gamma\right)^{-1}({y^{(i)}} - A{u^{(i)}})\]
Why is that the specific form of the EnKF? You may compare it to the form of the Kalman filter and see that this is essentially the "empirical" version of it. Alternatively, you can read the following sanity check box.
Sanity check: EnKF
Recall the Kalman update formula: \[u|y \sim \mathcal N(m^*, C^*)\] where \[\begin{aligned}m^* &= m + K(y-Am)\\C^* &= C-KAC\end{aligned}\] with \(K = CA^T\left(\Gamma + ACA^T\right)^{-1}\) and \(u\sim \mathcal N(m, C)\) In order for the EnKF to be an "empiricalization" of the Kalman filter, we would expect that we get something similar with the respective empirical versions.
Empirical mean:
\[\begin{aligned} \bar u^* &= \frac{1}{J} \sum_{i=1}^J {u^{(i)}}^* \\ &= \frac{1}{J} \sum_{i=1}^J \left[{u^{(i)}} + C^\text{emp}A^T\left(AC^\text{emp}A^T + \Gamma\right)^{-1}({y^{(i)}} - A{u^{(i)}})\right] \\ &= \bar u + C^\text{emp}A^T\left(AC^\text{emp}A^T + \Gamma\right)^{-1} \cdot (\bar y - A\bar u)\\ &= \bar u + K^\text{emp} (\bar y-A\bar u). \end{aligned}\]where \(K^\text{emp} = C^\text{emp} A^T \left(\Gamma + A C^\text{emp} A^T\right)^{-1}\) and \(\bar y = \frac{1}{J}\sum_{i=1}^J {y^{(i)}}.\)
This is essentially the "updating the mean" formula from above with empirical versions of the mean, covariance and Kalman gain.
Empirical covariance:
The following (lengthy) calculation shows that the empirical covariance gets mapped to
\[\begin{aligned} C^{*\text{,emp}} &=C^\text{emp} - K^{\text{emp}}AC^{\text{emp} } + R_J , \end{aligned}\]where \(R_J\) is an error term which vanishes in probability for \(J\to 0\).
In the calculation we will also see why we needed to randomize the data by defining \({y^{(i)}} = y + \epsilon_i\) instead of using the same data \(y\) for all ensemble particles \({u^{(i)}}\).
Calculation
Now all terms labeled [mixed] will be
Now we note that \(C^\text{emp}A^T K^{\text{emp,} T} = K^{\text{emp}}AC^{\text{emp} }\) and we additionally artificially include a term \(K^{\text{emp}}\Gamma K^{\text{emp,} T}\) and subtract it again:
\[\begin{aligned} C^{*\text{,emp}} &=C^\text{emp} + R_J \\ & \qquad\qquad- 2K^{\text{emp}}AC^{\text{emp} } \\ & \qquad\qquad+ K^{\text{emp}}\cdot \left[\Gamma + AC^{\text{emp} }A^T\right] \cdot K^{\text{emp,} T} \\ & \qquad\qquad- K^{\text{emp}}\left[Y - \Gamma\right] K^{\text{emp,} T} \end{aligned}\]By looking at the definition \(K^\text{emp} = C^\text{emp} A^T \left(\Gamma + A C^\text{emp} A^T\right)^{-1}\) we realize that the third line is actually
\(K^{\text{emp}}\cdot AC^{\text{emp} }\) and can be combined with the first line. The last line is
which is exactly the "empiricalized" version of the Kalman covariance update formula with error terms that stem from the randomization of the data \(\{{y^{(i)}}\}\) but vanish in probability (due to the central limit theorem and the law of large numbers) when \(J\) gets larger.
Note that we needed the randomization of the data \(\{{y^{(i)}}\}\) explicitly in order to create the term
\[ K^{\text{emp}}\cdot \left[\Gamma + AC^{\text{emp} }A^T+ Y - \Gamma\right ] \cdot K^{\text{emp,} T}\]If we hadn't randomized the \({y^{(i)}}\) by adding an error term, \(Y\) would have been zero and we would have been stuck with an error term of the form \(K^{\text{emp}}\left[0-\Gamma\right] K^{\text{emp,} T}\) which does not vanish for \(J\to \infty.\)
You can drag and drop the particles in the following figure in order to get a feeling how the ensemble constitutes a Gaussian measure with a mean and covariance matrix. The mean is drawn as a hollow circle and the covariance structure is shown by a cross centered at the mean along the matrix's two eigenvectors with a length of the square root of the matrix's eigenvalues (i.e. the standard deviation).
Interactive Box: Ensemble particles
EnKF "Interpolation"
We saw above that the Bayesian methodology for the inverse problem \(y = \mathcal G(u) + \eta\) works by mapping
Bayes: prior \(\mapsto\) posterior
\[ \mu_0 \mapsto \mu = \frac{1}{Z}\exp\left(-\Phi(u; y)\right) \mu_0\]where \(\mu_0\) is the prior on \(u\), and \(\mu\) is the posterior on \(u\) obtained by incorporating the data \(y\) with Bayes' theorem . A good idea is now to make this transition more smoothly by introducing \(N\) intermediate steps. Set \(h:= \frac{1}{N}\), then
Interpolate:
\[\begin{aligned} \mu_0 &\mapsto \mu_1 &\mapsto \cdots &\mapsto \mu_{N-1} &\mapsto \mu_N = \mu \end{aligned}\]where \(\mu_k = \frac{1}{Z_k}\exp\left(-k\cdot h\cdot \Phi(u; y)\right) \cdot \mu_0 = \frac{1}{Z_k}\exp\left(-h\cdot\Phi(u; y)\right) \cdot \mu_{k-1}\) and \(Z_k, Z_k'\) are normalization constants.
Here, again, we work in the linear setting \(\mathcal G(u) = Au.\)
Here, \(\mathcal G(u) = Au.\) again
This yields a sequence of iterations of particles in the following form:
For a detailed derivation, see Iglesias, Law, Stuart,
Then
\[{u_{n+1}^{(i)}} = {u_n^{(i)}} + h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}({y_n^{(i)}} - A{u_n^{(i)}}) \]Here, again, we perturb the observational replications \(y_n^{(i)} = y + \epsilon_n^{(i)}\) in each time step, where \(\{\epsilon_n^{(i)}\}_{n,i}\) are i.i.d. distributed according to \(\mathcal N(0, h^{-1}\Gamma)\). With that we can write \(y_n^{(i)} = y + h^{-1/2} \Gamma^{1/2} \zeta_n^{(i)}\) with \(\zeta_n^{(i)}\sim \mathcal N(0,1)\) and further:
with \(y_n^{(i)} = y + h^{-1/2} \Gamma^{1/2} \mathcal N(0,1),\) i.e.
\[\begin{aligned}{u_{n+1}^{(i)}} = {u_n^{(i)}} &+ h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}({y} - A{u_n^{(i)}})\\ &+ \sqrt h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}\Gamma^{1/2} \zeta_{n+1}^{(i)}\end{aligned}\]A remark on the specific form of the observational perturbation: This can be derived rigorously, but intuitively the reason for upscaling the noise for \(h\to 0\) is the following:
If we use the same observation more and more often, then this will increase its influence Just like telling your friends over and over again that you are not interested in yesterday's soccer match will increase the probability that you will talk abouth something else next time on the EnKF iteration. If you turn up the noise perturbation then the actual information content in the data \(y_n^{(i)}\) is diminished. The factor of \(h^{-1}\) in the variance is exactly the correct scaling.
From discrete EnKF iteration to continuous EnKF evolution (gist of [1])
We derived the discrete EnKF iteration as
\[\begin{aligned}{u_{n+1}^{(i)}} = {u_n^{(i)}} &+ h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}({y} - A{u_n^{(i)}})\\ &+ \sqrt h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}\Gamma^{1/2} \zeta_{n+1}^{(i)}\end{aligned}\]In this section we will talk about how this is a discretization scheme for the stochastic differential equation
is (formally) a discretization of
\[\mathrm d u^{(i)} = C^\text{emp}A^T\Gamma^{-1}(y-Au^{(i)}) \mathrm d t + C^\text{emp}A^T\Gamma^{-1/2}\mathrm d W^{(i)}_t\]and how it might be an
But first, we will recap a few important ideas from numerical discretizations of ODEs and SDEs.
ODE discretizations
Let's take an ordinary differential equation which is not solvable analytically, like
\[ y'(t) = \exp(-t^2).\]This has no closed form solution. Given that the \(\operatorname{erf}\) function is usually not deemed elementary but is defined by being the solution of above ODE.
Actually, differential equations are usually not analytically solvable, so a numerical approach is in order. Let's consider \[y'(t) = f(t, y(t)) ~\text{ and }~ y(0) = y_0.\]
Then the easiest way of approximating the ODE on some interval \([0,T]\) is by defining \(N+1\) discrete steps \(t_k = k\cdot \frac{T}{N}.\) The the deriative of \(y\) is approximated by the secant slope: \(f(t_k, y(t_k)) = y'(t_k) \approx \frac{y(t_{k+1})-y(t_k)}{h}\), or, iteratively \[y_{k+1} = y_k + h \cdot f(t_k, y(t_k)).\] Here, \(y_k\) is an approximation for \(y(t_k)\).
This is called the Euler discretization scheme for ODEs and it usually does not work very well.
Consider the ODE \(x' = -x^3\), which is actually analytically solvable (so we can compare solution and discretization). Try setting the parameters such that the discretization "breaks".
\(x' = -x^3\)
As you probably will have realized, for a given stepsize (i.e. a given number of steps) there is a "dangerous" level for the initial value \(x_0\). If the initial value is below this level, the discretization will be an alright approximation. If the initial value is higher than this level, runaway oscillations begin to develop. This is due too the strong amplification of the derivative \(x' = -x^3\).
A simple numerical example: Let's consider \(x_0 = 10\) and \(h = 1\). Then \[\begin{aligned}x_1 &= x_0 - h\cdot x_0^3 = 10 - 1000 = -990 \\ x_2 &= x_1 - h\cdot x_1^3 = -990 + 970299000 = 970298010\end{aligned}\] The next value will be even more strongly negative etc.
A remedy for this behaviour is choosing better discretizations, for example implicit or higher-order methods. This is not the focus here but we will need this image of "dangerous sets" which, when hit, kick off runaway iterations.
Just like for ODEs, there is a very simple discretization scheme for stochastic differential equations that you can try out first but which will fail for most interesting SDEs. This is the Euler-Maruyama scheme.
The Euler-Maruyama SDE discretization scheme
We consider an SDE
\[\mathrm d y_t = f(t, y_t)\mathrm d t + g(t, y_t) \mathrm d W_t ~\text{ and } ~ y_0 \in \R.\]Prepare a sequence of random numbers \(\epsilon_n\sim N(0,1)\), then the Euler-Maruyama scheme is the following iteration.
Simplest SDE discretization: Euler-Maruyama. Choose \(\epsilon_n\sim N(0,1)\):
\[y_{n+1} = y_n + h\cdot f(t_n, y_n) + \sqrt h \cdot f(t_n, y_n)\cdot \epsilon_{n+1}.\]Here, \(y_n\) is an approximation for \(y(t_n)\), where \(t_n = T\frac{n}{N}\) for some \(N\) and \(y(t)\) is a path (a realization) of the SDE.
The EnKF iteration as a discretization scheme
Now we consider the EnKF iteration scheme
\[\begin{aligned}{u_{n+1}^{(i)}} = {u_n^{(i)}} &+ h\cdot C^\text{emp}A^T\left({\color{blue}h\cdot AC^\text{emp}A^T} + \Gamma\right)^{-1}({y} - A{u_n^{(i)}})\\ &+ \sqrt h\cdot C^\text{emp}A^T\left({\color{blue}h\cdot AC^\text{emp}A^T} + \Gamma\right)^{-1}\Gamma^{1/2} \zeta_{n+1}^{(i)}\end{aligned}\]Now we drop the (asymptotically vanishing) terms \(h\cdot AC^\text{emp}A^T\) in the iteration and get
Drop terms \({\color{blue}h\cdot AC^\text{emp}A^T}\)
\[\begin{aligned}{u_{n+1}^{(i)}} = {u_n^{(i)}} &+ h\cdot C^\text{emp}A^T\Gamma^{-1}({y} - A{u_n^{(i)}})\\ &+ \sqrt h\cdot C^\text{emp}A^T \Gamma^{-1}\Gamma^{1/2} \zeta_{n+1}^{(i)}\end{aligned}\]This is the Euler-Maruyama discretization scheme for the
Now we can ask,
Is the EnKF iteration similar enough to the Euler-Maruyama iteration in order to be a numerical approximation scheme for the continuous EnKF evolution?
Is the EnKF iteration a "valid" discretization of the SDE (like Euler-Maruyama?)
But actually, it looks like
The EnKF iteration is a much better approximation scheme for the continuous EnKF evolution than the Euler-Maruyama scheme.
Here, "better" means that the EnKF iteration converges (supposedly, see below) in \(p\)th moments to the continuous EnKF evolution, whereas Euler-Maruyama does not.
Actually, it is better than Euler-Maruyama!
A numerical schemeWith discretization parameter \(h\). This is originally a discrete quantity but we interpolate it in a suitable way for continuous time. \(\hat u_h\) is said to converge strongly in \(p\)th moments to the solution of a continuous SDE \(u\) if \[\lim_{h\to 0}\mathbb E\sup_{t}|u(t)-\hat u(t)|^p = 0\]
Now beware that whether the EnKF iteration indeed does converge strongly in \(p\)th moments to the solution of the EnKF SDE is still an open problem but for a simplifying particular case of parameters we were able to prove this in "A strongly convergent numerical scheme from Ensemble Kalman inversion".
Claim (partially proven in [1] for a strongly simplified model):
EnKF iteration yields strongly converging discretization.
Why is Euler-Maruyama "bad"?
Of course, discretization schemes for SDEs have the same problems as the ones for ODEs but due to stochasticity we need to be even more careful (even more things can go wrong).
It is of major importance that not only almost surely (path-)element wise we have convergence of the numerical scheme to the solution of the continuous SDE solution, but that the
collection of discretization paths as a random variable "converges" to the collection of continuous paths. One notion of this kind of convergence is
For the same reasons as Euler for ODEs, but worse due to stochasticity
Now why is that a "good" notion of convergence? One reason is that convergence in \(p\)th moments is a classical convergence notion which implies convergence in probability and convergence in distribution. Second, This means that \(p\)th moments of the solution can be approximated by \(p\)th moments of the discretization. This is very important for many applications, for example, "What is the variance of this particular financial option?". If the SDE governing the evolution of the option price cannot be solved analytically and needs to be approximated, then we can only answer this question by computing the empirical variance of the discretization and using it as a surrogate.
What could go wrong? Let's consider the SDE \(\mathrm d u = -u^3\mathrm d t + \mathrm d W.\) It is as well behaved as one might wish from an SDE, and paths will "typically" never leave the set \([-5,5].\)
But if we use the Euler–Maruyama discretization scheme to approximate it, we sometimes get weird results. Convince yourself of the following observations:
- A discretization of only \(N=10\) steps in the interval \([0,10]\) is too coarse and the discretization will break for the same reason as above for the deterministic differential equation \(u'(t) = -u^3(t).\)
- \(N=20\) steps are
sometimes enough but sometimes they are not. - \(N=25\) steps are
usually enough. - \(N > 25\) steps are almost always enough.
\(\mathrm d u = -u^3\mathrm d t + \mathrm d W\)
So, we might say, "Aha, we could eliminate oscillations from the deterministic Euler scheme by taking at least [
There is some pragmatism to this line of thought but there is a huge caveat: While the probability of "oscillating path events" for a discretization sample converges (exponentially) to 0, the "badness" of divergence in those oscillating paths increases bi-exponentially. This means that we get "fewer but way, way worse outliers". As all paths "get a vote" in the calculation of the empirical \(p\)th moments, this leads to divergence of those quantities. For a proof of this, see "Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients" by Hutzenthaler, Jentzen, and Kloeden. See also "Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients" by Hutzenthaler and Jentzen.
Euler-Maruyama for non-Lipschitz drift and diffusion terms: Explosion probability: exponentially small in stepsize but bi-exponentially "bad".
This is the reason that the Euler-Maruyama discretization does not converge strongly in \(p\)th moments.
Taming schemes and the EnKF iteration
So what can we do about diverging outliers in discretization schemes? Obviously, we can just remove all "terrible" paths as outliers. Another option is "taming", see again "Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients" by Hutzenthaler and Jentzen. An example for a tamed scheme is the following:
\[y_{n+1} = y_n + \frac{h\cdot f(t_n, y_n) + \sqrt h \cdot f(t_n, y_n)\cdot \epsilon_{n+1}}{1 + h\cdot |h\cdot f(t_n, y_n) + \sqrt h \cdot f(t_n, y_n)\cdot \epsilon_{n+1}|}\]If the denominator is replaced by one, then this is exactly the Euler-Maruyama scheme. If the Euler-Maruyama increment (everything in the numerator) is too large (i.e. the start of a "bad" oscillation"), then the denominator term "tames" the whole increment and prevents it from becoming too large.
Those (and similar taming or truncation schemes) are strongly convergent in \(p\)th moments even for the class of "difficult" SDEs like the one considered above. See "Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients" by Hutzenthaler and Jentzen.
Strongly convergent in \(p\)th moments.
Now consider again the EnKF iteration scheme
\[\begin{aligned}{u_{n+1}^{(i)}} = {u_n^{(i)}} &+ h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}({y} - A{u_n^{(i)}})\\ &+ \sqrt h\cdot C^\text{emp}A^T\left(h\cdot AC^\text{emp}A^T + \Gamma\right)^{-1}\Gamma^{1/2} \zeta_{n+1}^{(i)}\end{aligned}\]Now we pretend that all matrices here are real quantities in order to compare this visually to a taming scheme:
\[''\begin{aligned}{u_{n+1}^{(i)}} = {u_n^{(i)}} &+ \frac{h\cdot C^\text{emp}A^T({y} - A{u_n^{(i)}})+ \sqrt h\cdot C^\text{emp}A^T\Gamma^{1/2} \zeta_{n+1}^{(i)}}{\left(h\cdot AC^\text{emp}A^T + \Gamma\right)}\end{aligned}''\]Our paper "A strongly convergent numerical scheme from Ensemble Kalman inversion" shows that the EnKF iteration scheme (albeit in a very simplified setting) also strongly converges in \(p\)th moments to the continuous SDE
\[\mathrm d u^{(i)} = C^\text{emp}A^T\Gamma^{-1}(y-Au^{(i)}) \mathrm d t + C^\text{emp}A^T\Gamma^{-1/2}\mathrm d W^{(i)}_t.\]This is remarkable: While usually one starts with an SDE and tries to devise a numerical scheme that will approximate it, here we started with an iteration and it turns out that it is a particularly nice discretization of its own "formal" continuum limit. This means that in order to capture the dynamical properties of the iteration (which are generally harder to obtain for non-continuous dynamical systems than for continuous ones like solutions of ODEs), we can justifiably work with the continuum limit instead. This motivates the paper we will discuss next.
Stability and Convergence properties of the continuous EnKF evolution (gist of [2])
We talked about how the EnKF iteration gives rise to a continuous "EnKF evolution" to which it (supposedly, still largely open problem) converges strongly in \(p\)th moments. This is good news because this means that instead of analyzing the dynamical properties of a (stochastic) discrete iteration, which is hard,Cf. the Collatz conjecture for an elementary deterministic iteration with unknown dynamical properties. we can analyse the dynamics of the solution of an SDE, which is comparatively simple.
Recall the form of the EnKF SDE:
\[\mathrm d u^{(i)} = C^\text{emp}(u)A^T\Gamma^{-1}(y-Au^{(i)}) \mathrm d t + C^\text{emp}(u)A^T\Gamma^{-1/2}\mathrm d W^{(i)}_t~ \text{ for } i=1,\ldots,J.\]Alternatively, the following form can be more revealing:
Alternatively
Let's record a few observations about this dynamics:
- The dynamics does not leave the subspace spanned by initial ensemble of particles \(\{u^{i}_0\}_{i=1}^J.\)
Proof: \[C^\text{emp}(u) = \frac{1}{J} \sum_{i=1}^J (u^{(i)}-\bar u) \cdot (u^{(i)}-\bar u)^T.\] This means that \(\mathrm d u^{(i)} = \frac{1}{J} \sum_{i=1}^J (u^{(i)}-\bar u) \cdot [\cdots \mathrm d t + \cdots \mathrm d W^{(i)}]\) and the increments are a linear combination of the \(u^{(j)},~j=1\ldots J\) and hence the particles themselves never cease to be a linear combination of the initial ensemble. - The mapped particles \(Au^{(i)}\) have a tendency to get closer to \(y\) (albeit this is perturbed by the noise).
Proof: \[\mathrm d A u^{(i)} = AC^\text{emp}(u)A^T\Gamma^{-1} \left(\left[y-Au^{(i)}\right]\mathrm d t +\Gamma^{1/2}\mathrm d W^{(i)}_t\right)\] As the preconditioner \(AC^\text{emp}(u)A^T\Gamma^{-1}\) is positive definite, the drift term seeks to minimize the distance between \(y\) and \(Au^{(i)}.\)
Interactive Box: Play with the EnKF
Consider \(A = (2, 0)\), i.e. \(\mathbb R \ni y = Au + \epsilon = 2\cdot u_1 + \epsilon.\)
This means that the observation operator only "sees" the first component of a two-dimensional parameter. Hence we can only hope to infer the first component.
Explanation of the figure's components- The green solid dot is the true parameter value, set constant to \((1,1).\)
- The red dashed line marks \(x_1 = \frac{1}{2}y\) where the observation \(y\) is set constant to 2.5. This dashed line is the maximum likelihood estimator of the inference problem.
- The blue dashed line marks \(x_1 = 1\), the "observable true parameter" manifold.
- The grey dashed lines mark a one (noise's) standard deviation interval around the true parameter.
- The "big" black solid dots (may be hidden behind red solid dots) are the initial ensemble.
- The "small" black solid dots show the trajectory of all particles throughout the EnKF dynamics
- You can scroll through the dynamics by using the slider abov the figure.
- The solid red dot marks the particles' positions at the specified by the slider.
Existence and Uniqueness
A classical way of proving existence and uniqueness of an SDE
\[\mathrm d x_t = b(t,x_t)\mathrm d t + \sigma(t,x_t)\mathrm d B_t\]is by checking the following two sufficient conditionsFrom Oksendal, "Stochastic differential equations":
has global unique solution if
\[ |b(t,x)| + |\sigma(t,x)| \leq C(1+|x|) ~\text{ (linear growth bound)} \]for some \(C\) and
\[ |b(t,x)-b(t,y)| + |\sigma(t,x)-\sigma(t,y)| \leq D|x-y| ~\text{ (global Lipschitz continuity)}.\]This (and the additional usual measurability condition) proves the existence of a unique (strong) solution for all times considered.
\(\mathrm d x = x^3 \mathrm d t+ x^2\mathrm d W\) does not fulfill these conditions and has no global solution a.s.
For example, \(\mathrm d x = x^3 \mathrm d t+ x^2\mathrm d W\) does not fulfill the linear growth condition and actually does not have a global solution because solutions tend to blow up in finite time.
Problem is, a lot of interesting SDEs do not fall in this category: \(\mathrm d x = -x^3 \mathrm d t + x^2\mathrm d W\) does not fulfill the conditionsFor the same reason as the previous counterexample, as the linear growth condition does not take into account the sign of the drift term. either, but it actually has unique global solutions! This is due to the fact that the dynamics forces trajectories towards 0.
\(\mathrm d x = -x^3 \mathrm d t + x^2\mathrm d W\) does not fulfill these conditions but has a global solution a.s.
The SDE we consider in the SDE setting,
Similarly, EnKF SDE:
\[\mathrm d u^{(i)} = C^\text{emp}(u)A^T\Gamma^{-1}(y-Au^{(i)}) \mathrm d t + C^\text{emp}(u)A^T\Gamma^{-1/2}\mathrm d W^{(i)}_t,\]actually is very similarFor a hand-wavy argument: Count all occurences of \(u\) (for example two in \(C^\text{emp}\) and ignore all other terms. to \(\mathrm d u = -u^3 \mathrm d t + u^2\mathrm d W\)
The general idea on how to prove existence and uniqueness in this setting is the following:
How do we prove existence and uniqueness for such "dissipative" SDEs?
- Define regularized SDEs by cutting off the problematic drift and diffusion terms outside of growing regions. This gives a sequence of SDEs which now have unique solutions which are local solutions of the original SDE (up to some stochastic exit time from the unchanged region).
- Show that the regularized SDEs do not blow up in finite time and that for any finite time we can use the solution of a regularized SDEs solution as the solution of the original SDE.
- This constructs a unique solution of the original SDE.
In order to prove that the regularized SDEs do not blow up in finite time, we use a Lyapunov function whose growth we can bound. For details, see Khasminskii, "Stochastic stability of differential equations".
- Regularize drift and diffusion terms (outside of large ball); regularized SDEs have unique solution
- Show that the regularized SDEs solve original SDE up to large stopping times (using Lyapunov functions)
- Construct original SDE's solution from regularized solutions
In order to study the dynamics of the EnKF evolution, we separate it into two natural parts:
Separate dynamics into
- Ensemble collapse: Convergence of the particles to their joint average.
- Mean convergence: Limit behavior of the average to some value (hopefully the true parameter).
These two types of convergence are at the same time "orthogonal" to each other and strongly intertwined.
In order to separate ensemble collapse from mean convergence, we define new quantities: Instead of considering the particles \(u^{(i)}\) individually, we define
New quantities:
\[\begin{aligned} \bar u &= \frac{1}{J} \sum_{i=1}^J u^{(i)} \quad \text{(ensemble mean)} \\ e^{(i)} &= u^{(i)} - \bar u \quad \text{(particle $i$'s difference vector to the ensemble mean)} \end{aligned} \]Of course, we now have one redundant \(e\) but this does not bother us.
We will see that we can quantify the ensemble collapse just by considering the \(e^{(i)},\) in the following sense:
Ensemble collapse \(\quad\Leftrightarrow\quad\) \(e^{(i)}\to 0 ~ \) for all \(~i.\)Similarly,
Mean convergence \(\quad\Leftrightarrow\quad\) \(\bar u \to u^\dagger.\)where \(u^\dagger\) is ground truth.
where \(u^\dagger\) is some "good" parameter value like the "true" one.
Actually, we will prove convergence to the true parameter value slightly differently, by defining the residuals \(r^{(j)} = u^{(j)} - u^\dagger\) as the difference vectors from each particle to ground truth. We then showalbeit only in observation space and with a "variance inflation" modification that all residuals converge to 0. This image of decomposition into \(e^{(j)}, \bar u\) is a useful mental tool nonetheless.
We start with deriving the SDE for the mean \(\bar u\) and the mean-neutral shifted particles \(e^{(i)}:\)
SDEs:
\[\begin{aligned} \mathrm d\bar u &= C^\text{emp}(e)A^T\Gamma^{-1}(y-A\bar u) \mathrm d t + C^\text{emp}(e)A^T\Gamma^{-1/2}\mathrm d \bar W_t\\ \mathrm d e^{(i)}&= -C^\text{emp}(e)A^T\Gamma^{-1}e^{(i)} \mathrm d t + C^\text{emp}(e)A^T\Gamma^{-1/2}\mathrm d (W^{(i)}_t-\bar W_t) \end{aligned}\]Note that \(e^{(i)}\) is independent from \(\bar u\), but not vice versa.
Interestingly, the dynamics of the \(e^{(i)}\) is independent from \(\bar u\) (and even from each individual \(u^{(i)}\)), but not vice versa! This means that we can prove ensemble collapse independently from "where we go", but ensemble collapse plays an important role for mean convergence. We will discuss this issue more in a later section.
There is another issue: The forward operator \(A\) is usually not one-to-one, so the observation \(y\) (which is the only information about the true parameter \(u\) that we have) lacks information in multiple dimensions of \(u\). Hence we can only expect convergence (both ensemble collapse and mean convergence) in those "informed" dimensions.
For this reason we define the mapped quantities
Can only expect convergence in observation space, hence define
\[ \begin{aligned} \mathfrak e^{(i)} &= \Gamma^{-1/2}A e^{(i)}\\ \bar \mathfrak u &= \Gamma^{-1/2}A \bar u. \end{aligned} \]The SDEs for those quantities are even more simple:
\[ \begin{aligned} \mathrm d \mathfrak e^{(i)} &= -C^\text{emp}(\mathfrak e) \mathfrak e^{(i)} \mathrm d t + C^\text{emp}(\mathfrak e)\mathrm d (W^{(i)}_t-\bar W_t)\\ \mathrm d \bar \mathfrak u &= -C^\text{emp}(\mathfrak e) (\mathfrak u^\dagger - \bar \mathfrak u) \mathrm d t + C^\text{emp}(\mathfrak e)\mathrm d \bar W_t\end{aligned}\]Ensemble collapse
Want to prove \(e^{(i)}\to 0.\) Which convergence?
Ensemble collapse is solely governed by the dynamics of the \(e^{(i)}.\) As we discussed, ensemble collapse means that all \(e^{(i)}\) (or rather, all observation space quantities \(\mathfrak e^{(i)}\)) converge to zero. As those are probabilistic quantities, "convergence" is ambiguous. We could have convergence in distribution, convergence in probability, almost sure convergence and convergence in \(p\)th moments. Those different notions of convergence depend on each other in the following way:
We see that by proving both almost sure convergence and convergence in \(p\)th moments, we get all other types of stochastic convergence. This is what we will do.
Without further ado, here is what we can prove:
The ensemble does not explode (in \(p\)th moments): \[ \begin{aligned} \mathbb E \left[\frac{1}{J} \sum_{j=1}^J | e_t^{(j)}|^p \right] \end{aligned} \text{ is monotonically decreasing}\]The ensemble collapses in observation space (in \(p\)th moments): \[ \begin{aligned} \mathbb E \left[\frac{1}{J} \sum_{j=1}^J |\mathfrak e_t^{(j)}|^p \right] \in O(t^{-p/2}) \end{aligned} \] for all \(p\in(2, \frac{J+3}{2}).\) Larger ensemble \(\Rightarrow\) more regular ensemble collapse. This means that larger ensembles (i.e. \(J\) is big) lead to more regular ensemble collapse (i.e. possible values of \(p\) are higher).The ensemble collapses in observation space (almost surely): All \(\mathfrak e_t^{(j)}\) converge to \(0\) with asymptotic rate \(t^{-1/2-\epsilon}\) where \(0<\epsilon<\frac{1}{2}\) arbitrary
The proofs are straightforward calculations with the SDEs of all quantities involved and a few justifications for things like "why is this integral a martingale" etc. Almost all proofs involve some kind of Lyapunov function like \(V = \frac{1}{J} \sum_{j=1}^J | e_t^{(j)}|^p.\) The idea is that we can use a one-dimensional quantity (the function \(V\)) which controls ensemble collapse: If \(V\to 0\)(in which probabilistic sense needs to be clarified), then all \(e_t^{(j)}\) converge to \(0\) and ensemble collapse needs to happen as well.
Proof method: Lyapunov functions like \(V = \frac{1}{J} \sum_{j=1}^J | e_t^{(j)}|^p.\)
Issues with mean convergence and ensemble collapse
We mentioned that ensemble collapse and mean convergence are tightly knit. Actually, the following holds true:
Why is that?
In order to understand this, we investigate the relationship between the dynamics of the \(\mathfrak e^{(j)}\) and \(\bar\mathfrak u = \Gamma^{-1/2}A \bar u.\)
Recall that \(\mathfrak e^{(j)} = \Gamma^{-1/2}A\cdot(u^{(j)}-\bar u)\) are the mean-shifted particles mapped into observation space. Their SDEs are:
\[ \begin{aligned} \mathrm d \mathfrak e^{(i)} &= -C^\text{emp}(\mathfrak e) \mathfrak e^{(i)} \mathrm d t &+& C^\text{emp}(\mathfrak e)\mathrm d (W^{(i)}_t-\bar W_t) \\ \mathrm d \bar \mathfrak u &= C^\text{emp}(\mathfrak e) \left(\Gamma^{-1/2}y - \bar \mathfrak u\right) \mathrm d t &+& C^\text{emp}(\mathfrak e)\mathrm d \bar W_t \end{aligned}\]The basic problem is this:
- If the ensemble collapses, all quantities \(\mathfrak e^{(i)}\) become small.
- Because of this, \(C^\text{emp}(\mathfrak e)\) becomes smaller (as a matrix, i.e. eigenvalues degenerate).
- Hence, all increment terms in the dynamics of \( \bar \mathfrak u\) start to disappear and the dynamics stops.
Now there is a gap in this argument: If the preconditioner matrix \(C^\text{emp}(\mathfrak e)\) becomes smaller, then shouldn't also the ensemble collapse governed by the \(\mathfrak e^{(j)}\) slow down, thus giving the dynamics of \( \bar \mathfrak u\) "just enough time" to converge after all?
This is, in fact, a subtle matter which we will illustrate with the following two-dimensional toy model.
A hopefully illustrative example
Consider a 2d system of ODEs:
\[\begin{aligned} w' &= -w^3\\ z' &= -w^2 \cdot z\end{aligned}\]In this model, \(w\) plays the role of the \(\mathfrak e\) (the ensemble collapse quantity) and \(z\) is supposed to symbolize \(\bar\mathfrak u\) (the quantity which ought to converge, in this case to \(0.\))
Note that \(w\) tends to \(0\) and thereby slows down the convergence of \(z\) to 0.
- If \(w\) is large, then \(z\) approaches \(0\) quickly.
- If \(w = 0\), the dynamics of \(z\) stops, as \(z' = 0\).
Luckily, in this case we can solve this ODE system analytically, with
\[\begin{aligned}w(t) &= \pm \frac{1}{\sqrt{C + 2t}} \text{ (sign depends on $w(0)$)} \\ z(t) &= \frac{D}{\sqrt{C+2t} } \end{aligned}\]where the constants \(C,D\) are fixed by the initial values \(w(0),z(0).\)
Hence, we see that all trajectories of \(w,z\) converge to 0
In the language of dynamical systems, this means that \(0\) is
This fact is also the reason why we cannot use Lyapunov theory here: There is no local neighborhood of the fix point \((0,0)\) such that it is asymptotically stable in this neighborhood.
We cannot use Lyapunov functions here: \((0,0)\) is asymptotically for no local neighborhood.
The following interactive figure plots the phase portrait for a few initial conditions on the unit circle. Note how those who are close to \(w=0\) have slower convergence to the origin.
A rate of \(t^{-1}\) is "only just" alright
In order to generalize the discussion, we abstract away the "collapse" ODE for a moment and consider ODEs of the form
\[z' = -\gamma(t) \cdot z(t)\]with \(\gamma(t) \to 0.\) Does \(z\to 0\)? Depends on rate of \(\gamma\)!
with \(\gamma(t) \to 0.\) If the factor \(\gamma\) converges to \(0\) too quickly, the dynamics of \(z\) will stop. If \(\gamma\) converges to \(0\) in a "mild" way, \(z\) will converge to 0.
It turns out that \(\gamma(t) \in O(t^{-1})\) is exactly the border case where convergence still happens. Note that this is the 2d ODE example from above (\(w\in O(t^{-1/2})\) and \(z'\) had a prefactor of \(w^2\)). Here, all trajectories converged (all but the one where the initial conditions were such that actually \(\gamma(t) = 0\) from the start).
For \(\gamma(t)\) converging faster than \(t^{-1}\) to \(0\), there is an open interval containing \(0\) such that the trajectory of \(z\) never enters this interval. Consider, for example, \[z' = -\frac{1}{t^2}\cdot z(t).\] This has unique solution \(Ce^{1/t}\) which never enters the interval \((-C, C).\)
Compare \(z' = -\frac{1}{t}\cdot z(t)\) with \(z' = -\frac{1}{t^2}\cdot z(t)\)
In the toy example we "learned" (this is just heuristics, of course) that for a dynamical system of the form \[z' = -\gamma(t) \cdot z(t),\]
- we expect \(z\to 0\) for \(\gamma \to 0\) slower than or equal to \(\frac{1}{t}\)
- we expect \(z\not\to 0\) for \(\gamma \to 0\) faster than \(\frac{1}{t}\)
Now we recall the dynamics which governs mean convergence,
\[ \mathrm d \bar \mathfrak u = C^\text{emp}(\mathfrak e) \left(\Gamma^{-1/2}y - \bar \mathfrak u\right) \mathrm d t + C^\text{emp}(\mathfrak e)\mathrm d \bar W_t\]In order to make the relevant terms more transparent, we define \(z:= \bar \mathfrak u - \Gamma^{-1/2}y.\) Then
Define \(z:= \bar \mathfrak u - \Gamma^{-1/2}y.\)
\[\mathrm d z = - C^\text{emp}(\mathfrak e) z + C^\text{emp}(\mathfrak e)\mathrm d \bar W_t\]Now convergence of \(z\) to \(0\) means mean convergence to ground truth.
Roughly (ignoring stochasticity): \(\mathfrak e^{(i)}\in O(t^{-1/2})\), hence \(C^\text{emp}(\mathfrak e)\in O(t^{-1})\)
We saw (see the corresponding section) that the collapse of the particles \(\mathfrak e^{(i)}\) happens with "velocity" \(t^{-1/2}.\) This means that the preconditioner \(C^\text{emp}\) has rate \(t^{-1}\)The rate of handwaving here is accelerating; What does it mean for a matrix to have a rate? How can we transfer from our 1d example to a nd matrix case? But we only want to motivate what the problem here is. The solution and its rigorous proof can be found in the paper. which is the border case for convergence. Now we claimed that for exactly this rate we still have convergence of \(z\) to \(0,\) but we cannot directly transfer this to the multi-dimensional setting and, a lot more importantly, we now have an SDE instead of an ODE. Hence we cannot (a priori) decide whether the ensemble collapse happens
- so fast that the dynamics of \(\bar \mathfrak u\) is slowed down so much that the quantity does not converge to the correct parameter or
- slowly enough that we still get convergence to ground truth.
This is exactly "border case". Do we have convergence? Obviously, everything is more complicated here than in the ODE case.
Note that the SDE's white noise terms make the analysis a lot more difficult and it is not obvious whether the effect is helping us or hindering ground truth convergence. For the same reason as in the 2d ODE example above we cannot use Lyapunov theory either: The origin cannot be an asymptotically stable fix point as there is a whole subspace of fixpoints which goes through it (which is the set \(C^\text{emp}(\mathfrak e) = 0\)) and hence we cannot find a neighborhood for which the origin is stable.
Proving mean convergence using variance inflation
As the ensemble collapse happens with such a rate that we only just cannot tell whether it is "too fast" in order for mean convergence to work or not, we modify the EnKF evolution. That way, we can actually prove convergence of the mean to the correct parameter. Or rather, convergence of the mean, mapped into observation space, to the correct parameter, mapped into observation space. Obviously we cannot expect convergence in parameter space where we might have too little information.
The modification is called
for some positive definite matrix \(B,\) a constant \(R\) and a rate \(\alpha\in (0,1).\)
With this modification we can prove convergence of the mean, using similar methods as for the proof of ensemble collapse.