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}:
 70.98836233316806
 81.8523311692993
 71.78662756848107
 64.01484983571878
 65.17205350312749
 65.50174502937183
 69.4726734799322
 77.93153736094685
 61.520455210387915
 63.78925201847574
  ⋮
 77.55265301781725
 71.4455534962671
 69.20371519966567
 70.64012918320782
 75.58595242655676
 77.59390669788792
 72.26788948777424
 81.3856776881302
 86.55071868869402

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.001 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.0 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}:
 28.554825391325572
 22.62721778567065
 29.44901562928446
 37.78338707663348
 28.03910878407212
 26.83647678045867
 30.59582317530933
 28.20921944302582
 26.926585560445254
 32.5931505392938
  ⋮
 28.224344258363832
 23.025048650962134
 31.45243725447051
 39.132902316768835
 28.104972156245303
 28.905607999295327
 31.97508837086241
 28.20336426741758
 26.68559885299182

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