# Linear inversion chrestomathy

Jump to navigation
Jump to search

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

## Contents

## The problem

Solve

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

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

- ↑ 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.
- ↑ 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.