Linear inversion chrestomathy

From SubSurfWiki
Jump to: navigation, search

Software chrestomathy (Lämmel 2015[1]) for linear inversion.

The problem

Solve

 \mathbf{G}{m} = {d}

using the minimum norm solution for the model estimate m-hat given by

 \hat{m} = \mathbf{G}^\mathrm{T} (\mathbf{G} \mathbf{G}^\mathrm{T})^{-1} d

where T denotes the transpose.

The challenge is to implement the solution in as many programming languages as you can think of... but at least Python (with NumPy I think), Julia, R, MATLAB, Java(?), Javascript(??), C, C++, and FORTRAN.

MATLAB

From Mauricio Sacchi's course notes.[2]

Forward problem

Define a model and kernel, and generate data.

% Define model m
M = 50;
m = zeros(M,1);
m(10:14,1) = 1.;
m(15:26,1) = -.3;
m(27:34,1) = 2.1;
 
% Discrete kernel G
N = 20;
L = 100;
alpha = .8
x = (0:1:M-1)*L/(M-1);
dx = L/(M-1);
r = (0:1:N-1)*L/(N-1);
for j=1:M
for k=1:N
G(k,j) = dx*exp(-alpha*abs(r(k)-x(j))^2);
end
end
 
% Compute data
d = G*m;

Inverse problem

The minimum norm solution is given by:

function [m_est,d_pred] = min_norm_sol(G,d);
%
% Given a discrete kernel G and the data d, we compute the
% minimum norm solution of the inverse problem d = Gm.
% m_est: solution (minimum norm solution)
% d_pred: predicted data
%
m_est = G’*inv(G*G’)*d;
d_pred = G* m_est;

References

  1. Lämmel, Ralf (2015). Software chrestomathies. Science of Computer Programming, Volume 97, Part 1, 1 January 2015, Pages 98–104, Special Issue on New Ideas and Emerging Results in Understanding Software. DOI:10.1016/j.scico.2013.11.014. Available online.
  2. Sacchi, MD. GEOPH431 and GEOPH531: Notes on linear inverse theory, minimum norm solutions, and quadratic regularization with MATLAB examples, PART 1. Physics, University of Alberta.

External links