Gm = d in lua/torch for LuaJIT.
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.
We'll start with the example from the torch7 documentation.
In [50]:
G = torch.Tensor({{6.80, -2.11, 5.66, 5.97, 8.23},
{-6.05, -3.30, 5.36, -4.44, 1.08},
{-0.45, 2.58, -2.70, 0.27, 9.04},
{8.32, 2.71, 4.35, -7.17, 2.14},
{-9.67, -5.14, -7.26, 6.08, -6.87}}):t() -- transpose
In [51]:
d = torch.Tensor({{4.02, 6.19, -8.22, -7.57, -3.03},
{-1.56, 4.00, -8.67, 1.75, 2.86},
{9.81, -4.09, -4.57, -8.61, 8.99}}):t()
Later we can use lower-level linear algebraic functions, but for now here's the solution using LAPACK's gesv, which is described like this:
Computes the solution to a real system of linear equations
A*X=B(simple driver). LU decomposition with partial pivoting and row interchanges is used to solve the system.
In [52]:
timer = torch.Timer()
m_est = torch.gesv(d, G)
print(timer:time().real)
Out[52]:
In [53]:
d:dist(G*m_est)
Out[53]:
In [54]:
-- pseudo-inverse
torch.inverse(G:t()*G) * G:t()
Out[54]:
In [55]:
-- Define model
M = 50
m = torch.zeros(M, 1)
m[{{11,15}, {}}] = 1.0
m[{{16,27}, {}}] = -0.3
m[{{28,35}, {}}] = 2.1
In [56]:
Plot = require 'itorch.Plot'
In [57]:
x = torch.range(1,m:size(1))
plot = Plot():line(x, m[{{},1}]):circle(x, m[{{},1}]):title('Model m'):draw()
In [58]:
-- Discrete kernel G
N = 20
L = 100
alpha = 0.8
In [59]:
M = 50
x = torch.range(0, M, 1) * L/(M-1)
dx = L/(M-1)
r = torch.range(0, N, 1) * L/(N-1)
In [60]:
G = torch.zeros(N, M)
In [61]:
for j = 1,M do
for k = 1,N do
t = torch.abs(r[{k}] - x[{j}])
G[{{k},{j}}] = dx * torch.exp(-alpha * torch.pow(t, 2))
end
end
In [62]:
d = G * m
In [63]:
Plot = require 'itorch.Plot'
x = torch.range(1,d:size(1))
plot = Plot():line(x, d[{{},1}]):circle(x, d[{{},1}]):title('Data d'):draw()
In [64]:
-- Compute data
d = G * m
In [65]:
-- Minimum norm solution
m_est = G:t() * torch.inverse(G * G:t()) * d
d_pred = G * m_est
In [66]:
-- Plot as before
Plot = require 'itorch.Plot'
x = torch.range(1,d_pred:size(1))
plot = Plot():line(x, d_pred[{{},1}]):circle(x, d_pred[{{},1}]):title('Predicted data d_pred'):draw()
In [67]:
-- Use torch.gels()
m_est = torch.gels(d, G)
d_pred = G * m_est
-- Plot
Plot = require 'itorch.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 [68]:
d:dist(G * m_est)
Out[68]:
In [69]:
-- Add noise
dc = G * m
s = 1
d = dc + s * torch.DoubleTensor(d:size()):uniform(0,1)
In [21]:
-- Compute using the second form
I = torch.eye(N)
mu = 2.5
m_est = G:t() * torch.inverse(G*G:t() + mu*I) * d
d_pred = G * m_est
In [22]:
-- Plot
Plot = require 'itorch.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()
It's not a good idea to use the matrix inverse. So let's try factorization:
In [23]:
chol = torch.potrf(G * G:t() + mu * I)
In [24]:
m_est = torch.potrs(G:t() * d, chol)
In [25]:
d_pred = G * m_est
In [26]:
-- Plot
Plot = require 'itorch.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()
OK, I guess I screwed that up.
We need to implement convmtx, the matrix equivalent of the convolutional operator. This helps us make the first derivative.
There's no Toeplitz helper function that I can find, so we'll have to do it a slightly different way from the one we used in the NumPy notebook.
In [27]:
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 [28]:
h = torch.Tensor({1,-1})
n = 8 -- number of rows in conv. matrix
In [29]:
convmtx(h, n)
Out[29]:
As far as I can tell, there's no simple way to index off the last column, as with the -1 trick in Python. But we know the number of columns is 8 + h.size()[1] - 1 so we can just substract one from that...
In [30]:
cols = n + h:size()[1] - 1
print(convmtx(h, n)[{ {},{1, cols-1} }])
Out[30]:
In [31]:
cols = M + h:size()[1] - 1
W = convmtx(h, M)[{ {},{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 [32]:
m_est = torch.inverse(G:t()*G + mu*W:t()*W) * G:t() * d
d_pred = G * m_est
In [33]:
-- Plot
Plot = require 'itorch.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 [34]:
-- Plot
Plot = require 'itorch.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()
From the iTorch readme.
In [35]:
x1 = torch.randn(40):mul(100)
y1 = torch.randn(40):mul(100)
x2 = torch.randn(40):mul(100)
y2 = torch.randn(40):mul(100)
x3 = torch.randn(40):mul(200)
y3 = torch.randn(40):mul(200)
-- scatter plots
plot = Plot():circle(x1, y1, 'red', 'hi'):circle(x2, y2, 'blue', 'bye'):draw()
plot:circle(x3,y3,'green', 'yolo'):redraw()
plot:title('Scatter Plot Demo'):redraw()
plot:xaxis('length'):yaxis('width'):redraw()
plot:legend(true)
plot:redraw()
In [36]:
batchSize, inputSize, outputSize = 4000, 2000, 3000
input = torch.FloatTensor(batchSize, inputSize):uniform(0,1)
weight = torch.FloatTensor(outputSize, inputSize):uniform(0,1)
output = torch.FloatTensor(batchSize, outputSize) -- Must give size here, AFAICT
In [37]:
timer = torch.Timer()
output:addmm(0, output, 1, input, weight:t()) -- transpose of weight
print(timer:time().real)
Out[37]:
How about CudaTensors?
Note that you need an nVidia graphics card with a recent version of CUDA installed, and you will need to install cutorch with luarocks.
In [38]:
require 'cutorch'
input = torch.CudaTensor(batchSize, inputSize):uniform(0,1)
weight = torch.CudaTensor(outputSize, inputSize):uniform(0,1)
output = torch.CudaTensor(batchSize, outputSize)
In [39]:
time = timer:reset()
output:addmm(0, output, 1, input, weight:t())
print(timer:time().real)
Out[39]:
In [35]:
-- output