Merge #74

74: Substitute unnecessary matrix inversion, use linear solve instead. r=ilopezgp a=ilopezgp

This PR solves issue #73 . Inverting the observational covariance matrix is not necessary, a linear solve can be used instead. The linear solve is faster and more accurate, which can be important when the observational covariance matrix is large or ill-conditioned. 

Co-authored-by: ilopezgp <ilopezgp@gmail.com>

src/MCMC.jl

@@ -33,8 +33,6 @@ struct MCMCObj{FT<:AbstractFloat, IT<:Int}
     obs_sample::Vector{FT}
     "covariance of the observational noise"
     obs_noise_cov::Array{FT, 2}
-    "inverse of obs_noise_cov"
-    obs_noise_covinv::Array{FT, 2}
     "array of length N_parameters with the parameters' prior distributions"
     prior::Array{Prior, 1}
     "MCMC step size"
@@ -96,14 +94,12 @@ function MCMCObj(
     param = param_init_copy
     log_posterior = [nothing]
     iter = [1]
-    obs_noise_covinv = inv(obs_noise_cov)
     accept = [0]
     if algtype != "rwm"
         error("only random walk metropolis 'rwm' is implemented so far")
     end
     MCMCObj{FT,IT}(obs_sample,
                    obs_noise_cov,
-                   obs_noise_covinv,
                    priors,
                    [step],
                    burnin,
@@ -167,12 +163,11 @@ function log_likelihood(mcmc::MCMCObj{FT},
     log_rho = FT[0]
     if gvar == nothing
         diff = g - mcmc.obs_sample
-        log_rho[1] = -FT(0.5) * diff' * mcmc.obs_noise_covinv * diff
+        log_rho[1] = -FT(0.5) * diff' * (mcmc.obs_noise_cov \ diff)
     else
-        gcov_inv = inv(Diagonal(gvar))
         log_gpfidelity = -FT(0.5) * log(det(Diagonal(gvar))) # = -0.5 * sum(log.(gvar))
         diff = g - mcmc.obs_sample
-        log_rho[1] = -FT(0.5) * diff' * gcov_inv * diff + log_gpfidelity
+        log_rho[1] = -FT(0.5) * diff' * (Diagonal(gvar) \ diff) + log_gpfidelity
     end
     return log_rho[1]
 end

test/MCMC/runtests.jl

@@ -69,7 +69,6 @@ using CalibrateEmulateSample.GPEmulator
     # 1.0/sqrt(0.05) * obs_sample ≈ 4.472
     @test mcmc.obs_sample ≈ [4.472] atol=1e-2
     @test mcmc.obs_noise_cov == σ2_y
-    @test mcmc.obs_noise_covinv == inv(σ2_y)
     @test mcmc.burnin == burnin
     @test mcmc.algtype == mcmc_alg
     @test mcmc.iter[1] == max_iter + 1