torch.manualSeed(1234)

-- choose a dimension
N = 5

-- create a random NxN matrix
A = torch.rand(N, N)

-- make it symmetric positive
A = A*A:t()

-- make it definite
A:add(0.001, torch.eye(N))

-- add a linear term
b = torch.rand(N)

-- create the quadratic form
function J(x)
  return 0.5*x:dot(A*x)-b:dot(x)
end

-- print the function value
print(J(torch.rand(N)) )

-- inverse the matrix
xs = torch.inverse(A)*b
print(string.format('J(x^*) = %g', J(xs)))


-- define the gradient
function dJ(x)
   return A*x-b
end

-- define current solution
x = torch.rand(N)


-- apply gradient descent
lr = 0.01
for i=1,20000 do
  x = x - dJ(x)*lr
  -- we print the value of the objective function at each iteration
  print(string.format('at iter %d J(x) = %f', i, J(x)))
end