Large forward maps make Julia crash

[adelinehillier]

When trying to calibrate to LESbrary simulations, the forward map typically has over half a million elements, making it impossible to store a covariance matrix for EKI. How do you guys want to handle this? We could

1. Look into using sparse matrices.
2. Store fewer time steps in observations. This would make it harder to animate the solution downstream.
3. Add a time_range attribute to ConcatenatedOutputMap that specifies which time steps to track--defaults to all time steps.
4. Delete OceanTurbulenceParameterEstimation.jl and sip guava juice from Costco while brainstorming alternative research topics.

@glwagner @navidcy What do you guys think?

[navidcy]

Oh! That's why Julia was unexpectedly quitting?!

I think I was running onto similar problem when I tried to calibrate the 3D model....

[navidcy]

Let's do 4. I don't have Costco near me though. Problem!

[navidcy]

How about 3?

Sparse matrices sounds good but also sounds bit like a mini-research endeavour?

[adelinehillier]

I just realized that EnsembleKalmanProcesses.jl can't take non-Matrix types for the covariance matrix, so I'm also leaning towards 3.

[glwagner]

I'd vote for using existing tools to solve this problem without changing source code, which is option 2, so:

i. Restrict comparison times using the times keyword argument when constructing OneDimensionalTimeSeries;
ii. When evaluating a parameterization / animating results, add more entries to times and build a new InverseProblem from the existing one (just changing out the new observations); or
iii. Add output writers to ensemble_simulation (this can be done just by appending an output writers to an existing InverseProblem.simulation)

We'll want to do iii anyways, especially when the simulations get more expensive.

[glwagner]

Another option is to add a new output map that computes some kind of loss (perhaps the easiest is just the norm of the difference between forward output and observations)?

[glwagner]

I just realized that EnsembleKalmanProcesses.jl can't take non-Matrix types for the covariance matrix, so I'm also leaning towards 3.

Eg we can't use sparse matrices? Damn the maps.

[glwagner]

Let's do 4. I don't have Costco near me though. Problem!

We can fix that

[glwagner]

We can also add a convenience utility that rebuilds an inverse problem with more times, eg something like

new_ip = with_observation_times(old_ip, new_times)

[eviatarbach]

Just a note: we implemented ensemble transform Kalman inversion in EKP: CliMA/EnsembleKalmanProcesses.jl#329

This variant of EKI not have the problem described here. In particular, the covariance matrix never has to be formed, and the time complexity will scale linearly with observation size.

[glwagner]

That sounds useful @eviatarbach. It looks like the noise covariance is used here:

ParameterEstimocean.jl/src/EnsembleKalmanInversions.jl

Lines 500 to 505 in a9dbf70

	function step_parameters(X, G, y, Γy, process; Δt=1.0)
	ekp = EnsembleKalmanProcess(X, y, Γy, process; Δt)
	update_ensemble!(ekp, G)
	return get_u_final(ekp)
	end

Can you show us how we would use ensemble transform Kalman inversion instead?

Unfortunately, looking at the code, it seems that when using any adaptive stepping schemes the use of the noise covariance during the EKI update is a bit tangled up with the time-step estimation. The code could use a little refactoring to create a separation of concerns there.

(The interface for "no adaptive timestepping" may be broken right now --- I'm not sure its tested)

[eviatarbach]

Ah yes, you still need to provide the noise covariance (in fact, its inverse); ETKI will only have better time complexity than EKI in the case that this has a simple structure, ideally diagonal. If this is satisfied, then I would try, using the latest dev version of EKP (haven't created a new release yet), to just create an EKI struct with argument ensemble_kalman_process = TransformInversion(inv_Γ) and see if it works automatically.

If this runs into problems I would be happy to help troubleshoot.

[glwagner]

Ah ok, so this might work now with user input here:

ParameterEstimocean.jl/src/EnsembleKalmanInversions.jl

Line 168 in a9dbf70

ensemble_kalman_process = Inversion(),

PS the distinction between the "process" and EnsembleKalmanProcess is a little confusing (same name for different things?) how should we think about this?

[eviatarbach]

Yes, exactly.

I agree, the naming is confusing: in one case, it refers to the struct which indicates what type of inversion to use, e.g. Inversion() or Unscented(). In the other case, it refers to the struct holding the entire history of the ensembles, etc. I think it might be good to rename the former to inversion_type or something like that.

[glwagner]

Mm yeah! We avoid names with "type" in them for julia code (too easy to confused with julia types) but something like that might work. Unscented is a kind of Inversion right? But Inversion and Unscented are different?

[eviatarbach]

Good point on the naming. Yes, Unscented is a variant of EKI, whereas Inversion is regular EKI. There is also TransformInversion, Sampler, etc.