Implement options for pseudo-time stepping

[adelinehillier]

This branch builds on glw/pseudotime and adds several options for pseudo-time stepping schemes in the PseudoSteppingSchemes module.

To do:

• Add tests

[glwagner]

Recommdnation: generalize noise_covariance to include storage of precomputed matrice, then define an interface with a few lines:

noise_covariance(eki::EnsembleKalmanInversion) = noise_covariance(eki.noise_covariance)
inverse_noise_covariance(eki::EnsembleKalmanInversion) = inverse_noise_covariance(eki.noise_covariance)
sqrt_inverse_noise_covariance(eki::EnsembleKalmanInversion) = sqrt_inverse_noise_covariance(eki.noise_covariance)


noise_covariance(Γy) = Γy
inverse_noise_covariance(Γy) = inv(Γy)
sqrt_inv_noise_covariance(Γy) = sqrt(inv(Γy))

struct PrecomputedInverseNoiseCovariance
    covariance
    inverse_covariance
    sqrt_inverse_covariance
end

PrecomputedInverseNoiseCovariance(Γy) = PrecomputedInverseNoiseCovariance(Γy, inv(Γy), sqrt(inv(Γy)))

noise_covariance(pinc::PrecomputedInverseNoiseCovariance) = pinc.covariance
inverse_noise_covariance(pinc::PrecomputedInverseNoiseCovariance) = pinc.inverse_covariance
sqrt_inv_noise_covariance(pinc::PrecomputedInverseNoiseCovariance) = pinc.sqrt_inverse_covariance

[adelinehillier]

Thanks for the suggestion! I had implemented precomputed_arrays before seeing your reply to #245. My concern is that it might confuse users who probably have no use for inverse_covariance or sqrt_inverse_covariance. Also, I have since decided to add Γθ, μθ, and inv_sqrt_Γθ, which all correspond to EKI prior, to the list of precomputed arrays. We could add a PrecomputedPriors object to store this information, but I don't want to drown users or developers in a jungle of made-up names and abstractions.

[glwagner]

Certainly, there's no way around using words and/or symbols to describe things. I prefer names to symbols both as a developer and user, but if others prefer Γθ to noise_covariance then we should go with that.

I agree that we often won't have use for inverse_covariance. So I suggest that when we don't need it, noise_covariance is just Matrix, UniformScaling, etc. When we do need it, noise_covariance is the more complicated object PrecomputedInverseNoiseCovariance. This will solve the problems you've brought up:

1. Users aren't drowned in excess names, because the objects they build only contain the information they need.
2. Developers aren't drowned in excess names, because each "noise covariance" has its own implementation. Thus you only need to read the part of the code that's specific to your problem.

[glwagner]

What are "precomputed_arrays"?

[adelinehillier]

I agree that we should move the fixes to a new PR.

[adelinehillier]

It's all developer stuff; it doesn't affect any of the user scripts.

[glwagner]

Can you clarify when this is used? Does this pertain to adaptive stepping (eg to a particular scheme?) or does this pertain to information stored in the iteration summaries?

[adelinehillier]

tikhonov determines whether Tikhonov regularization (corresponding to the prior misfit term $\frac{1}{2}{\left| {\Gamma_\theta}^{-\frac{1}{2}}(\theta - \mu_{\theta}) \right|}^2$ in the EKI objective function $\frac{1}{2}{\left| {\Gamma_y}^{-\frac{1}{2}}(y - G(\theta)) \right|}^2 + \frac{1}{2}{\left| {\Gamma_\theta}^{-\frac{1}{2}}(\theta - \mu_{\theta}) \right|}^2$) happens. If it is false then EKI only considers the data misfit $\frac{1}{2}{\left| {\Gamma_y}^{-\frac{1}{2}}(y - G(\theta)) \right|}^2$. Certain arrays have to be augmented if one wishes to incorporate the Tikhonov term, hence these lines in PseudoSteppingSchemes.jl:

observations(eki) = eki.tikhonov ? eki.precomputed_arrays[:y_augmented] : eki.mapped_observations
obs_noise_covariance(eki) = eki.tikhonov ? eki.precomputed_arrays[:Σ] : eki.noise_covariance
inv_obs_noise_covariance(eki) = eki.tikhonov ? eki.precomputed_arrays[:inv_Σ] : eki.precomputed_arrays[:inv_Γy]

Ideally these augmented arrays should only be computed once, hence the need to store some of the precomputed_arrays.

[glwagner]

Can you add a comment that documents what each of these properties mean? When we have a type that represents this data, we can add that to the docstring.

[adelinehillier]

If I were to implement an object with non-symbol attributes to store all the information in precomputed_arrays, it would look something like this:

struct PrecomputedArrays{M, T}
        inverse_observation_noise_covariance :: M
        inverse_sqrt_observation_noise_covariance :: M
        prior_covariance :: M
        inverse_sqrt_prior_covariance :: M
        prior_means :: T
        augmented_mapped_observations :: T
        augmented_observation_noise_mean :: M
        augmented_observation_noise_covariance :: M
        inverse_augmented_observation_noise_covariance :: M
        inverse_sqrt_augmented_observation_noise_covariance :: M
end

[glwagner]

Can you also report which pseudo-time-stepping methods require which data?

[glwagner]

Alignment is off here

[glwagner]

Alignment is off in this docstring

[glwagner]

We just need a bit more documentation so that other developers can contribute, modify, and refactor this code in the future. Ideally, we don't use dictionaries for important properties. Something that isn't clear to me is whether precomputed_arrays and precomputed_augmented_arrays are needed for all pseudostepping schemes, or just some. It would seem that they can't always be needed, since they didn't exist previously.