| Title: | Geophysical Inverse Theory and Optimization |
|---|---|
| Description: | Several functions introduced in Aster et al.'s book on inverse theory. The functions are often translations of MATLAB code developed by the authors to illustrate concepts of inverse theory as applied to geophysics. Generalized inversion, tomographic inversion algorithms (conjugate gradients, 'ART' and 'SIRT'), non-linear least squares, first and second order Tikhonov regularization, roughness constraints, and procedures for estimating smoothing parameters are included. |
| Authors: | Jonathan M. Lees [aut, cre] |
| Maintainer: | Jonathan M. Lees <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.2-5 |
| Built: | 2024-11-13 03:57:51 UTC |
| Source: | https://github.com/cran/PEIP |
Auxilliary functions and routines for running the examples and excersizes described in the book on inverse theory.
These functions are used in conjunction with the example described in the PEIP book.
There is one C-code routine, interp2grid. This is introduced to replicate the MATLAB code interp2. It does not work exactly as the matlab code prescribes.
In the PEIP library one LAPACK routine is called: dggsvd. In R, LAPACK routines are stored in slightly different locations on Linux, Windows and Mac computers. Be aware. This will come up in examples from Chapter 4.
Almost all examples work as scripts run with virtually no user input, e.g.
Jonathan M. Lees<jonathan.lees.edu> Maintainer:Jonathan M. Lees<jonathan.lees.edu>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
QR decomposition solution to Ax=b
Ainv(GAB, x, tol = 1e-12)Ainv(GAB, x, tol = 1e-12)
GAB |
design matrix |
x |
right hand side |
tol |
tolerance for singularity |
I needed something to make up for the lame-o matlab code that does this h = G\x to get the inverse
Inverse Solution
Jonathan M. Lees<[email protected]>
set.seed(2015) GAB = matrix(runif(36), ncol=6) truex =rnorm(ncol(GAB)) rhs = GAB %*% truex rhs = as.vector(rhs ) tout = Ainv(GAB, rhs, tol = 1e-12) tout - truexset.seed(2015) GAB = matrix(runif(36), ncol=6) truex =rnorm(ncol(GAB)) rhs = GAB %*% truex rhs = as.vector(rhs ) tout = Ainv(GAB, rhs, tol = 1e-12) tout - truex
ART algorythm for solving sparse linear inverse problems
art(A, b, tolx, maxiter)art(A, b, tolx, maxiter)
A |
Constraint matrix |
b |
right hand side |
tolx |
difference tolerance for successive iterations (stopping criteria) |
maxiter |
maximum iterations (stopping criteria). |
Alpha is a damping factor. If alpha<1, then we won't take full steps in the ART direction. Using a smaller value of alpha (say alpha=.75) can help with convergence on some problems.
x |
solution |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
set.seed(2015) G = setDesignG() ### % Setup the true model. mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### % reshape the true model to be a vector mtruev=as.vector(mtruem); ### % Compute the data. dtrue=G %*% mtruev; ### % Add the noise. d=dtrue+0.01*rnorm(length(dtrue)); mkac<-art(G,d,0.01,200) par(mfrow=c(1,2)) imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); imagesc(matrix(mkac,16,16) , asp=1 , main="ART Solution" );set.seed(2015) G = setDesignG() ### % Setup the true model. mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### % reshape the true model to be a vector mtruev=as.vector(mtruem); ### % Compute the data. dtrue=G %*% mtruev; ### % Add the noise. d=dtrue+0.01*rnorm(length(dtrue)); mkac<-art(G,d,0.01,200) par(mfrow=c(1,2)) imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); imagesc(matrix(mkac,16,16) , asp=1 , main="ART Solution" );
Bartlett (triangle) window of length m
bartl(m)bartl(m)
m |
integer, length of vector |
vector
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
bartl(11)bartl(11)
Given a linear inverse problem Gm=d, a prior mean mprior and covariance matrix covm, data d, and data covariance matrix covd, this function computes the MAP solution and the corresponding covariance matrix.
bayes(G, mprior, covm, d, covd)bayes(G, mprior, covm, d, covd)
G |
Design Matrix |
mprior |
vector, prior model |
covm |
vector, model covariance |
d |
vector, right hand side |
covd |
vector, data covariance |
vector model
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
## Not run: set.seed(2015) G = setDesignG() ### mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### mtruev=as.vector(mtruem); imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); matrix(mtruem,16,16) , asp=1 , main="True Model" ) ### dtrue=G %*% mtruev; ### d=dtrue+0.01*rnorm(length(dtrue)); covd = 0.1*diag( nrow=length(d) ) covm = 1*diag( nrow=dim(G)[2] ) ## End(Not run)## Not run: set.seed(2015) G = setDesignG() ### mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### mtruev=as.vector(mtruem); imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); matrix(mtruem,16,16) , asp=1 , main="True Model" ) ### dtrue=G %*% mtruev; ### d=dtrue+0.01*rnorm(length(dtrue)); covd = 0.1*diag( nrow=length(d) ) covm = 1*diag( nrow=dim(G)[2] ) ## End(Not run)
Bounded least squares
blf2(A, b, c, delta, l, u)blf2(A, b, c, delta, l, u)
A |
Design Matrix |
b |
Right hand side |
c |
matrix weight on x |
delta |
tolerance |
l |
lower bound |
u |
upper bound |
Solves the problem: min/max c'*x where || Ax-b || <= delta and l <= x <= u.
x |
solution |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
Stark, P.B. , and R. L. Parker, Bounded-Variable Least-Squares: An Algorithm and Applications, Computational Statistics 10:129-141, 1995.
### set up an inverse problem:Shaw problem n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 spiken = G %*% spike wts = rep(1, n) delta = 1e-03 set.seed(2015) dspiken = spiken + 6e-6 *rnorm(length(spiken)) lb = spike - (.2) * wts ub = spike + (.2) * wts dspiken = dspiken blf2(G, dspiken, wts , delta, lb, ub)### set up an inverse problem:Shaw problem n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 spiken = G %*% spike wts = rep(1, n) delta = 1e-03 set.seed(2015) dspiken = spiken + 6e-6 *rnorm(length(spiken)) lb = spike - (.2) * wts ub = spike + (.2) * wts dspiken = dspiken blf2(G, dspiken, wts , delta, lb, ub)
Conjugate gradient Least squares
cgls(Gmat, dee, niter)cgls(Gmat, dee, niter)
Gmat |
input matrix |
dee |
right hand side |
niter |
max number of iterations |
Performs niter iterations of the CGLS algorithm on the least squares problem min norm(G*m-d). Gmat should be a sparse matrix.
X |
matrix of models |
rho |
misfit norms |
eta |
model norms |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
set.seed(11) #### perfect data with no noise n <- 5 A <- matrix(runif(n*n),nrow=n) B <- runif(n) ### get right-hand-side (data) trhs = as.vector( A %*% B ) Lout = cgls(A, trhs , 15) ### solution is Lout$X[,15] Lout$X[,15] - Bset.seed(11) #### perfect data with no noise n <- 5 A <- matrix(runif(n*n),nrow=n) B <- runif(n) ### get right-hand-side (data) trhs = as.vector( A %*% B ) Lout = cgls(A, trhs , 15) ### solution is Lout$X[,15] Lout$X[,15] - B
Chi function
chi(x, n)chi(x, n)
x |
value |
n |
degrees of freedom |
function evaluated
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
x = seq(0, 10, length=100) n = 5 y=chi(x, n) plot(x, y)x = seq(0, 10, length=100) n = 5 y=chi(x, n) plot(x, y)
Computes the Chi^2 CDF, using a transformation to N(0,1) on page 333 of Thistead, Elements of Statistical Computing.
chi2cdf(x, n)chi2cdf(x, n)
x |
end value of chi^2 pdf to integrate to. (scalar) |
n |
degrees of freedom (scalar) |
Note that x and m must be scalars.
p |
probability that Chi^2 random variable is less than or equal to x (scalar). |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
x= seq(from=0.1, to=0.9, length=20) chi2cdf(x , 3)x= seq(from=0.1, to=0.9, length=20) chi2cdf(x , 3)
Inverse Chi-Sq
chi2inv(x, n)chi2inv(x, n)
x |
probability that Chi^2 random variable is less than or equal to x (scalar). |
n |
degrees of freedom(scalar) |
Computes the inverse Chi^2 distribution corresponding to a given probability that a Chi^2 random variable with the given degrees of freedom is less than or equal to x. Uses chi2cdf.m.
corresponding value of x for given probability.
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
chi, chi2cdf
x = seq(from=0.1, to=0.9, length=10) h = chi2cdf(x, 3) chi2inv(h, 3)x = seq(from=0.1, to=0.9, length=10) h = chi2cdf(x, 3) chi2inv(h, 3)
Computes the column-by-column discrete cosine transform of X.
dcost(X)dcost(X)
X |
Time series matrix |
cosine transformaed data
Jonathan M. Lees<[email protected]>
x <- 1:4 ### compare fft with cosine transform fft(x) dcost(x)x <- 1:4 ### compare fft with cosine transform fft(x) dcost(x)
Plot Error Bar
error.bar(x, y, lo, hi, pch = 1, col = 1, barw = 0.1, add = FALSE, ...)error.bar(x, y, lo, hi, pch = 1, col = 1, barw = 0.1, add = FALSE, ...)
x |
X-values |
y |
Y-values |
lo |
Lower limit of error bars |
hi |
Upper limit of error bars |
pch |
plotting character |
col |
color |
barw |
width of the bar |
add |
logical, add=FALSE starts a new plot |
... |
other plotting parameters |
graphical side effects
Jonathan M. Lees<[email protected]>
x = 1:10 y = 2*x+5 zup = rnorm(10) zup = zup-min(zup)+.5 zdown = rnorm(10) zdown = zdown-min(zdown)+.2 #### example with same error on either side: error.bar(x, y, y-zup, y+zup, pch = 1, col = 'brown' , barw = 0.1, add = FALSE) #### example with different error on either side: error.bar(x, y, y-zdown, y+zup, pch = 1, col = 'brown' , barw = 0.1, add = FALSE)x = 1:10 y = 2*x+5 zup = rnorm(10) zup = zup-min(zup)+.5 zdown = rnorm(10) zdown = zdown-min(zdown)+.2 #### example with same error on either side: error.bar(x, y, y-zup, y+zup, pch = 1, col = 'brown' , barw = 0.1, add = FALSE) #### example with different error on either side: error.bar(x, y, y-zdown, y+zup, pch = 1, col = 'brown' , barw = 0.1, add = FALSE)
Flip (reverse order) output of GSVD
flipGSVD(vs, d1 = c(50, 50), d2 = c(48, 50))flipGSVD(vs, d1 = c(50, 50), d2 = c(48, 50))
vs |
list output of GSVD |
d1 |
dimensionals of A |
d2 |
dimensions of B |
This flipping of the matrix is done to agree with the Matlab code.
Same as GSVD, but order of eigenvectors is reversed.
U |
m by m orthogonal matrix |
V |
p by p orthogonal matrix, p=rank(B) |
X |
n by n nonsingular matrix |
C |
singular values, m by n matrix with diagonal elements shifted from main diagonal |
S |
singular values, p by n diagonal matrix |
The GSVD routines are from LAPACK.
Jonathan M. Lees<[email protected]>
GSVD
set.seed(12) n <- 5 A <- matrix(runif(n*n),nrow=n) B <- matrix(runif(n*n),nrow=n) VS = GSVD(A, B) FVS = flipGSVD(VS, d1 = dim(A) , d2 = dim(B) ) ## see that order of eigen vectors is reversed diag(VS$S) diag(FVS$S)set.seed(12) n <- 5 A <- matrix(runif(n*n),nrow=n) B <- matrix(runif(n*n),nrow=n) VS = GSVD(A, B) FVS = flipGSVD(VS, d1 = dim(A) , d2 = dim(B) ) ## see that order of eigen vectors is reversed diag(VS$S) diag(FVS$S)
Auxiliary routine for GCV calculations
gcv_function(alpha, gamma2, beta)gcv_function(alpha, gamma2, beta)
alpha |
parameter |
gamma2 |
square of the gamma from the gsvd |
beta |
projected data to fit |
vector, g - || Gm_(alpha,L) - d ||^2 / (Tr(I - GG#)^2
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
Extract the smallest regularization parameter.
gcval(U, s, b, npoints)gcval(U, s, b, npoints)
U |
U matrix from gsvd(G, L) |
s |
[diag(C) diag(S)] which are the lambdas and mus from the gsvd |
b |
the data to try and match |
npoints |
number of alphas to estimate |
Evaluate the GCV function gcv_function at npoints points.
List:
reg_min |
alpha with the minimal g (scalar) |
g |
|| Gm_(alpha,L) - d ||^2 / (Tr(I - GG#)^2 |
alpha |
alpha for the corresponding g |
Jonathan M. Lees<[email protected]>
gcv_function
set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S lam=sqrt(diag(t(Lam1 %*% Lam1))); mu=sqrt(diag(t(M1)%*%M1)); p=rnk(L1); sm1=cbind(lam[1:p],mu[1:p]) ### % get the gcv values varying alpha ### ngcvpoints=1000; HI = gcval(U1,sm1,t,ngcvpoints);set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S lam=sqrt(diag(t(Lam1 %*% Lam1))); mu=sqrt(diag(t(M1)%*%M1)); p=rnk(L1); sm1=cbind(lam[1:p],mu[1:p]) ### % get the gcv values varying alpha ### ngcvpoints=1000; HI = gcval(U1,sm1,t,ngcvpoints);
Returns a 1D differentiating matrix operating on a series with n points.
get_l_rough(n, deg)get_l_rough(n, deg)
n |
integer, number of data points |
deg |
order of the derivative to approximate |
Used to get first and 2nd order roughening matrices for 1-D problems
Matrix:discrete differentiation matrix
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
### first order roughening matrix for a 10 by 10 model: a sparse matrix N = 10 L1 = get_l_rough(10,1); ### second order roughening matrix for a 10 by 10 model N = 10 L2 = get_l_rough(10,2);### first order roughening matrix for a 10 by 10 model: a sparse matrix N = 10 L1 = get_l_rough(10,1); ### second order roughening matrix for a 10 by 10 model N = 10 L2 = get_l_rough(10,2);
Get inverse of matrisx or solve Ax=b.
ginv(G, x, tol = 1e-12)ginv(G, x, tol = 1e-12)
G |
Design Matrix |
x |
right hand side |
tol |
tolerance |
This function used as alternative to matlab code that does this h = G\x to get the inverse
inverse as a N by 1 matrix.
Be careful about the usage of tolerance
Jonathan M. Lees<[email protected]>
solve, Ainv
set.seed(2015) GAB = matrix(runif(36), ncol=6) truex =rnorm(ncol(GAB)) rhs = GAB %*% truex rhs = as.vector(rhs ) tout = ginv(GAB, rhs, tol = 1e-12) tout - truexset.seed(2015) GAB = matrix(runif(36), ncol=6) truex =rnorm(ncol(GAB)) rhs = GAB %*% truex rhs = as.vector(rhs ) tout = ginv(GAB, rhs, tol = 1e-12) tout - truex
Wrapper for generalized svd from LAPACK
GSVD(A, B)GSVD(A, B)
A |
Matrix, see below |
B |
Matrix, see below |
The A and B matrices will be, A=U*C*t(X) and B=V*S*t(X), respectively.
Since PEIP is based on a book, which is iteslef based on MATLAB routines, the convention here follows the book. The R implementation uses LAPACK and wraps the function so the output will comply with the book. See page 104 of the second edition of the Aster book cited below. That said, the purpose is to find an inversion of the form Y = t(A aB), where a is a regularization parameter, B is smoothing matrix and A is the design matrix for the forward problem. The input matrices A and B are assumed to have full rank, and p = rank(B). The generalized singular values are then gamma = lambda/mu, where lambda = sqrt(diag(t(C)*C) ) and mu = sqrt(diag(t(S)*S) ).
U |
m by m orthogonal matrix |
V |
p by p orthogonal matrix, p=rank(B) |
X |
n by n nonsingular matrix |
C |
singular values, m by n matrix with diagonal elements shifted from main diagonal |
S |
singular values, p by n diagonal matrix |
Requires R version of LAPACK. The code is a wrapper for the dggsvd function in LAPACK. The author thanks Berend Hasselman for advice and help preparing this function.
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
flipGSVD
# Example from NAG F08VAF A <- matrix(1:15, nrow=5,ncol=3) B <- matrix(c(8,1,6, 3,5,7, 4,9,2), nrow=3, byrow=TRUE) z <- GSVD(A,B) C <- z$C S <- z$S sqrt(diag(t(C) %*% C)) / sqrt(diag(t(S) %*% S)) testA = A - z$U %*% C %*% t(z$X) testB = B - z$V %*% S %*% t(z$X) print(testA) print(testB)# Example from NAG F08VAF A <- matrix(1:15, nrow=5,ncol=3) B <- matrix(c(8,1,6, 3,5,7, 4,9,2), nrow=3, byrow=TRUE) z <- GSVD(A,B) C <- z$C S <- z$S sqrt(diag(t(C) %*% C)) / sqrt(diag(t(S) %*% S)) testA = A - z$U %*% C %*% t(z$X) testB = B - z$V %*% S %*% t(z$X) print(testA) print(testB)
Takes the column-by-column inverse discrete cosine transform of Y.
idcost(Y)idcost(Y)
Y |
Input cosine transform |
Time series
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
dcost
x <- 1:4 ### compare fft with cosine transform fft(x) zig = dcost(x) zag = idcost(zig)x <- 1:4 ### compare fft with cosine transform fft(x) zig = dcost(x) zag = idcost(zig)
Display image in matlab format, i.e. flip and transpose.
imagesc(G, col = grey((1:99)/100), ...) contoursc(G, ...)imagesc(G, col = grey((1:99)/100), ...) contoursc(G, ...)
G |
Image matrix |
col |
color scale |
... |
graphical parameters |
Program flips image and transposes prior to plotting. The contour version does the same and can be used to add contours.
graphical side effects
Jonathan M. Lees<[email protected]>
mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; imagesc(mtruem, asp=1)mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; imagesc(mtruem, asp=1)
This code includes a bicubic interpolation and a bilinear interpolation adapted from Numerical Recipes in C: The art of scientific computing (chapter 3... bicubic interpolation) and a bicubic interpolation from in java code.
Inputs are a list of points to interpolate to and from raster objects of class 'asc' (adehabitat package), 'RasterLayer' (raster package) or 'SpatialGridDataFrame' (sp package).
interp2grid(mat,xout,yout,xin=NULL,yin=NULL,type=2)interp2grid(mat,xout,yout,xin=NULL,yin=NULL,type=2)
mat |
a matrix of data that can be a raster matrix of class 'asc' (adehabitat package), 'RasterLayer' (raster package) or 'SpatialGridDataFrame' (sp package) NA values are not permitted.. data must be complete. |
xout |
a vector of data representing x coordinates of the output grid. Resulting grid must have square cell sizes if mat is of class 'asc', 'RasterLayer' or 'SpatialGridDataFrame'. |
yout |
a vector of data representing x coordinates of the output grid. Resulting grid must have square cell sizes if mat is of class 'asc', 'RasterLayer' or 'SpatialGridDataFrame'. |
xin |
a vector identifying the locations of the columns of the input data matrix. These are automatically populated if mat is of class 'asc', 'RasterLayer' or 'SpatialGridDataFrame'. |
yin |
a vector identifying the locations of the rows of the input data matrix. These are automatically populated if mat is of class 'asc', 'RasterLayer' or 'SpatialGridDataFrame'. |
type |
an integer value representing the type of interpolation method used. 1 - bilinear adapted from Numerical Recipes in C 2 - bicubic adapted from Numerical Recipes in C 3 - bicubic adapted from online java code |
Returns a matrix of the originating class.
Jeremy VanDerWal [email protected]
tx = seq(0,3,0.1) ty = seq(0,3,0.1) tmat = matrix(runif(16,1,16),nrow=4) txin = seq(0,3,length=4) tyin = seq(0,3,length=4) bilinear1 = interp2grid(tmat,tx,ty,txin, tyin, type=1) bicubic2 = interp2grid(tmat,tx,ty,txin, tyin, type=2) bicubic3 = interp2grid(tmat,tx,ty,txin, tyin, type=3) par(mfrow=c(2,2),cex=1) image(tmat,main='base',zlim=c(0,16),col=heat.colors(100)) image(bilinear1,main='bilinear',zlim=c(0,16),col=heat.colors(100)) image(bicubic2,main='bicubic2',zlim=c(0,16),col=heat.colors(100)) image(bicubic3,main='bicubic3',zlim=c(0,16),col=heat.colors(100))tx = seq(0,3,0.1) ty = seq(0,3,0.1) tmat = matrix(runif(16,1,16),nrow=4) txin = seq(0,3,length=4) tyin = seq(0,3,length=4) bilinear1 = interp2grid(tmat,tx,ty,txin, tyin, type=1) bicubic2 = interp2grid(tmat,tx,ty,txin, tyin, type=2) bicubic3 = interp2grid(tmat,tx,ty,txin, tyin, type=3) par(mfrow=c(2,2),cex=1) image(tmat,main='base',zlim=c(0,16),col=heat.colors(100)) image(bilinear1,main='bilinear',zlim=c(0,16),col=heat.colors(100)) image(bicubic2,main='bicubic2',zlim=c(0,16),col=heat.colors(100)) image(bicubic3,main='bicubic3',zlim=c(0,16),col=heat.colors(100))
Uses the iteratively reweight least squares strategy to find an approximate L_p solution to Ax=b.
irls(A, b, tolr, tolx, p, maxiter)irls(A, b, tolr, tolx, p, maxiter)
A |
Matrix of the system of equations. |
b |
Right hand side of the system of equations |
tolr |
Tolerance below which residuals are ignored |
tolx |
Stopping tolerance. Stop when (norm(newx-x)/(1+norm(x)) < tolx) |
p |
Specifies which p-norm to use (most often, p=1.) |
maxiter |
Limit on number of iterations of IRLS |
Use to get L-1 norm solution of inverse problems.
x |
Approximate L_p solution |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
t = 1:10 y=c(109.3827,187.5385,267.5319,331.8753,386.0535, 428.4271,452.1644,498.1461,512.3499,512.9753) sigma = rep(8, length(y)) N=length(t); ### % Introduce the outlier y[4]=y[4]-200; G = cbind( rep(1, N), t, -1/2*t^2 ) ### % Apply the weighting yw = y/sigma; Gw = G/sigma m2 = solve( t(Gw) %*% Gw , t(Gw) %*% yw, tol=1e-12 ) ### Solve for the 1-norm solution m1 = irls(Gw,yw,1.0e-5,1.0e-5,1,25) m1t = 1:10 y=c(109.3827,187.5385,267.5319,331.8753,386.0535, 428.4271,452.1644,498.1461,512.3499,512.9753) sigma = rep(8, length(y)) N=length(t); ### % Introduce the outlier y[4]=y[4]-200; G = cbind( rep(1, N), t, -1/2*t^2 ) ### % Apply the weighting yw = y/sigma; Gw = G/sigma m2 = solve( t(Gw) %*% Gw , t(Gw) %*% yw, tol=1e-12 ) ### Solve for the 1-norm solution m1 = irls(Gw,yw,1.0e-5,1.0e-5,1,25) m1
Solves the system Gm=d using sparsity regularization on Lm. Solves the L1 regularized least squares problem: min norm(G*m-d,2)^2+alpha*norm(L*m,1)
irlsl1reg(G, d, L, alpha, maxiter = 100, tolx = 1e-04, tolr = 1e-06)irlsl1reg(G, d, L, alpha, maxiter = 100, tolx = 1e-04, tolr = 1e-06)
G |
design matrix |
d |
right hand side |
L |
regularization matrix |
alpha |
regularization parameter |
maxiter |
Maximum number of IRLS iterations |
tolx |
Tolerance on successive iterates |
tolr |
Tolerance below which we consider an element of L*m to be effectively zero |
m |
model vector |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 spiken = G %*% spike wts = rep(1, n) delta = 1e-03 set.seed(2015) dspiken = spiken + 6e-6 *rnorm(length(spiken)) L1 = get_l_rough(n,1); alpha = 0.001 k = irlsl1reg(G, dspiken, L1, alpha, maxiter = 100, tolx = 1e-04, tolr = 1e-06) plotconst(k,-pi/2,pi/2, ylim=c(-.2, 0.5), xlab="theta", ylab="Intensity" );n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 spiken = G %*% spike wts = rep(1, n) delta = 1e-03 set.seed(2015) dspiken = spiken + 6e-6 *rnorm(length(spiken)) L1 = get_l_rough(n,1); alpha = 0.001 k = irlsl1reg(G, dspiken, L1, alpha, maxiter = 100, tolx = 1e-04, tolr = 1e-06) plotconst(k,-pi/2,pi/2, ylim=c(-.2, 0.5), xlab="theta", ylab="Intensity" );
Implements Kaczmarz's algorithm to solve a system of equations iteratively
kac(A, b, tolx, maxiter)kac(A, b, tolx, maxiter)
A |
Constraint matrix |
b |
right hand side |
tolx |
difference tolerence for successive iterations (stopping criteria) |
maxiter |
maximum iterations (stopping criteria) |
x |
solution |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
set.seed(2015) G = setDesignG() ### % Setup the true model. mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### % reshape the true model to be a vector mtruev=as.vector(mtruem); ### % Compute the data. dtrue=G %*% mtruev; ### % Add the noise. d=dtrue+0.1*rnorm(length(dtrue)); mkac<-kac(G,d,0.0,200) par(mfrow=c(1,2)) imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); imagesc(matrix(mkac,16,16) , asp=1 , main="Kacz Solution" );set.seed(2015) G = setDesignG() ### % Setup the true model. mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### % reshape the true model to be a vector mtruev=as.vector(mtruem); ### % Compute the data. dtrue=G %*% mtruev; ### % Add the noise. d=dtrue+0.1*rnorm(length(dtrue)); mkac<-kac(G,d,0.0,200) par(mfrow=c(1,2)) imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); imagesc(matrix(mkac,16,16) , asp=1 , main="Kacz Solution" );
Retrieve corner of L-curve
l_curve_corner(rho, eta, reg_param)l_curve_corner(rho, eta, reg_param)
rho |
misfit |
eta |
model norm or seminorm |
reg_param |
regularization parameter |
reg_corner |
the value of reg_param with maximum curvature |
ireg_corner |
the index of the value in reg_param with maximum curvature |
kappa |
the curvature for each reg_param |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tgsvd(U1,t,X1,Lam1,G,L1); rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### % store where the corner is (from visual inspection) vcorn = l_curve_corner(rho1, eta1, reg_param1) ireg_corner1=vcorn$reg_corner rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste('1st order reg corner is: ',ireg_corner1)); plot(rho1,eta1,type="b", log="xy" , xlim=c(1e-4, 1e-2) , ylim=c(6e-6, 2e-4) , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" ); points(rho_corner1, eta_corner1, col='red', cex=2 )#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tgsvd(U1,t,X1,Lam1,G,L1); rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### % store where the corner is (from visual inspection) vcorn = l_curve_corner(rho1, eta1, reg_param1) ireg_corner1=vcorn$reg_corner rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste('1st order reg corner is: ',ireg_corner1)); plot(rho1,eta1,type="b", log="xy" , xlim=c(1e-4, 1e-2) , ylim=c(6e-6, 2e-4) , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" ); points(rho_corner1, eta_corner1, col='red', cex=2 )
L curve parematers and models for truncated gsvd regularization.
l_curve_tgsvd(U, d, X, Lam, G, L)l_curve_tgsvd(U, d, X, Lam, G, L)
U |
U, output of GSVD |
d |
output of GSVD |
X |
output of GSVD |
Lam |
output of GSVD |
G |
output of GSVD |
L |
output of GSVD |
List:
eta |
the solution seminorm ||Lm|| |
rho |
the residual norm ||G m - d|| |
reg_param |
corresponding regularization parameters |
m |
corresponding suite of models for truncated GSVD |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tgsvd(U1,t,X1,Lam1,G,L1); rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### % store where the corner is (from visual inspection) ireg_corner1=8; rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste('1st order reg corner is: ',ireg_corner1)); plot(rho1,eta1,type="b", log="xy", xlim=c(1e-4, 1e-2) , ylim=c(6e-6, 2e-4) , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" );#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tgsvd(U1,t,X1,Lam1,G,L1); rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### % store where the corner is (from visual inspection) ireg_corner1=8; rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste('1st order reg corner is: ',ireg_corner1)); plot(rho1,eta1,type="b", log="xy", xlim=c(1e-4, 1e-2) , ylim=c(6e-6, 2e-4) , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" );
L-curve tikh gsvd
l_curve_tikh_gsvd(U, d, X, Lam, Mu, G, L, npoints, varargin = NULL)l_curve_tikh_gsvd(U, d, X, Lam, Mu, G, L, npoints, varargin = NULL)
U |
from the gsvd |
d |
data vector for the problem G*m=d |
X |
from the gsvd |
Lam |
from the gsvd |
Mu |
from the gsvd |
G |
system matrix |
L |
roughening matrix |
npoints |
Number of points |
varargin |
alpha_min, alpha_max: if specified, constrain the logrithmically spaced regularization parameter range, otherwise an attempt is made to estimate them from the range of generalized singular values |
Uses output of GSVD
eta |
- the solution seminorm ||Lm|| |
rho |
- the residual norm ||G m - d|| |
reg_param |
- corresponding regularization parameters |
m |
- corresponding suite of models for truncated GSVD |
Jonathan M. Lees<[email protected]>
#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tikh_gsvd(U1,t,X1,Lam1,M1, G,L1, 25); rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### store where the corner is (from visual inspection) ireg_corner1=8; rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste('1st order reg corner is: ',ireg_corner1)); plot(rho1,eta1,type="b", log="xy", xlim=c(1e-4, 1e-2) , ylim=c(6e-6, 2e-4) , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" );#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tikh_gsvd(U1,t,X1,Lam1,M1, G,L1, 25); rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### store where the corner is (from visual inspection) ireg_corner1=8; rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste('1st order reg corner is: ',ireg_corner1)); plot(rho1,eta1,type="b", log="xy", xlim=c(1e-4, 1e-2) , ylim=c(6e-6, 2e-4) , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" );
L-curve for Tikhonov regularization
l_curve_tikh_svd(U, s, d, npoints, varargin = NULL)l_curve_tikh_svd(U, s, d, npoints, varargin = NULL)
U |
matrix of data space basis vectors from the svd |
s |
vector of singular values |
d |
the data vector |
npoints |
the number of logarithmically spaced regularization parameters |
varargin |
alpha_min, alpha_max: if specified, constrain the logrithmically spaced regularization parameter range, otherwise an attempt is made to estimate them from the range of singular values |
Calculates the L-curve
eta |
the solution norm ||m|| or seminorm ||Lm|| |
rho |
the residual norm ||G m - d|| |
reg_param |
corresponding regularization parameters |
Jonathan M. Lees<[email protected]>
#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tikh_svd(U1, diag(M1) , X1, 25, varargin = NULL) rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### store where the corner is (from visual inspection) ireg_corner1=8; rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste("1st order reg corner is: ",ireg_corner1)); plot(rho1,eta1,type="b", log="xy" , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" );#### Vertical Seismic Profile example set.seed(2015) VSP = vspprofile() t = VSP$t2 G = VSP$G M = VSP$M N = VSP$N L1 = get_l_rough(N,1); littleU = PEIP::GSVD(as.matrix(G), as.matrix(L1) ); BIGU = flipGSVD(littleU, dim(G), dim(L1) ) U1 = BIGU$U V1 =BIGU$V X1=BIGU$X Lam1=BIGU$C M1=BIGU$S K1 = l_curve_tikh_svd(U1, diag(M1) , X1, 25, varargin = NULL) rho1 =K1$rho eta1 =K1$eta reg_param1 =K1$reg_param m1s =K1$m ### store where the corner is (from visual inspection) ireg_corner1=8; rho_corner1=rho1[ireg_corner1]; eta_corner1=eta1[ireg_corner1]; print(paste("1st order reg corner is: ",ireg_corner1)); plot(rho1,eta1,type="b", log="xy" , xlab="Residual Norm ||Gm-d||_2", ylab="Solution Seminorm ||Lm||_2" );
Add to plotting model in piecewise constant form over n subintervals, where n is the length of x.
linesconst(x, l, r, ...)linesconst(x, l, r, ...)
x |
model to be plotted |
l |
left endpoint of plot |
r |
right endpoint of plot |
... |
graphical parameters |
Used for plotting vector models
graphical side effects
Jonathan M. Lees<[email protected]>
plotconst
zip = runif(25) plotconst(zip, 0, 1 ) linesconst(runif(25) , 0, 1 , col='red' )zip = runif(25) plotconst(zip, 0, 1 ) linesconst(runif(25) , 0, 1 , col='red' )
Use the Levenberg-Marquardt algorithm to minimize f(p)=sum(F_i(p)^2)
lmarq(afun, ajac, p0, tol, maxiter)lmarq(afun, ajac, p0, tol, maxiter)
afun |
name of the function F(x) |
ajac |
name of the Jacobian function J(x) |
p0 |
initial guess |
tol |
stopping tolerance |
maxiter |
maximum number of iterations allowed |
pstar |
best solution found. |
iter |
Iteration count. |
Jonathan M. Lees<[email protected]>
fun<-function(p){ ### Compute the function values. fvec=rep(0,length(TM)) fvec=(Q*exp(-D^2*p[1]/(4*p[2]*TM))/(4*pi*p[2]*TM) - H)/SIGMA return(fvec) } jac <-function( p) { ### use known formula for the derivatives in the Jacobian n=length(TM) J= matrix(0,nrow=n,ncol=2) J[,1]=(-Q*D^2*exp(-D^2*p[1]/(4*p[2]*TM))/(16*pi*p[2]^2*TM^2))/SIGMA J[,2]=(Q/(4*pi*p[2]^2*TM))* ((D^2*p[1])/(4*p[2]*TM)-1)*exp(-D^2*p[1]/(4*p[2]*TM))/SIGMA return(J) } H=c(0.72, 0.49, 0.30, 0.20, 0.16, 0.12) TM=c(5.0, 10.0, 20.0, 30.0, 40.0, 50.0) ### Fixed parameter values. D=60 Q=50 ### We'll use sigma=1cm. SIGMA=0.01*rep(1,length(H)) ### The unknown/estimated parameters are S=p(1) and T=p(2). p0=c(0.001, 1.0) ### Solve the least squares problem with LM. PEST = lmarq('fun','jac',p0,1.0e-12,100)fun<-function(p){ ### Compute the function values. fvec=rep(0,length(TM)) fvec=(Q*exp(-D^2*p[1]/(4*p[2]*TM))/(4*pi*p[2]*TM) - H)/SIGMA return(fvec) } jac <-function( p) { ### use known formula for the derivatives in the Jacobian n=length(TM) J= matrix(0,nrow=n,ncol=2) J[,1]=(-Q*D^2*exp(-D^2*p[1]/(4*p[2]*TM))/(16*pi*p[2]^2*TM^2))/SIGMA J[,2]=(Q/(4*pi*p[2]^2*TM))* ((D^2*p[1])/(4*p[2]*TM)-1)*exp(-D^2*p[1]/(4*p[2]*TM))/SIGMA return(J) } H=c(0.72, 0.49, 0.30, 0.20, 0.16, 0.12) TM=c(5.0, 10.0, 20.0, 30.0, 40.0, 50.0) ### Fixed parameter values. D=60 Q=50 ### We'll use sigma=1cm. SIGMA=0.01*rep(1,length(H)) ### The unknown/estimated parameters are S=p(1) and T=p(2). p0=c(0.001, 1.0) ### Solve the least squares problem with LM. PEST = lmarq('fun','jac',p0,1.0e-12,100)
Load a Matlab matfile, rename the internal parameters to get R-objects
loadMAT(fn, pos=1)loadMAT(fn, pos=1)
fn |
file name of MATfile |
pos |
integer, position in search path, default=1 |
Program reads in previously saved mat-files and extracts the data, and renames the variables to match the book.
Whatever is in the MATfile
Matfiles are created using the matlab2R routines
Jonathan M. Lees<[email protected]>
Maximum likelihood Models
mcmc(alogprior, aloglikelihood, agenerate, alogproposal, m0, niter)mcmc(alogprior, aloglikelihood, agenerate, alogproposal, m0, niter)
alogprior |
Name of a function that computes the log of the prior distribution. |
aloglikelihood |
Name of a function the computes the log of the likelihood. |
agenerate |
Name of a function that generates a random model from the current model using the |
alogproposal |
Name of a function that computes the log of the proposal distribution r(x,y). |
m0 |
Initial model |
niter |
Number of iterations to perform |
mout |
MCMC samples |
mMAP |
Best model found in the MCMC simulation. |
accrate |
Acceptance rate |
Jonathan M. Lees<[email protected]>
fun <-function(m,x) { y=m[1]*exp(m[2]*x)+m[3]*x*exp(m[4]*x) return(y) } generate <-function( x) { y=x+step*rnorm(4) return(y) } logprior <-function(m) { if( (m[1]>=0) & (m[1]<=2) & (m[2]>=-0.9) & (m[2]<=0) & (m[3]>=0) & (m[3]<=2) & (m[4]>=-0.9) & (m[4]<=0) ) { lp=0 } else { lp= -Inf } return(lp) } loglikelihood <-function(m) { fvec=(y-fun(m,x))/sigma L=(-1/2)*sum(fvec^2) return(L) } logproposal <-function(x,y) { LR=(-1/2)*sum((x-y)^2/step^2) return(LR) } ### Generate the data set. x=seq(from=1, by=0.25, to=7.0) mtrue=c(1.0, -0.5, 1.0, -0.75) ytrue=fun(mtrue,x) sigma=0.01*rep(1, times= length(ytrue) ) y=ytrue+sigma*rnorm(length(ytrue) ) ### set the MCMC parameters ### number of skips to reduce autocorrelation of models skip=100 ### burn-in steps BURNIN=1000 ### number of posterior distribution samples N=4100 ### MVN step size step = 0.005*rep(1,times=4) ### We assume flat priors here m0 = c(0.9003, -0.5377, 0.2137, -0.0280) alogprior='logprior' aloglikelihood='loglikelihood' agenerate='generate' alogproposal='logproposal' ### ### initialize model at a random point on [-1,1] ### m0=(runif(4)-0.5)*2 ### this is the matlab initialization: m0 = c(0.9003, -0.5377, 0.2137, -0.0280) MM = mcmc('logprior','loglikelihood','generate','logproposal',m0,N) mout = MM[[1]] mMAP= MM[[2]] pacc= MM[[3]]fun <-function(m,x) { y=m[1]*exp(m[2]*x)+m[3]*x*exp(m[4]*x) return(y) } generate <-function( x) { y=x+step*rnorm(4) return(y) } logprior <-function(m) { if( (m[1]>=0) & (m[1]<=2) & (m[2]>=-0.9) & (m[2]<=0) & (m[3]>=0) & (m[3]<=2) & (m[4]>=-0.9) & (m[4]<=0) ) { lp=0 } else { lp= -Inf } return(lp) } loglikelihood <-function(m) { fvec=(y-fun(m,x))/sigma L=(-1/2)*sum(fvec^2) return(L) } logproposal <-function(x,y) { LR=(-1/2)*sum((x-y)^2/step^2) return(LR) } ### Generate the data set. x=seq(from=1, by=0.25, to=7.0) mtrue=c(1.0, -0.5, 1.0, -0.75) ytrue=fun(mtrue,x) sigma=0.01*rep(1, times= length(ytrue) ) y=ytrue+sigma*rnorm(length(ytrue) ) ### set the MCMC parameters ### number of skips to reduce autocorrelation of models skip=100 ### burn-in steps BURNIN=1000 ### number of posterior distribution samples N=4100 ### MVN step size step = 0.005*rep(1,times=4) ### We assume flat priors here m0 = c(0.9003, -0.5377, 0.2137, -0.0280) alogprior='logprior' aloglikelihood='loglikelihood' agenerate='generate' alogproposal='logproposal' ### ### initialize model at a random point on [-1,1] ### m0=(runif(4)-0.5)*2 ### this is the matlab initialization: m0 = c(0.9003, -0.5377, 0.2137, -0.0280) MM = mcmc('logprior','loglikelihood','generate','logproposal',m0,N) mout = MM[[1]] mMAP= MM[[2]] pacc= MM[[3]]
Matrix Norm
Mnorm(X, k = 2)Mnorm(X, k = 2)
X |
matrix |
k |
norm number |
returns the largest singular value of the matrix or vector
Scalar Norm
if k=1, absolute value; k=2 2-norm (rms); k>2, largest singular value.
Jonathan M. Lees<[email protected]>
x = runif(10) Mnorm(x, k = 2)x = runif(10) Mnorm(x, k = 2)
Number of non-zero elements in a vector
nnz(h)nnz(h)
h |
vector |
integer
Jonathan M. Lees<[email protected]>
zip<-rnorm(15) nnz(zip)zip<-rnorm(15) nnz(zip)
Occam's inversion
occam(afun, ajac, L, d, m0, delta)occam(afun, ajac, L, d, m0, delta)
afun |
character, function handle that computes the forward problem |
ajac |
character, function handle that computes the Jacobian of the forward problem |
L |
regularization matrix |
d |
data that should be fit |
m0 |
guess at the model |
delta |
cutoff to use for the discrepancy principle portion |
vector, model found
This is a simple brute force way to do the line search. Much more sophisticated methods are available. Note: we've restricted the line search to the range from 1.0e-20 to 1. This seems to work well in practice, but might need to be adjusted for a particular problem.
Jonathan M. Lees<[email protected]>
bayes
normal distribution and returns the value of the integral
phi(x)phi(x)
x |
endpoint of integration (scalar) |
value of integral
Jonathan M. Lees<[email protected]>
erf
x <- 1.0 ## pracma::erf(x) phi(x) phiinv( phi(x) )x <- 1.0 ## pracma::erf(x) phi(x) phiinv( phi(x) )
Calculates the inverse normal distribution from the value of the integral
phiinv(x)phiinv(x)
x |
endpoint value of integration (scalar) |
value of integral (scalar)
Jonathan M. Lees<[email protected]>
phi
x <- 1.0 ## pracma::erf(x) phi(x) phiinv( phi(x) )x <- 1.0 ## pracma::erf(x) phi(x) phiinv( phi(x) )
Picard plot parameters for subsequent plotting.
picard_vals(U, sm, d)picard_vals(U, sm, d)
U |
the U matrix from the SVD or GSVD |
sm |
singular values in decreasing order, or the GSVD lambdas divided by the mus in decreasing order |
d |
data to fit, right hand side |
The Picard plot is a method of helping to determine regularization schemes.
List:
utd |
the columns of U transposed times d |
utd_norm |
utd./sm |
Jonathan M. Lees<[email protected]>
GSVD
#### n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 dspiken = G set.seed(2015) dspiken = dspiken + 6e-6 *rnorm(length(dspiken)) Utube=svd(G); U = Utube$u V = Utube$v S = Utube$d s=Utube$d R3 = picard_vals(U,s,dspiken); utd = R3$utd utd_norm= R3$utd_norm ### Produce the Picard plot. x_ind=1:length(s); ## plot( range(x_ind) , range(c(s ,abs(utd),abs(utd_norm))), type='n', log='y', xlab="i", ylab="" ) lines(x_ind,s, col='black') points(x_ind,abs(utd), pch=1, col='red') points(x_ind,abs(utd_norm), pch=2, col='blue') title("Picard Plot for Shaw Problem")#### n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 dspiken = G set.seed(2015) dspiken = dspiken + 6e-6 *rnorm(length(dspiken)) Utube=svd(G); U = Utube$u V = Utube$v S = Utube$d s=Utube$d R3 = picard_vals(U,s,dspiken); utd = R3$utd utd_norm= R3$utd_norm ### Produce the Picard plot. x_ind=1:length(s); ## plot( range(x_ind) , range(c(s ,abs(utd),abs(utd_norm))), type='n', log='y', xlab="i", ylab="" ) lines(x_ind,s, col='black') points(x_ind,abs(utd), pch=1, col='red') points(x_ind,abs(utd_norm), pch=2, col='blue') title("Picard Plot for Shaw Problem")
Plots a model in piecewise constant form over n subintervals, where n is the length of x.
plotconst(x, l, r, ...)plotconst(x, l, r, ...)
x |
model to be plotted |
l |
left endpoint of plot |
r |
right endpoint of plot |
... |
graphical parameters |
Used for plotting vector models
graphical side effects
Jonathan M. Lees<[email protected]>
linesconst
zip = runif(25) plotconst(zip, 0, 1 ) linesconst(runif(25) , 0, 1 , col='red' )zip = runif(25) plotconst(zip, 0, 1 ) linesconst(runif(25) , 0, 1 , col='red' )
Quadratic Linearization
quadlin(Q, A, b)quadlin(Q, A, b)
Q |
positive definite symmetric matrix |
A |
matrix with linearly independent rows |
b |
data vector |
Solves the problem: min (1/2) t(x)*Q*x with Ax = b. using the Lagrange multiplier technique, where Q is assumed to be symmetric and positive definite and the rows of A are linearly independent.
list:
x |
vector of solution values |
lambda |
Lagrange multiplier |
Jonathan M. Lees<[email protected]>
###% Radius of the Earth (km) Re=6370.8; rad = 5000 ri=rad/Re; q=c(1.083221147, 1.757951474) H = matrix(rep(0, 4), ncol=2, nrow=2) H[1,1]=1.508616069 - 3.520104161*ri + 2.112062496*ri^2; H[1,2]=3.173750352 - 7.140938293*ri + 4.080536168*ri^2; H[2,1]=H[1,2]; H[2,2]=7.023621326 - 15.45196692*ri + 8.584426066*ri^2; A1 =quadlin(H,t(q), 1.0 );###% Radius of the Earth (km) Re=6370.8; rad = 5000 ri=rad/Re; q=c(1.083221147, 1.757951474) H = matrix(rep(0, 4), ncol=2, nrow=2) H[1,1]=1.508616069 - 3.520104161*ri + 2.112062496*ri^2; H[1,2]=3.173750352 - 7.140938293*ri + 4.080536168*ri^2; H[2,1]=H[1,2]; H[2,2]=7.023621326 - 15.45196692*ri + 8.584426066*ri^2; A1 =quadlin(H,t(q), 1.0 );
Return the rank of a matrix. Not to be confused with the R function rank.
rnk(G, tol = 1e-14)rnk(G, tol = 1e-14)
G |
Matrix |
tol |
machine tolerance for small numbers |
Number of singular values greater than tol.
integer, number of non-zero singular values
duplicate the matlab function rank
Jonathan M. Lees<[email protected]>
svd
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } X <- hilbert(9)[,1:6] rnk(X)hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } X <- hilbert(9)[,1:6] rnk(X)
Creata design matrix for simulating a tomographic inversion on a simple grid.
setDesignG()setDesignG()
Set up a simple design matrix for tomographic in version. This is used in examples and illustrations of tomographics and matrix inversion methods.
Jonathan M. Lees<[email protected]>
G = setDesignG() ### show the 56-th row g = matrix( G[56,] , ncol=16, nrow=16) imagesc(g) ## Not run: ### show total coverage zim = matrix(0 , ncol=16, nrow=16) for(i in 1:dim(G)[1]) { g = matrix( G[i,] , ncol=16, nrow=16) zim =zim + g } image(zim) ## End(Not run)G = setDesignG() ### show the 56-th row g = matrix( G[56,] , ncol=16, nrow=16) imagesc(g) ## Not run: ### show total coverage zim = matrix(0 , ncol=16, nrow=16) for(i in 1:dim(G)[1]) { g = matrix( G[i,] , ncol=16, nrow=16) zim =zim + g } image(zim) ## End(Not run)
Creates the design matrix for the Shaw inverse problem of diffraction through a narrow slot.
shawG(m, n)shawG(m, n)
m |
integer, number of rows |
n |
integern number of columns |
See Aster's book for a details explaination.
Matrix used for creating data and inversion.
Jonathan M. Lees<[email protected]>
C. B. Shaw, Jr., "Improvements of the resolution of an instrument by numerical solution of an integral equation", J. Math. Anal. Appl. 37: 83-112, 1972.
n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 dspiken = G %*% spike plot(dspiken)n = 20 G = shawG(n,n) spike = rep(0,n) spike[10] = 1 dspiken = G %*% spike plot(dspiken)
Row action method for inversion of matrices, using SIRT algorithm.
sirt(A, b, tolx, maxiter)sirt(A, b, tolx, maxiter)
A |
Design Matrix |
b |
vector, Right hand side |
tolx |
numeric, tolerance for stopping |
maxiter |
integer, Maximum iterations |
Iterates until conversion
Solution vector
Jonathan M. Lees<[email protected]>
Lees, J. M. and R. S. Crosson (1989): Tomographic inversion for three-dimensional velocity structure at Mount St. Helens using earthquake data, J. Geophys. Res., 94(B5), 5716-5728.
art, kac
set.seed(2015) G = setDesignG() ### Setup the true model. mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### reshape the true model to be a vector mtruev=as.vector(mtruem); ### Compute the data. dtrue=G %*% mtruev; ### Add the noise. d=dtrue+0.01*rnorm(length(dtrue)); msirt<-sirt(G,d,0.01,200) par(mfrow=c(1,2)) imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); imagesc(matrix(msirt,16,16) , asp=1 , main="SIRT Solution" );set.seed(2015) G = setDesignG() ### Setup the true model. mtruem=matrix(rep(0, 16*16), ncol=16,nrow=16); mtruem[9,9]=1; mtruem[9,10]=1; mtruem[9,11]=1; mtruem[10,9]=1; mtruem[10,11]=1; mtruem[11,9]=1; mtruem[11,10]=1; mtruem[11,11]=1; mtruem[2,3]=1; mtruem[2,4]=1; mtruem[3,3]=1; mtruem[3,4]=1; ### reshape the true model to be a vector mtruev=as.vector(mtruem); ### Compute the data. dtrue=G %*% mtruev; ### Add the noise. d=dtrue+0.01*rnorm(length(dtrue)); msirt<-sirt(G,d,0.01,200) par(mfrow=c(1,2)) imagesc(matrix(mtruem,16,16) , asp=1 , main="True Model" ); imagesc(matrix(msirt,16,16) , asp=1 , main="SIRT Solution" );
Inverse T-distribution, qt
tinv(p, nu)tinv(p, nu)
p |
P-value |
nu |
degrees of freedom |
Wrapper for qt
Quantile for T-distribution
Jonathan M. Lees<[email protected]>
qt
tinv(.4, 10)tinv(.4, 10)
Singular Value Decomposition
USV(G)USV(G)
G |
Matrix |
returns matrices U, S, V according to matlab convention.
list:
U |
Matrix |
S |
Matrix, singular values |
V |
Matrix |
Jonathan M. Lees<[email protected]>
svd
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } X <- hilbert(9)[,1:6] h = USV(X) print( h$U )hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } X <- hilbert(9)[,1:6] h = USV(X) print( h$U )
Vector 2-Norm.
Vnorm(X)Vnorm(X)
X |
numeric vector |
Numeric scale norm
This function is intended to duplicated the matlab function norm.
Jonathan M. Lees<[email protected]>
V = Vnorm(rnorm(10))V = Vnorm(rnorm(10))
Example vertical 1-dimensional seismic profile used for setting up examples for inverse theory.
vspprofile(M = 50, N = 50, maxdepth = 1000, deltobs = 20, noise = 2e-04, M1 = c(9000, -6, 0.001))vspprofile(M = 50, N = 50, maxdepth = 1000, deltobs = 20, noise = 2e-04, M1 = c(9000, -6, 0.001))
M |
integer, number of rows in in design matrix G, default=50 |
N |
integer, number of columns in design matrix G, default=50 |
maxdepth |
Maximum depth of model, default = 1000 |
deltobs |
integer, sampling interval in depth, default=20 |
noise |
gausian noise multiplier, default=2e-04 |
M1 |
3-vector, linear model for velocity versus depth model |
Vertical seismic profile in 1D dimension used for setting up examples in PEIP. Given a simple velocity profile, defined by input parameter M1 create the travel times and designe matrix used for solving an inverse problem. The velocity model is defined as depth versus velocity, and the function inverts that from the slowness. Any model could be used to replace this model. The default model here is taken from an inversion in the Aster book.
list:
G |
M by N design matrix |
tee |
true travel times from model |
t2 |
travel times with noise added |
depth |
depth samples of model |
vee |
velocity at the depths indicated |
M |
input M |
N |
input N |
maxdepth |
input maxdepth |
deltobs |
input delta observation |
noise |
input noise |
M1 |
True model used for depth versus velocity |
Jonathan M. Lees<[email protected]>
Aster, R.C., C.H. Thurber, and B. Borchers, Parameter Estimation and Inverse Problems, Elsevier Academic Press, Amsterdam, 2005.
V = vspprofile() ### plot quadratic velocity profile plot(V$vee, -V$depth, main="VSP: velocity increasing with depth") dobs = seq(from=V$deltobs, to=V$maxdepth, by=V$deltobs) ### plotdepth versus time (not linear) plot(dobs, V$t2) abline(lm(V$t2 ~ dobs) )V = vspprofile() ### plot quadratic velocity profile plot(V$vee, -V$depth, main="VSP: velocity increasing with depth") dobs = seq(from=V$deltobs, to=V$maxdepth, by=V$deltobs) ### plotdepth versus time (not linear) plot(dobs, V$t2) abline(lm(V$t2 ~ dobs) )