KronLinInv

Quick overview

Kronecker-product-based linear inversion of geophysical (or other kinds of) data under Gaussian and separability assumptions. The code computes the posterior mean model and the posterior covariance matrix (or subsets of it) in an efficient manner (parallel algorithm) taking into account 3-D correlations both in the model parameters and in the observed data.

If you use this code for research or else, please cite the related paper:

Andrea Zunino, Klaus Mosegaard (2019), An efficient method to solve large linearizable inverse problems under Gaussian and separability assumptions, Computers & Geosciences. ISSN 0098-3004, https://doi.org/10.1016/j.cageo.2018.09.005.

See the above mentioned paper for a detailed description.

Installation

To install the package simple enter into the package manager mode in Julia by typing "]" at the REPL prompt and then use add, i.e.,

(v1.2) pkg> add KronLinInv

The package will be automatically downloaded from the web and installed.

Theoretical background

KronLinInv solves the linear inverse problem with Gaussian uncertainties represented by the following objective function

\[S( \mathbf{m}) = \frac{1}{2} ( \mathbf{G} \mathbf{m} - \mathbf{d}_{\sf{obs}} )^{\sf{T}} \mathbf{C}^{-1}_{\rm{D}} ( \mathbf{G} \mathbf{m} - \mathbf{d}_{\sf{obs}} ) + \frac{1}{2} ( \mathbf{m} - \mathbf{m}_{\sf{prior}} )^{\sf{T}} \mathbf{C}^{-1}_{\rm{M}} ( \mathbf{m} - \mathbf{m}_{\sf{prior}} )\]

under the following separability conditions (for a 3-way decomposition):

\[\mathbf{C}_{\rm{M}} = \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{C}_{\rm{M}}^{\rm{z}} , \quad \mathbf{C}_{\rm{D}} = \mathbf{C}_{\rm{D}}^{\rm{x}} \otimes \mathbf{C}_{\rm{D}}^{\rm{y}} \otimes \mathbf{C}_{\rm{D}}^{\rm{z}} \quad \textrm{ and } \quad \mathbf{G} = \mathbf{G}^{\rm{x}} \otimes \mathbf{G}^{\rm{y}} \otimes \mathbf{G}^{\rm{z}} \, .\]

From the above, the posterior covariance matrix is given by

\[ \mathbf{\widetilde{C}}_{\rm{M}} = \left( \mathbf{G}^{\sf{T}} \, \mathbf{C}^{-1}_{\rm{D}} \, \mathbf{G} + \mathbf{C}^{-1}_{\rm{M}} \right)^{-1}\]

and the center of posterior gaussian is

\[ \mathbf{\widetilde{m}} = \mathbf{m}_{\rm{prior}}+ \mathbf{\widetilde{C}}_{\rm{M}} \, \mathbf{G}^{\sf{T}} \, \mathbf{C}^{-1}_{\rm{D}} \left(\mathbf{d}_{\rm{obs}} - \mathbf{G} \mathbf{m}_{\rm{prior}} \right) \, .\]

KronLinInv solves the inverse problem in an efficient manner, with a very low memory imprint, suitable for large problems where many model parameters and observations are involved.

The paper describes how to obtain the solution to the above problem as shown hereafter. First the following matrices are computed

\[ \mathbf{U}_1 \mathbf{\Lambda}_1 \mathbf{U}_1^{-1} = \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \mathbf{G}^{\rm{x}}\]

\[\mathbf{U}_2 \mathbf{\Lambda}_2 \mathbf{U}_2^{-1} = \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \mathbf{G}^{\rm{y}}\]

\[\mathbf{U}_3 \mathbf{\Lambda}_3 \mathbf{U}_3^{-1} = \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \mathbf{G}^{\rm{z}} \, .\]

The posterior covariance is then expressed as

\[\mathbf{\widetilde{C}}_{\rm{M}} = \left( \mathbf{U}_1 \otimes \mathbf{U}_2 \otimes \mathbf{U}_3 \right) \big( \mathbf{I} + \mathbf{\Lambda}_1 \! \otimes \! \mathbf{\Lambda}_2 \! \otimes \! \mathbf{\Lambda}_3 \big)^{-1} \big( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} \big) \, .\]

and the posterior mean model as

\[\mathbf{\widetilde{m}} = \mathbf{m}_{\rm{prior}} + \Big[ \! \left( \mathbf{U}_1 \otimes \mathbf{U}_2 \otimes \mathbf{U}_3 \right) \big( \mathbf{I} + \mathbf{\Lambda}_1\! \otimes \! \mathbf{\Lambda}_2 \! \otimes\! \mathbf{\Lambda}_3 \big)^{-1} \\ \times \Big( \left( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \right) \! \otimes \left( \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \right) \! \\ \otimes \left( \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \right) \Big) \Big] \\ \times \Big( \mathbf{d}_{\rm{obs}} - \big( \mathbf{G}^{\rm{x}} \otimes \mathbf{G}^{\rm{y}} \otimes \mathbf{G}^{\rm{z}} \big) \, \mathbf{m}_{\rm{prior}} \Big) \, .\]

These last two formulae are those used by the KronLinInv algorithm.

Several function are exported by the module KronLinInv:

kli_calcfactors: Computes the factors necessary to solve the inverse problem
kli_posteriormean: Computes the posterior mean model using the previously computed "factors" with kli_calcfactors.
kli_blockposteriorcov: Computes a block (or all) of the posterior covariance using the previously computed "factors" with kli_calcfactors.
kli_bandpostcov: NOT YET IMPLEMENTED! Computes a band of the posterior covariance the previously computed "factors" with kli_calcfactors.

Usage examples

The input needed is represented by the set of three covariance matrices of the model parameters, the three covariances of the observed data, the three forward model operators, the observed data (a vector) and the prior model (a vector). The first thing to compute is always the set of "factors" using the function [`klicalcfactors](@ref). Finally, the posterior mean (see [kliposteriormean`](@ref KronLinInv.kliposteriormean)) and/or covariance (or part of it) can be computed (see kli_blockposteriorcov).

2D example
3D example

2D example

An example of how to use the code for 2D problems is shown in the following. Notice that the code is written for a 3D problem, however, by setting some of the matrices as identity matrices with size of 1$\times$1, a 2D problem can be solved without much overhead.

Creating a test problem

First, we create some input data to simulate a real problem.

# set the sizes of the problem
nx = 1
ny = 20
nz = 30
nxobs = 1
nyobs = 18
nzobs = 24

We then construct some covariance matrices and a forward operator. The "first" covariance matrices for model parameters ($\mathbf{C}_{\rm{M}}^{\rm{x}} \, , \mathbf{C}_{\rm{M}}^{\rm{y}} \, , \mathbf{C}_{\rm{M}}^{\rm{z}}$) and observed data ($\mathbf{C}_{\rm{D}}^{\rm{x}} \, , \mathbf{C}_{\rm{D}}^{\rm{y}} \, , \mathbf{C}_{\rm{D}}^{\rm{z}}$) are simply an identity matrix of shape 1$\times$1, since it is a 2D problem. The forward relation (forward model) is created from three operators ($\mathbf{G}^{\rm{x}} \, , \mathbf{G}^{\rm{y}} \, , \mathbf{G}^{\rm{z}}$). Remark: the function mkCovSTAT used in the following example is not part of KronLinInv.

function mkCovSTAT(sigma::Array{Float64,1},nx::Integer,ny::Integer,nz::Integer,
                   corrlength::Array{Float64,1},kind::String)
    function cgaussian(dist,corrlength)
        if maximum(dist)==0.0
            return 1.0
        else
            @assert(corrlength>0.0)
            return exp.(-(dist./corrlength).^2)
        end
    end
    function cexponential(dist,corrlength)
        if maximum(dist)==0.0
            return 1.0
        else
            @assert(corrlength>0.0)
            return exp.(-(dist./corrlength))
        end
    end

    npts = nx*ny*nz
    x = [float(i) for i=1:nx]
    y = [float(i) for i=1:ny]
    z = [float(i) for i=1:nz]
    covmat_x = zeros(nx,nx)
    covmat_y = zeros(ny,ny)
    covmat_z = zeros(nz,nz)

    if kind=="gaussian"
        calccovfun = cgaussian
    elseif kind=="exponential"
        calccovfun = cexponential
    else
        println("Error, no or wrong cov 'kind' specified")
        exit()
    end

    for i=1:nx
        covmat_x[i,:] .= sigma[1]^2 .*
            calccovfun(sqrt.((x.-x[i]).^2),corrlength[1])
    end
    for i=1:ny
        covmat_y[i,:] .= sigma[2]^2 .*
            calccovfun(sqrt.(((y.-y[i])).^2),corrlength[2])
    end
    for i=1:nz
        covmat_z[i,:] .= sigma[3]^2 .*
            calccovfun(sqrt.(((z.-z[i])).^2),corrlength[3])
    end

    return covmat_x,covmat_y,covmat_z
 end

# standard deviations
sigmaobs  = [1.0, 0.1, 0.1] # notice the 1.0 as first element (2D problem)
sigmam    = [1.0, 0.8, 0.8] # notice the 1.0 as first element (2D problem)
# correlation lengths
corlenobs = [0.0, 1.4, 1.4] # notice the 0.0 as first element (2D problem)
corlenm   = [0.0, 2.5, 2.5] # notice the 0.0 as first element (2D problem)
# create the covariance matrices on observed data
Cd1,Cd2,Cd3 = mkCovSTAT(sigmaobs,nxobs,nyobs,nzobs,corlenobs,"gaussian")
# create the covariance matrices on model parameters
Cm1,Cm2,Cm3 = mkCovSTAT(sigmam,nx,ny,nz,corlenm,"gaussian")

# forward model operator
G1 = rand(nxobs,nx) # notice that nx=1 and nxobs=1 (2D problem)
G2 = rand(nyobs,ny)
G3 = rand(nzobs,nz)

Finally, a "true/reference" model, in order to compute some synthetic "observed" data and a prior model.

# create a reference model
refmod = rand(nx*ny*nz)

# create a prior model
mprior = copy(refmod) .+ 0.3*randn(length(refmod))

# create some "observed" data
dobs = kron(G1,kron(G2,G3)) * refmod
# add some noise to make it more realistic
dobs = dobs .+ 0.02.*randn(length(dobs))

432-element Vector{Float64}:
 13.582396486944235
 12.886410417844642
 13.162398850397112
 13.206216683898653
 16.257313687529642
 11.53005096074554
 14.190278459742004
 14.02106127530625
 15.42546236371982
 15.821884025818386
  ⋮
 18.030561479393697
 13.256722166875553
 13.175051103722899
 16.0322085975692
 14.108922764534961
 15.537000064376317
 15.906462669544686
 16.88650301675155
 14.812576862547717

Now we have create a synthetic example to play with, which we can solve as shown in the following.

Solving the 2D problem

In order to solve the inverse problem using KronLinInv, we first need to compute the "factors" using the function kli_calcfactors, which takes as inputs two structs containing the covariance matrices and the forward operators.

using InverseAlgos.KronLinInv

# create the covariance matrix structure
Covs = CovMats(Cd1,Cd2,Cd3,Cm1,Cm2,Cm3)

# forward model operator
Gfwd = FwdOps(G1,G2,G3)

# calculate the required factors
klifac = kli_calcfactors(Gfwd,Covs)

Now the inverse problem can be solved. We first compute the posterior mean and then a subset of the posterior covariance.

# calculate the posterior mean model
postm = kli_posteriormean(klifac,Gfwd,mprior,dobs)

# calculate the posterior covariance
npts = nx*ny*nz
astart, aend = 1,div(npts,3) # set of rows to be computed
bstart, bend = 1,div(npts,3) # set of columns to be computed

# compute the block of posterior covariance
postC = kli_blockposteriorcov(klifac,astart,aend,bstart,bend)

posteriormean_serial(): Serial version
posteriormean(): loop 1/3, 2 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 30 of 432; ETA: 0.005 min  
posteriormean(): loop 1/3, 60 of 432; ETA: 0.002 min  
posteriormean(): loop 1/3, 90 of 432; ETA: 0.002 min  
posteriormean(): loop 1/3, 120 of 432; ETA: 0.001 min  
posteriormean(): loop 1/3, 150 of 432; ETA: 0.001 min  
posteriormean(): loop 1/3, 180 of 432; ETA: 0.001 min  
posteriormean(): loop 1/3, 210 of 432; ETA: 0.001 min  
posteriormean(): loop 1/3, 240 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 270 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 300 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 330 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 360 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 390 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 420 of 432; ETA: 0.0 min  
posteriormean(): loop 2/3, 2 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 30 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 60 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 90 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 120 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 150 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 180 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 210 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 240 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 270 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 300 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 330 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 360 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 390 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 420 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 450 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 480 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 510 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 540 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 570 of 600; ETA: 0.0 min  
posteriormean(): loop 2/3, 600 of 600; ETA: 0.0 min  
posteriormean(): loop 3/3, 2 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 30 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 60 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 90 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 120 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 150 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 180 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 210 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 240 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 270 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 300 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 330 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 360 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 390 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 420 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 450 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 480 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 510 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 540 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 570 of 600; ETA: 0.0 min   
posteriormean(): loop 3/3, 600 of 600; ETA: 0.0 min
blockpostcov_serial(): Serial version
blockpostcov():  100 of 1 to 200 
blockpostcov():  200 of 1 to 200

3D example

An example of how to use the code for 3D problems is shown in the following. It follows closely the 3D example.

Creating a test problem

First, we create some input data to simulate a real problem.

# set the sizes of the problem
nx = 7
ny = 9
nz = 7
nxobs = 6
nyobs = 8
nzobs = 9

We then construct some covariance matrices for model parameters ($\mathbf{C}_{\rm{M}}^{\rm{x}} \, , \mathbf{C}_{\rm{M}}^{\rm{y}} \, , \mathbf{C}_{\rm{M}}^{\rm{z}}$) and observed data ($\mathbf{C}_{\rm{D}}^{\rm{x}} \, , \mathbf{C}_{\rm{D}}^{\rm{y}} \, , \mathbf{C}_{\rm{D}}^{\rm{z}}$). The forward relation (forward model) is created from three operators ($\mathbf{G}^{\rm{x}} \, , \mathbf{G}^{\rm{y}} \, , \mathbf{G}^{\rm{z}}$). Remark: the function mkCovSTAT used in the following example is not part of KronLinInv.

function mkCovSTAT(sigma::Array{Float64,1},nx::Integer,ny::Integer,nz::Integer,
                   corrlength::Array{Float64,1},kind::String)
    function cgaussian(dist,corrlength)
        if maximum(dist)==0.0
            return 1.0
        else
            @assert(corrlength>0.0)
            return exp.(-(dist./corrlength).^2)
        end
    end
    function cexponential(dist,corrlength)
        if maximum(dist)==0.0
            return 1.0
        else
            @assert(corrlength>0.0)
            return exp.(-(dist./corrlength))
        end
    end

    npts = nx*ny*nz
    x = [float(i) for i=1:nx]
    y = [float(i) for i=1:ny]
    z = [float(i) for i=1:nz]
    covmat_x = zeros(nx,nx)
    covmat_y = zeros(ny,ny)
    covmat_z = zeros(nz,nz)

    if kind=="gaussian"
        calccovfun = cgaussian
    elseif kind=="exponential"
        calccovfun = cexponential
    else
        println("Error, no or wrong cov 'kind' specified")
        exit()
    end

    for i=1:nx
        covmat_x[i,:] .= sigma[1]^2 .*
            calccovfun(sqrt.((x.-x[i]).^2),corrlength[1])
    end
    for i=1:ny
        covmat_y[i,:] .= sigma[2]^2 .*
            calccovfun(sqrt.(((y.-y[i])).^2),corrlength[2])
    end
    for i=1:nz
        covmat_z[i,:] .= sigma[3]^2 .*
            calccovfun(sqrt.(((z.-z[i])).^2),corrlength[3])
    end

    return covmat_x,covmat_y,covmat_z
 end

# standard deviations
sigmaobs  = [0.1, 0.1, 0.1]
sigmam    = [0.7, 0.8, 0.8]
# correlation lengths
corlenobs = [1.3, 1.4, 1.4]
corlenm   = [2.5, 2.5, 2.5]
# create the covariance matrices on observed data
Cd1,Cd2,Cd3 = mkCovSTAT(sigmaobs,nxobs,nyobs,nzobs,corlenobs,"gaussian")
# create the covariance matrices on model parameters
Cm1,Cm2,Cm3 = mkCovSTAT(sigmam,nx,ny,nz,corlenm,"gaussian")

# Forward model operator
G1 = rand(nxobs,nx)
G2 = rand(nyobs,ny)
G3 = rand(nzobs,nz)

Finally, a "true/reference" model, in order to compute some synthetic "observed" data and a prior model.

# create a reference model
refmod = rand(nx*ny*nz)

# create a prior model
mprior = copy(refmod) .+ 0.3.*randn(length(refmod))

# create some "observed" data
dobs = kron(G1,kron(G2,G3)) * refmod
# add some noise to make it more realistic
dobs = dobs .+ 0.02.*randn(length(dobs))

432-element Vector{Float64}:
 13.186362472290332
 21.411276845283595
  9.655590587505674
 19.865652488395416
 13.686890642611504
 16.448891508856548
 15.539324187144373
 14.474353751331472
 17.758019563690176
 26.46443238247835
  ⋮
 29.50626465207921
 48.9483471122704
 21.432546072263346
 45.232099539448285
 31.795064686995502
 35.913065724329876
 35.803110142396456
 31.441289797478397
 38.570507438138605

Now we have create a synthetic example to play with, which we can solve as shown in the following.

Solving the 3D problem

using InverseAlgos.KronLinInv

# create the covariance matrix structure
Covs = CovMats(Cd1,Cd2,Cd3,Cm1,Cm2,Cm3)

# forward model operator
Gfwd = FwdOps(G1,G2,G3)

# Calculate the required factors
klifac = kli_calcfactors(Gfwd,Covs)

Now the inverse problem can be solved. We first compute the posterior mean and then a subset of the posterior covariance.

# Calculate the posterior mean model
postm = kli_posteriormean(klifac,Gfwd,mprior,dobs)

# Calculate the posterior covariance
npts = nx*ny*nz
astart, aend = 1,div(npts,3) # set of rows to be computed
bstart, bend = 1,div(npts,3) # set of columns to be computed

# compute the block of posterior covariance
postC = kli_blockposteriorcov(klifac,astart,aend,bstart,bend)

posteriormean_serial(): Serial version
posteriormean(): loop 1/3, 2 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 22 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 44 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 66 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 88 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 110 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 132 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 154 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 176 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 198 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 220 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 242 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 264 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 286 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 308 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 330 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 352 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 374 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 396 of 432; ETA: 0.0 min  
posteriormean(): loop 1/3, 418 of 432; ETA: 0.0 min  
posteriormean(): loop 2/3, 2 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 22 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 44 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 66 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 88 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 110 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 132 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 154 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 176 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 198 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 220 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 242 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 264 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 286 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 308 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 330 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 352 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 374 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 396 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 418 of 441; ETA: 0.0 min  
posteriormean(): loop 2/3, 440 of 441; ETA: 0.0 min  
posteriormean(): loop 3/3, 2 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 22 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 44 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 66 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 88 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 110 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 132 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 154 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 176 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 198 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 220 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 242 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 264 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 286 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 308 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 330 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 352 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 374 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 396 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 418 of 441; ETA: 0.0 min   
posteriormean(): loop 3/3, 440 of 441; ETA: 0.0 min
blockpostcov_serial(): Serial version
blockpostcov():  100 of 1 to 147

Public API

InverseAlgos.KronLinInv — Module

A module to perform Kronecker-product-based linear inversion of geophysical (or other kinds of) data under Gaussian and separability assumptions. The code computes the posterior mean model and the posterior covariance matrix (or subsets of it) in an efficient manner (parallel algorithm) taking into account 3-D correlations both in the model parameters and in the observed data.

source

InverseAlgos.KronLinInv.FwdOps — Type

FwdOps

A structure containing the three forward model matrices G1, G2, G3, where $\mathbf{G} = \mathbf{G_1} \otimes \mathbf{G_2} \otimes \mathbf{G_3}$

source

InverseAlgos.KronLinInv.CovMats — Type

CovMats

A structure containing the six covariance matrices Cm1, Cm2, Cm3, Cd1, Cd2, Cd3, where $\mathbf{C}_{\rm{M}} = \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{C}_{\rm{M}}^{\rm{z}}$ and $\quad \mathbf{C}_{\rm{D}} = \mathbf{C}_{\rm{D}}^{\rm{x}} \otimes \mathbf{C}_{\rm{D}}^{\rm{y}} \otimes \mathbf{C}_{\rm{D}}^{\rm{z}}$

source

InverseAlgos.KronLinInv.KLIFactors — Type

KLIFactors

A structure containing all the factors necessary to perform further calculations with KronLinInv, as, for instance, computations of the posterior mean model or the posterior covariance matrix. The structure includes:

U1, U2, U3: $\mathbf{U}_1$, $\mathbf{U}_2$, $\mathbf{U}_3$ of $F_{\sf{A}}$
invlambda: $\big( \mathbf{I} + \mathbf{\Lambda}_1\! \otimes \! \mathbf{\Lambda}_2 \! \otimes\! \mathbf{\Lambda}_3 \big)^{-1}$ of $F_{\sf{B}}$
iUCm1, iUCm2, iUCm3: $\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}}$ of $F_{\sf{C}}$
iUCmGtiCd1, iUCmGtiCd1, iUCmGtiCd1: $\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}}(\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1}$, $\mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1}$ of $F_{\sf{D}}$

source

InverseAlgos.KronLinInv.kli_calcfactors — Function

 kli_calcfactors(Gfwd::FwdOps,Covs::CovMats)

Computes the factors necessary to solve the inverse problem.

The factors are the ones to be stored to subsequently calculate posterior mean and covariance. First an eigen decomposition is performed, to get

\[ \mathbf{U}_1 \mathbf{\Lambda}_1 \mathbf{U}_1^{-1} = \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \mathbf{G}^{\rm{x}}\]

\[ \mathbf{U}_2 \mathbf{\Lambda}_2 \mathbf{U}_2^{-1} = \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \mathbf{G}^{\rm{y}}\]

\[ \mathbf{U}_3 \mathbf{\Lambda}_3 \mathbf{U}_3^{-1} = \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \mathbf{G}^{\rm{z}} \]

The principal factors involved in the computation of the posterior covariance and mean are:

\[ F_{\sf{A}} = \mathbf{U}_1 \otimes \mathbf{U}_2 \otimes \mathbf{U}_3 \]

\[ F_{\sf{B}} = \big( \mathbf{I} + \mathbf{\Lambda}_1 \! \otimes \! \mathbf{\Lambda}_2 \! \otimes \! \mathbf{\Lambda}_3 \big)^{-1} \]

\[ F_{\sf{C}} = \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} \]

\[ F_{\sf{D}} = \left( \mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1} \right) \! \otimes \left( \mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1} \right) \! \otimes \left( \mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1} \right)\]

Uses LAPACK.sygvd!(), see http://www.netlib.org/lapack/lug/node54.html. Reduces a real symmetric-definite generalized eigenvalue problem to the standard form. \n $B A z = \lambda z$ B = LLT C = LT A L z = L y

A is symmetric
B is symmetric, positive definite

Arguments

Gfwd: a FwdOps structure containing the three forward model matrices G1, G2 and G3, where $\mathbf{G} = \mathbf{G_1} \otimes \mathbf{G_2} \otimes \mathbf{G_3}$
Covs: a CovMats structure containing the six covariance matrices $\mathbf{C}_{\rm{M}} = \mathbf{C}_{\rm{M}}^{\rm{x}} \otimes \mathbf{C}_{\rm{M}}^{\rm{y}} \otimes \mathbf{C}_{\rm{M}}^{\rm{z}}$ and $\mathbf{C}_{\rm{D}} = \mathbf{C}_{\rm{D}}^{\rm{x}} \otimes \mathbf{C}_{\rm{D}}^{\rm{y}} \otimes \mathbf{C}_{\rm{D}}^{\rm{z}}$

Returns

A KLIFactors structure containing all the "factors" necessary to perform further calculations with KronLinInv, as, for instance, computations of the posterior mean model or the posterior covariance matrix. The structure includes:
- U1, U2, U3: $\mathbf{U}_1$, $\mathbf{U}_2$, $\mathbf{U}_3$ of $F_{\sf{A}}$
- invlambda: $F_{\sf{B}}$
- iUCm1, iUCm2, iUCm3: $\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}}$ of $F_{\sf{C}}$
- iUCmGtiCd1, iUCmGtiCd1, iUCmGtiCd1: $\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}}(\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1}$, $\mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1}$ of $F_{\sf{D}}$

source

InverseAlgos.KronLinInv.kli_posteriormean — Function

kli_posteriormean(klifac::KLIFactors,Gfwd::FwdOps,mprior::Array{Float64,1},
                  dobs::Array{Float64,1}; parall::Symbol=:serial)

Computes the posterior mean model for a distributed memory setup.

Arguments

klifac: a structure containing the required "factors" previously computed with the function kli_calcfactors(). It includes
- U1, U2, U3 $\mathbf{U}_1$, $\mathbf{U}_2$, $\mathbf{U}_3$ of $F_{\sf{A}}$
- diaginvlambda $F_{\sf{B}}$
- iUCmGtiCd1, iUCmGtiCd2, iUCmGtiCd3 $\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}} (\mathbf{G}^{\rm{x}})^{\sf{T}}(\mathbf{C}_{\rm{D}}^{\rm{x}})^{-1}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}} (\mathbf{G}^{\rm{y}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{y}})^{-1}$, $\mathbf{U}_3^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}} (\mathbf{G}^{\rm{z}})^{\sf{T}} (\mathbf{C}_{\rm{D}}^{\rm{z}})^{-1}$ of $F_{\sf{D}}$
FwdOps: a structure containing the three forward matrices
- G1, G2, G3 $\mathbf{G} = \mathbf{G_1} \otimes \mathbf{G_2} \otimes \mathbf{G_3}$
mprior: prior model (vector)
dobs: observed data (vector)
parall: type of parallelization: :serial, ;sharedmem, :distrmem

Returns

The posterior mean model (vector)

source

InverseAlgos.KronLinInv.kli_blockposteriorcov — Function

kli_blockpostcov(klifac::KLIFactors,astart::Int64,aend::Int64,
             bstart::Int64,bend::Int64; parall::Symbol=:serial)

Computes a block of the posterior covariance, purely serial version.

Arguments

klifac: a structure containing the required "factors" previously computed with the function kli_calcfactors(). It includes
- U1,U2,U3 $\mathbf{U}_1$, $\mathbf{U}_2$, $\mathbf{U}_3$ of $F_{\sf{A}}$
- diaginvlambda $F_{\sf{B}}$
- iUCm1, iUCm2, iUCm3 $\mathbf{U}_1^{-1} \mathbf{C}_{\rm{M}}^{\rm{x}}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{y}}$, $\mathbf{U}_2^{-1} \mathbf{C}_{\rm{M}}^{\rm{z}}$ of $F_{\sf{C}}$
Gfwd: a structure containing the three forward model matrices G1,G2,G3, where $\mathbf{G} = \mathbf{G_1} \otimes \mathbf{G_2} \otimes \mathbf{G_3}$
astart, aend: indices of the first and last rowa of the requested block
bstart, bend: indices of the first and last columns of the requested block
parall: type of parallelization: :serial, ;sharedmem, :distrmem

Returns

The requested block of the posterior covariance.

source