iTorch / Torch7 implementation

Gm = d in lua/torch for LuaJIT.

Prerequisites

You need to install torch and (if you want to run this notebook) iTorch. Works for me on Ubuntu 16.04 and MacOS 10.11. Not gonna say the Mac installation is easy though. Everything was very smooth on Ubuntu.



In [1]:

    
Plot = require 'itorch.Plot'



In [2]:

    
-- Implements MATLAB's convmtx.
function convmtx (h, n)
    local m = h:size()[1] - 1
    local c = torch.Tensor(n, n+m):fill(0)
    for i = 1, n do
        c[{ i, {i, i+m} }] = h
    end
    return c
end



In [22]:

    
convmtx(torch.Tensor({1, -1}), 5)









    Out[22]:





 1 -1  0  0  0  0
 0  1 -1  0  0  0
 0  0  1 -1  0  0
 0  0  0  1 -1  0
 0  0  0  0  1 -1
[torch.DoubleTensor of size 5x6]



In [3]:

    
-- Implements bruges's ricker.
function ricker (f, length, dt)
    -- Defaults.
    f = f or 25
    length = length or 0.128
    dt = dt or 0.001

    -- Time basis.
    local t = torch.range(0, length, dt) - length/2
    
    -- Amplitudes.
    local x = 1 - 2 * math.pi * f^2 * torch.pow(t, 2)
    local y = torch.exp(-math.pi^2 * f^2 * torch.pow(t, 2))
    local a = torch.cmul(x, y)

    return t, a
end

Construct the model



In [4]:

    
--                         VP    RHO
imp = torch.ones(51, 1) * 2550 * 2650
imp[{{11, 15}, {}}] =     2700 * 2750   
imp[{{16, 27}, {}}] =     2400 * 2450
imp[{{28, 35}, {}}] =     2800 * 3000



In [5]:

    
x = torch.range(1, imp:size(1))
plot = Plot():line(x, imp[{{},1}])
             :title('Impedance')
             :draw()



In [6]:

    
numerator = imp[{{2, imp:size(1)}}] - imp[{{1, imp:size(1) - 1}}]
denominator = imp[{{2, imp:size(1)}}] + imp[{{1, imp:size(1) - 1}}]
m = torch.cdiv(numerator, denominator)



In [7]:

    
x = torch.range(1, m:size(1))
plot = Plot():line(x, m[{{}, 1}], 'blue')
             :title('Reflectivity')
             :draw()

Forward operator: convolution with wavelet

Now we make the kernel matrix G, which represents convolution.



In [35]:

    
t, wavelet = ricker(100, 0.02, 0.001)



In [36]:

    
plot = Plot():line(t, wavelet)
             :title('Wavelet')
             :draw()



In [89]:

    
G_ = convmtx(wavelet, m:size(1))[{{}, {11,60}}]



In [83]:

    
-- If anyone knows a better way to get every other row of G,
-- I'd love to hear it!
G = torch.Tensor(G_:size(1)/2, G_:size(2))
for i = 1, G_:size(1) do
    if i % 2 == 1 then
        G[{(i+1)/2}] = G_[{i}]
    end
end



In [84]:

    
d = G * m



In [85]:

    
x = torch.range(1, d:size(1))
plot = Plot():circle(x, d[{{}, 1}])
             :title('Data d')
             :draw()

Noise free: minimum norm



In [86]:

    
m_est = G:t() * torch.inverse(G * G:t()) * d
d_pred = G * m_est



In [87]:

    
x = torch.range(1, m_est:size(1))
plot = Plot():line(x, m_est[{{}, 1}])
             :title('Predicted model m_est')
             :draw()



In [46]:

    
x = torch.range(1, d_pred:size(1))
plot = Plot():circle(x, d_pred[{{}, 1}])
             :title('Predicted data d_pred')
             :draw()

Least squares solution using GESV

This uses LAPACK. It proceeds via factorization instead of the inverse of (G * G.T). It's the equivalent of calling NumPy or SciPy's linalg.solve function.



In [90]:

    
-- Use torch.gels()
m_est = torch.gels(d, G)
d_pred = G * m_est

-- Plot
x = torch.range(1,d_pred:size(1))
plot = Plot():circle(x, d[{{},1}], 'blue')
             :circle(x, d_pred[{{},1}], 'red')
             :title('blue:d   red:d_pred')
             :draw()



In [93]:

    
x = torch.range(1, m_est:size(1))
plot = Plot():line(x, m_est[{{}, 1}])
             :title('Predicted model m_est')
             :draw()



In [91]:

    
d:dist(d_pred)









    Out[91]:





9.978417090315e-17

This is the L2 norm; Mauricio used the square of this as the misfit:



In [92]:

    
(d_pred - d):t() * (d_pred - d)









    Out[92]:





1e-33 *
  9.9569
[torch.DoubleTensor of size 1x1]

QR factorization

Again, avoiding finding inverses. Uses LAPACK.

As with any QR factorization, G must be square.



In [ ]:

    
q, r = torch.qr(G)
p = q:t() * d
m_est = torch.inverse(r) * p
d_pred = G * m_est

-- Plot
x = torch.range(1,d_pred:size(1))
plot = Plot():circle(x, d[{{},1}], 'blue')
             :circle(x, d_pred[{{},1}], 'red')
             :title('blue:d   red:d_pred')
             :draw()



In [ ]:

    
d:dist(d_pred)

With noise: damped least squares



In [52]:

    
-- Add noise
dc = G * m
s = 0.025
d = dc + s * (torch.DoubleTensor(d:size()):uniform(0,1) - 0.5)



In [53]:

    
-- Compute using Mauricio's second form
I = torch.eye(G:size(1))
mu = 2.5
m_est = G:t() * torch.inverse(G*G:t() + mu*I) * d
d_pred = G * m_est



In [54]:

    
-- Plot
x = torch.range(1,d_pred:size(1))
plot = Plot():line(x, d[{{},1}], 'blue'):circle(x, d[{{},1}], 'blue')
             :line(x, d_pred[{{},1}]):circle(x, d_pred[{{},1}])
             :title('blue:d   red:d_pred')
             :draw()

Try using Cholesky decomposition

It's not a good idea to use the matrix inverse. So let's try factorization:



In [55]:

    
chol = torch.potrf(G * G:t() + mu * I)



In [56]:

    
m_est = torch.potrs(G:t() * d, chol)



In [57]:

    
d_pred = G * m_est



In [58]:

    
-- Plot
x = torch.range(1,d_pred:size(1))
plot = Plot():line(x, d[{{},1}], 'blue'):circle(x, d[{{},1}], 'blue')
             :line(x, d_pred[{{},1}]):circle(x, d_pred[{{},1}])
             :title('blue:d   red:d_pred')
             :draw()

With noise: damped LS with 1st derivative regularization



In [59]:

    
h = torch.Tensor({1,-1})
cols = G:size(1) + h:size()[1] - 1
W = convmtx(h, G:size(1))[{ {},{1, cols-1} }]

Now we solve:

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



In [60]:

    
m_est = torch.inverse(G:t()*G + mu*W:t()*W) * G:t() * d
d_pred = G * m_est



In [61]:

    
-- Plot
x = torch.range(1,d_pred:size(1))
plot = Plot():line(x, d[{{},1}], 'blue'):circle(x, d[{{},1}], 'blue')
             :line(x, d_pred[{{},1}]):circle(x, d_pred[{{},1}])
             :title('blue:d   red:d_pred')
             :draw()



In [62]:

    
-- Plot
x = torch.range(1,m_est:size(1))
plot = Plot():line(x, m[{{},1}], 'blue', 'blue')
             :line(x, m_est[{{},1}]):title('blue:m   red:m_est')
             :draw()

Damped LS with 2nd derivative regularization



In [63]:

    
h = torch.Tensor({1,-2, 1})
W = convmtx(h, G:size(1))[{{}, {2, G:size(1)+1}}]

Q = torch.inverse(W * W:t())
m_est = Q * G:t() * torch.inverse(G * Q * G:t()) * d
d_pred = G * m_est

-- Plot
x = torch.range(1, m_est:size(1))
plot = Plot():line(x, m[{{}, 1}], 'blue', 'blue')
             :line(x, m_est[{{}, 1}]):title('blue:m   red:m_est')
             :draw()



In [ ]: