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 variable y and corresponding linear equality constraint x = A * y. Most of the algorithms in Optim If your objective function in terms of x and y is linear, then you can use Clp, GLPK, etc. If it's still some complicated nonlinear thing even in terms of x and y, use JuMP with Ipopt. NLopt might work too, but it's not tied into the convenient JuMP reverse-mode autodiff framework yet.
-Tony On Monday, June 2, 2014 11:43:32 AM UTC-7, Thomas Covert wrote: > > I'm really excited to use DualNumbers and the autodiff option in Optim.jl. > However, the functions I would like to optimize often make use of linear > algebra, including cholesky factorizations. At present, DualNumbers.jl > does not implement a cholesky factorization. My first pass at a workaround > for this was to just write out a standard cholesky solver (i.e., just > writing out the equations listed in wikipedia). This seems to work fine - > i.e., my "choldn()" function factors a DualNumber array into upper and > lower triangular matrices whose product is the input. > > However, my "choldn" function listed below is MUCH slower than the LAPACK > call built into chol(). I can't say I'm surprised by this, but the speed > hit makes the whole enterprise of using autodiff with DualNumbers > considerably less useful. > > I was hoping to find some neat linear algebra trick that would let me > compute a DualNumber cholesky factorization without having to resort to > non-LAPACK code, but I haven't found it yet. That is, I figured that I > could compute the cholesky factorization of the real part of the matrix and > then separately compute a matrix which represented the DualNumber part. > However, I've not yet found a math text describing such a beast. > > Any ideas? Am I wrong in thinking that a Cholesky decomposition of a > square array with DualNumbers is a well-defined mathematical object? > > I've included my code for "choldn()" below. > > Thanks in advance. > > > CODE > > function choldn!(A::Array) > > N = size(A,1); > > T = eltype(A); > > tmp = zero(T); > > for j = 1:N # cols > > > for i = j:N # rows > > > if i == j > > if i == 1 > > A[i,j] = sqrt(A[i,j]); > > else > > tmp = zero(T); > > for k = 1:j-1 > > tmp += A[j,k] ^ 2; > > end > > A[i,j] = sqrt(A[i,j] - tmp); > > end > > else > > if j == 1 > > A[i,j] = A[i,j] / A[j,j]; > > else > > tmp = zero(T); > > for k = 1:j-1 > > tmp += A[i,k] * A[j,k]; > > end > > A[i,j] = (A[i,j] - tmp) / A[j,j]; > > end > > end > > end > > end > > # then zero out the non-cholesky stuff > > > for i = 1:N > > if i < N > > for j = i+1:N > > A[i,j] = zero(T); > > end > > end > > end > > end > >
