# Minimize a function with sparse Hessian

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.

• 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. – Kirill Jan 28 at 1:05
• 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? – noctilux Jan 28 at 1:14
• I don't know, sounds strange, that needs some investigating. – Kirill Jan 28 at 1:54
• 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. – whpowell96 Jan 29 at 2:54

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!