# 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. Jan 28, 2019 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? Jan 28, 2019 at 1:14
• I don't know, sounds strange, that needs some investigating. Jan 28, 2019 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. Jan 29, 2019 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!