Linear inversion chrestomathy

Software chrestomathy (Lämmel 2015) 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.

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;