GeostatInversion.jl
This package provides methods for inverse analysis using parameter fields that are represented using geostatistical (stochastic) methods. Currently, two geostatistical methods are implemented. One is the Principal Component Geostatistical Approach (PCGA) proposed by Kitanidis & Lee. The other utilizes a Randomized Geostatistical Approach (RGA) that builds on PCGA.
Randomized Geostatistical Approach (RGA) references:
- O'Malley, D., Le, E., Vesselinov, V.V., Fast Geostatistical Inversion using Randomized Matrix Decompositions and Sketchings for Heterogeneous Aquifer Characterization, AGU Fall Meeting, San Francisco, CA, December 14–18, 2015.
- Lin, Y, Le, E.B, O'Malley, D., Vesselinov, V.V., Bui-Thanh, T., Large-Scale Inverse Model Analyses Employing Fast Randomized Data Reduction, 2016.
Two versions of PCGA are implemented in this package
pcgadirect
, which uses full matrices and direct solvers during iterationspcgalsqr
, which uses low rank representations of the matrices combined with iterative solvers during iterations
The RGA method, rga
, can use either of these approaches using the keyword argument. That is, by doing rga(...; pcgafunc=GeostatInversion.pcgadirect)
or rga(...; pcgafunc=GeostatInversion.pcgalsqr)
.
GeostatInversion.jl module functions:
#
GeostatInversion.getxis
— Function.
Get the parameter subspace that will be explored during the inverse analysis
getxis(samplefield::Function, numfields::Int, numxis::Int, p::Int, q::Int=3, seed=nothing) getxis(Q::Matrix, numxis::Int, p::Int, q::Int=3, seed=nothing)
Arguments:
- samplefield : a function that takes no arguments and returns a sample of the field
- Q : the covariance matrix of the parameter field
- numfields : the number of fields that will be used to find the subspace
- numxis : the dimension of the subspace
- p : oversampling parameter when estimating the range of the covariance matrix (see Halko et al, SIAM Rev., 2011)
- q : number of power iterations when estimating the range of the covariance matrix (see Halko et al, SIAM Rev., 2011)
- seed : an optional seed to use when doing the randomized matrix factorization
#
GeostatInversion.pcgadirect
— Method.
Direct principal component geostatistical approach
pcgadirect(forwardmodel::Function, s0::Vector, X::Vector, xis::Array{Array{Float64, 1}, 1}, R, y::Vector; maxiters::Int=5, delta::Float64=sqrt(eps(Float64)), xtol::Float64=1e-6, callback=(s, obs_cal)->nothing)
Arguments:
- forwardmodel : param to obs map h(s)
- s0 : initial guess
- X : mean of parameter prior (replace with B*X drift matrix later for p>1)
- xis : K columns of Z = randSVDzetas(Q,K,p,q) where Q is the parameter covariance matrix
- R : covariance of measurement error (data misfit term)
- y : data vector
- maxiters : maximum # of PCGA iterations
- delta : the finite difference step size
- xtol : convergence tolerence for the parameters
- callback : a function of the form
(params, observations)->...
that is called during each iteration
#
GeostatInversion.pcgalsqr
— Method.
Iterative principal component geostatistical approach
pcgalsqr(forwardmodel::Function, s0::Vector, X::Vector, xis::Array{Array{Float64, 1}, 1}, R, y::Vector; maxiters::Int=5, delta::Float64=sqrt(eps(Float64)), xtol::Float64=1e-6)
Arguments:
- forwardmodel : param to obs map h(s)
- s0 : initial guess
- X : mean of parameter prior (replace with B*X drift matrix later for p>1)
- xis : K columns of Z = randSVDzetas(Q,K,p,q) where Q is the parameter covariance matrix
- R : covariance of measurement error (data misfit term)
- y : data vector
- maxiters : maximum # of PCGA iterations
- delta : the finite difference step size
- xtol : convergence tolerence for the parameters
#
GeostatInversion.rga
— Method.
Randomized (principal component) geostatistical approach
Example:
function rga(forwardmodel::Function, s0::Vector, X::Vector, xis::Array{Array{Float64, 1}, 1}, R, y::Vector, S; maxiters::Int=5, delta::Float64=sqrt(eps(Float64)), xtol::Float64=1e-6, pcgafunc=pcgadirect, callback=(s, obs_cal)->nothing)
Arguments:
- forwardmodel : param to obs map h(s)
- s0 : initial guess
- X : mean of parameter prior (replace with B*X drift matrix later for p>1)
- xis : K columns of Z = randSVDzetas(Q,K,p,q) where Q is the parameter covariance matrix
- R : covariance of measurement error (data misfit term)
- y : data vector
- S : sketching matrix
- maxiters : maximum # of PCGA iterations
- delta : the finite difference step size
- xtol : convergence tolerance for the parameters
- callback : a function of the form
(params, observations)->...
that is called during each iteration
Module GeostatInversion.FDDerivatives
GeostatInversion.FDDerivatives module functions:
#
GeostatInversion.FDDerivatives.makegradient
— Function.
Create Gradient function
#
GeostatInversion.FDDerivatives.makejacobian
— Function.
Create Jacobian function
Module GeostatInversion.RandMatFact
GeostatInversion.RandMatFact module functions:
#
GeostatInversion.RandMatFact.randsvd
— Method.
Random SVD based on algorithm 5.1 from Halko et al.
Module GeostatInversion.FFTRF
GeostatInversion.FFTRF module functions:
#
GeostatInversion.FFTRF.reducek
— Method.
Reduce k