4
$\begingroup$

The problem I am trying to solve involves minimising a function with respect to a large number (probably 10,000+) of parameters. I can cheaply compute both its Jacobian and its Hessian. The Hessian is very helpful; I am not entirely sure the function is convex, but using the Newton-CG algorithm in scipy.optimize:

result = minimize(to_minimize, f, jac=jac, hess=hess, method='Newton-CG')

converges quickly.

So far so good, but unfortunately, scipy.optimize.minimize is unable to deal with sparse Hessians, so I have to convert the extremely sparse Hessian to a dense matrix. Internally, I believe the Newton-CG method multiplies that Hessian with other things, which I assume will also be much slower than using its sparseness.

In short, my question is: does anyone know of a library that can exploit the sparse Hessian, preferably using the Newton-CG algorithm, or do I have to write one myself? It would be nice if it were in python, but I'm happy to use C++ if that is more fruitful.

Improve this question
$\endgroup$
4
  • $\begingroup$ Use the hessp argument (docs.scipy.org/doc/scipy/reference/generated/…) to compute the sparse matrix-vector product $x\mapsto Hx$ yourself? Then it doesn't need to know or care that the Hessian is sparse. $\endgroup$
    – Kirill
    Jan 28, 2019 at 1:05
  • $\begingroup$ Thanks Kirill, that sounds good, but for some reason, providing only the hessp seems to make the algorithm much slower than giving the full dense Hessian. Is that expected? $\endgroup$
    – noctilux
    Jan 28, 2019 at 1:14
  • $\begingroup$ I don't know, sounds strange, that needs some investigating. $\endgroup$
    – Kirill
    Jan 28, 2019 at 1:54
  • $\begingroup$ It might not be all that helpful, but for sparse direct solvers, nothing really comes close to the stuff that Tim Davis and others have packed into the backend of MATLAB. $\endgroup$
    – whpowell96
    Jan 29, 2019 at 2:54

1 Answer 1

Reset to default
1
$\begingroup$

Thanks for everyone's comments!

I managed to solve this for my problem by writing my own Newton optimisation:

import numpy as np
from sksparse.cholmod import cholesky

f = np.zeros(n_data)

difference = 999

while difference > 1e-5:

    cur_hess = hess(f)
    cur_jac = jac(f)

    cur_chol = cholesky(cur_hess)
    sol = cur_chol.solve_A(cur_jac)

    new_f = f - sol
    difference = np.linalg.norm(f - new_f)

    f = new_f

I'm using scikit-sparse to do a sparse Cholesky decomposition with CHOLMOD here. This works really well for my problem and typically converges within a few iterations. Hope it's of use to some!

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.