PrISM.jl

RTO-TKO

RTO-TKO is a stochastic framework introduced in the electromagnetic geophysical community by Blatter et al., 2022 (a and b)

Similar to fixed discretization scheme, the grid sizes do not change. RTO-TKO tries to obtain the samples of the posterior distribution by unrolling the optimization for the misfit function, that is, it perturbs the model space obtained after one iteration before optimizing for the next. It does this iteratively for both the values of the model space and the regularization coefficient in alternative steps.

Therefore, if we have

m = [m_{1}, m_{2}, m_{3}, . . ., m_{N}]

and

d = F (m)

Then for $C_{d}$ , the inverse of data covariance matrix, and $C_{m}$ the inverse of model covariance matrix, we want to explore uncertainty for the following misfit function:

J (m) = [F (m) - d]^{T} C_{d} [F (m) - d] + μ m^{T} C_{m} m

where $μ$ is the regularization weight and $L$ is the derivative matrix. Probabilistically, the above equation implies that the a priori distribution of $m$ is $N (0, C_{m})$ . RTO-TKO explores the uncertainty in $m$ , as well as the $μ$ -space instead of fixing it, and gives a family of models that fit the data.

The algorithm was proposed for $C_{m}$ constructed with $L^{'} L$ , where $L$ is the discrete derivative matrix. The algorithm works as:

\begin{aligned} Solving for \\ J (m) = [F (m) - d]^{T} C_{d} [F (m) - d] + μ (L m)^{T} (L m) \\ 1) Solve for m_{i + 1} \\ Sample \tilde{d} \sim N (d, C_{d}) and \tilde{m} \sim N (0, \frac{1}{μ} (L^{T} L)) \\ Solve \\ m_{i + 1} = \argmin_{m_{i + 1}} [F (m_{i + 1}) - \tilde{d}]^{T} C_{d} [F (m_{i + 1}) - \tilde{d}] + μ_{i} [L (m_{i + 1} - \tilde{m})]^{T} [L (m_{i + 1} - \tilde{m})] \\ 2) Solve for μ_{i + 1} \\ Sample \tilde{d} \sim N (d, C_{d}) \\ Solve \\ μ_{i + 1} = \argmin_{μ_{i + 1}} [F (m_{i + 1}) - \tilde{d}]^{T} C_{d} [F (m_{i + 1}) - \tilde{d}] - l o g (p (μ_{i + 1})) \end{aligned}

Note

Usually, the prior of $μ$ is a uniform distribution and we do not have to compute the corresponding log pdf term
Implementing the above from scratch might not be trivial because of $L^{'} L$ being non-invertible, and we do the optimization in the domain defined by $ξ = \sqrt{μ} L m$ . Such a variable will then have a standard normal distribution $N (0, I)$ when $m \sim N (0, \frac{1}{μ} (L^{T} L))$

In the following section, we demonstrate RTO-TKO for a 6-layered earth, including the half-space being imaged using the MT method. The prior distribution is composed of 26 layers, with the initial model being a half-space.

Demo

Let's create a synthetic dataset first, with 10% error floors:

julia

m_test = MTModel(log10.([100.0, 10.0, 1000.0]), [1e3, 1e3])
f = 10 .^ range(-4; stop=1, length=25)
ω = vec(2π .* f)

r_obs = forward(m_test, ω)

err_phi = asin(0.02) * 180 / π .* ones(length(ω))
err_appres = 0.1 * r_obs.ρₐ
err_resp = MTResponse(err_appres, err_phi)

We don't need to construct the a priori the same way as before. Instead, we need to define the optimiser we'd use. Here, we use Occam, and an initial model composed of half-space of 100 $Ω m$ .

julia

z = collect(0:100:2.5e3)
h = diff(z)
m_rto = MTModel(2 .* ones(length(z)), vec(h))

n_samples = 100

r_cache = rto_cache(
    m_rto, [1e-2, 1e4], Occam(), n_samples, n_samples, 1.0, [:ρₐ, :ϕ], false)

rto_chain = stochastic_inverse(
    r_obs, err_resp, ω, r_cache; model_trans_utils=(; m=sigmoid_tf))

Chains MCMC chain (99×27×1 reshape(adjoint(::Matrix{Float64}), 99, 27, 1) with eltype Float64):

Iterations        = 1:1:99
Number of chains  = 1
Samples per chain = 99
parameters        = m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14], m[15], m[16], m[17], m[18], m[19], m[20], m[21], m[22], m[23], m[24], m[25], m[26], m[27]

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

        m[1]    1.9974    0.2860    0.1114     8.0542    15.4299    1.1028     ⋯
        m[2]    1.9489    0.3051    0.1605     3.6712    15.7316    1.2158     ⋯
        m[3]    1.8874    0.3385    0.1742     3.6710    15.9556    1.2107     ⋯
        m[4]    1.8479    0.3163    0.1677     3.6193    14.7835    1.2173     ⋯
        m[5]    1.7780    0.3141    0.1758     3.2150    72.1426    1.2531     ⋯
        m[6]    1.6959    0.3624    0.1766     3.8941    55.7318    1.2175     ⋯
        m[7]    1.5703    0.3186    0.1520     4.6473    17.2789    1.1703     ⋯
        m[8]    1.5049    0.2909    0.1826     2.4493    19.0725    1.3915     ⋯
        m[9]    1.3807    0.2762    0.1144     7.0299    31.4173    1.1190     ⋯
       m[10]    1.2876    0.2787    0.0368    58.8741    45.6792    1.0098     ⋯
       m[11]    1.1714    0.2886    0.0430    43.7584   110.1771    1.0283     ⋯
       m[12]    1.1086    0.2738    0.0224   152.3711    98.9271    1.0286     ⋯
       m[13]    1.0476    0.2707    0.0356    63.6664    89.6049    1.0224     ⋯
       m[14]    1.0433    0.3161    0.0811    17.1313    49.8052    1.0536     ⋯
       m[15]    1.0537    0.2821    0.0805    17.4127   102.2516    1.0774     ⋯
       m[16]    1.1517    0.2722    0.1289     4.3633    69.5186    1.1703     ⋯
       m[17]    1.2825    0.3242    0.1408     5.1783    29.5006    1.1541     ⋯
      ⋮           ⋮         ⋮         ⋮         ⋮          ⋮          ⋮        ⋱
                                                    1 column and 10 rows omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

        m[1]    1.2839    1.8869    2.0236    2.1575    2.4350
        m[2]    1.2201    1.7603    1.9803    2.1750    2.4502
        m[3]    1.1764    1.7139    1.9160    2.1140    2.3926
        m[4]    1.0790    1.6840    1.8561    2.0544    2.4202
        m[5]    1.1185    1.5892    1.7632    1.9663    2.4354
        m[6]    0.9188    1.5089    1.6865    1.9067    2.3553
        m[7]    0.9183    1.3804    1.5795    1.8044    2.0896
        m[8]    0.9324    1.3074    1.5192    1.6818    2.0197
        m[9]    0.7679    1.2082    1.4095    1.5472    1.8932
       m[10]    0.7481    1.1236    1.2856    1.4831    1.8315
       m[11]    0.6199    0.9918    1.1335    1.3493    1.7272
       m[12]    0.6061    0.9306    1.0600    1.3063    1.6223
       m[13]    0.6657    0.8478    0.9869    1.2311    1.6789
       m[14]    0.4823    0.8206    1.0032    1.2200    1.6967
       m[15]    0.6201    0.8428    1.0186    1.2570    1.6560
       m[16]    0.6697    0.9420    1.1218    1.3083    1.7481
       m[17]    0.7629    1.0297    1.2505    1.5236    1.9400
      ⋮           ⋮         ⋮         ⋮         ⋮         ⋮
                                                 10 rows omitted

Since RTO-TKO also samples the regularization coefficient along with model parameters, we exclude it to obtain another chain as:

julia

mt_chain = Chains((rto_chain.value.data[:, 1:(end - 1), :]), [Symbol("m[$i]")
                                                              for i in 1:length(z)])

Chains MCMC chain (99×26×1 Array{Float64, 3}):

Iterations        = 1:1:99
Number of chains  = 1
Samples per chain = 99
parameters        = m[1], m[2], m[3], m[4], m[5], m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14], m[15], m[16], m[17], m[18], m[19], m[20], m[21], m[22], m[23], m[24], m[25], m[26]

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

        m[1]    1.9974    0.2860    0.1114     8.0542    15.4299    1.1028     ⋯
        m[2]    1.9489    0.3051    0.1605     3.6712    15.7316    1.2158     ⋯
        m[3]    1.8874    0.3385    0.1742     3.6710    15.9556    1.2107     ⋯
        m[4]    1.8479    0.3163    0.1677     3.6193    14.7835    1.2173     ⋯
        m[5]    1.7780    0.3141    0.1758     3.2150    72.1426    1.2531     ⋯
        m[6]    1.6959    0.3624    0.1766     3.8941    55.7318    1.2175     ⋯
        m[7]    1.5703    0.3186    0.1520     4.6473    17.2789    1.1703     ⋯
        m[8]    1.5049    0.2909    0.1826     2.4493    19.0725    1.3915     ⋯
        m[9]    1.3807    0.2762    0.1144     7.0299    31.4173    1.1190     ⋯
       m[10]    1.2876    0.2787    0.0368    58.8741    45.6792    1.0098     ⋯
       m[11]    1.1714    0.2886    0.0430    43.7584   110.1771    1.0283     ⋯
       m[12]    1.1086    0.2738    0.0224   152.3711    98.9271    1.0286     ⋯
       m[13]    1.0476    0.2707    0.0356    63.6664    89.6049    1.0224     ⋯
       m[14]    1.0433    0.3161    0.0811    17.1313    49.8052    1.0536     ⋯
       m[15]    1.0537    0.2821    0.0805    17.4127   102.2516    1.0774     ⋯
       m[16]    1.1517    0.2722    0.1289     4.3633    69.5186    1.1703     ⋯
       m[17]    1.2825    0.3242    0.1408     5.1783    29.5006    1.1541     ⋯
      ⋮           ⋮         ⋮         ⋮         ⋮          ⋮          ⋮        ⋱
                                                     1 column and 9 rows omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

        m[1]    1.2839    1.8869    2.0236    2.1575    2.4350
        m[2]    1.2201    1.7603    1.9803    2.1750    2.4502
        m[3]    1.1764    1.7139    1.9160    2.1140    2.3926
        m[4]    1.0790    1.6840    1.8561    2.0544    2.4202
        m[5]    1.1185    1.5892    1.7632    1.9663    2.4354
        m[6]    0.9188    1.5089    1.6865    1.9067    2.3553
        m[7]    0.9183    1.3804    1.5795    1.8044    2.0896
        m[8]    0.9324    1.3074    1.5192    1.6818    2.0197
        m[9]    0.7679    1.2082    1.4095    1.5472    1.8932
       m[10]    0.7481    1.1236    1.2856    1.4831    1.8315
       m[11]    0.6199    0.9918    1.1335    1.3493    1.7272
       m[12]    0.6061    0.9306    1.0600    1.3063    1.6223
       m[13]    0.6657    0.8478    0.9869    1.2311    1.6789
       m[14]    0.4823    0.8206    1.0032    1.2200    1.6967
       m[15]    0.6201    0.8428    1.0186    1.2570    1.6560
       m[16]    0.6697    0.9420    1.1218    1.3083    1.7481
       m[17]    0.7629    1.0297    1.2505    1.5236    1.9400
      ⋮           ⋮         ⋮         ⋮         ⋮         ⋮
                                                  9 rows omitted

Note that the chain contains fewer samples than we had asked for. This is because a few unstable samples were filtered out. The obtained mt_chain contains the a posteriori distributions that can be saved using JLD2.jl.

julia

using JLD2
JLD2.@save "file_path.jld2" mt_chain

Since RTO-TKO is closely related to Occam, we also include occam results in our figures below. Also, note that we did not require an a priori distribution to obtain samples. However, we would need to create the same to plot our posterior samples. This can just be something that encapsulates the whole prior space (or an envelope of the same), e.g., in the case of magnetotelluric imaging, we know that the resistivity values will always be in [-2, 5] on the log-scale.

julia

modelD = MTModelDistribution(Product([Uniform(-1.0, 5.0) for i in eachindex(z)]), vec(h))

Code for this figure

julia

fig = Figure()
ax = Axis(fig[1, 1])
hm = get_kde_image!(ax, mt_chain, modelD; kde_transformation_fn=log10,
    colormap=:binary, colorrange=(-3.0, 0.0), trans_utils=(m=no_tf, h=no_tf))
Colorbar(fig[1, 2], hm; label="log pdf")

mean_kws = (; color=:seagreen3, linewidth=2)
std_kws = (; color=:red, linewidth=1.5)
# get_mean_std_image!(ax, mt_chain, modelD; confidence_interval=0.99, trans_utils=(m=no_tf, h=no_tf), mean_kwargs=mean_kws,
#     std_plus_kwargs=std_kws, std_minus_kwargs=std_kws)
ylims!(ax, [2500, 0])

plot_model!(ax, m_test; color=:black, linestyle=:dash, label="true", linewidth=2)
Legend(fig[2, :], ax; orientation=:horizontal)

The list of models can then be obtained from chains using

julia

model_list = get_model_list(mt_chain, modelD)

We can then easily check the fit of the response curves

julia

fig = Figure()
ax1 = Axis(fig[1, 1])
ax2 = Axis(fig[1, 2])

resp_post = forward(model_list[1], ω);
for i in 1:(length(model_list) > 100 ? 100 : length(model_list))
    forward!(resp_post, model_list[i], ω)
    plot_response!([ax1, ax2], ω, resp_post; alpha=0.4, color=:gray)
end
plot_response!([ax1, ax2], ω, r_obs; errs=err_resp, plt_type=:errors, whiskerwidth=10)
plot_response!([ax1, ax2], ω, r_obs; plt_type=:scatter, label="true")
ylims!(ax1, exp10.([1.2, 3.1]))
ylims!(ax2, [12, 70])

fig

Warning

It is recommended that the samples from RTO-TKO can be filtered out to reject the samples that have a poor fit on the data.

RTO-TKO ​

Demo ​

RTO-TKO

Demo