Thomas,
If you're trying to invert (factorize) a matrix as part of your objective
function, it's generally a much better idea to reformulate your
optimization problem and use a different algorithm. Anywhere that you're
essentially trying to use inv(A)*x, you should try introducing a new
Whoops, didn't finish writing a sentence there... was saying that most of
the algorithms in Optim don't implement equality constraints, but you've
got other choices of optimization solvers.
On Tuesday, June 3, 2014 3:39:44 AM UTC-7, Tony Kelman wrote:
Thomas,
If you're trying to invert
Thanks for the suggestion, Tony. The problem I'm working with is that of
fitting a gaussian process to data, so my objective function (the log
likelihood of a GP given some data x) includes both an x' inv(K(theta)) x term
and a log determinant of K(theta). Here K() is a big dense positive
Alright if K is dense then Ipopt may not be your best choice. How do the
theta parameters enter into K(theta)? Is the problem convex? Would the dual
problem be easier to work with? Have you read Boyd and Vandenberghe more
recently than I have? This sounds suspiciously like it can be expressed
In the GP setting I'm using, K is the set of pair wise gaussian kernel
evaluations between all the points in my data. So if a point is x_i,z_i (x is
the dependent variable and the z's are effectively explanatory variables), K_ij
is more or less exp(-.5 * (z_i - z_j)^2 * theta).
I don't
Hi Thomas,
This is definitely a reasonable thing to want to do, and there are a number
of methods to differentiate through a matrix factorization. Julia doesn't
yet have a generic chol() method
(see https://github.com/JuliaLang/julia/issues/1629), but it does have a
generic LU decomposition
Thanks for the quick response. I'll check and see if an LU decomposition
will do the trick
As far as my own cholesky code, it is quite a bit slower, even after your
edits. For a 1000x1000 Float64 matrix (i.e., X = randn(1000,1000); X = X'
* X), here are the timing results I get:
3 Calls of
