CUTEst returns wrong sparsity patterns for a large number of lower-dimensional problems

[johannahaffner]

For a large number of lower-dimensional problems I have tested, the sparsity patterns returned from CUTEst do seem to be incorrect. I have not checked for the high-dimensional problems yet. Verification was done analytically/by taking the derivatives by hand.

I first raised this issue on the pycutest repository: jfowkes/pycutest#105.

Our use case is benchmarking of a sparsity detection tool, and we were hoping to be able to do that on CUTEst problems.

If this issue persists at scale, do you think that there is an option to run a sparsity detection algorithm on the CUTEst collection, perhaps the one currently implemented in Casadi? That one seems to be battle-tested enough that it could serve as a ground truth.

[nimgould]

You are correct, all that CUTEst currently does is provide the required
(objective or constraint) Hessian within the sparsity pattern of the
Hessian of the Lagrangian. Clearly this is not optimal, but with extra
programming it would be easy to do this properly. For HS88, the Hessian of the Lagrangian is dense (because the constraint Hessian is), and hence what you see.

I'll add it to the "to do" list

[nimgould]

Just a quick update here. There are two issues:

1. In the constrained case, Hessians of individual functions are currently returned within the sparsity pattern of the Hessian of the Lagrangian. This was done for expediency, and is clearly bad in some cases (such as HS92 as you mention). Since this a structural issue, I believe that I can do this better, but it will take a careful bit of reprogramming, so it may take a few days.

2. The second issue relates to numerical, not structural, values in the way that CUTEst builds derivatives from a group-partially-separable representation of a function. If you recall, a (group) function is expressed as f_j(x) = g_j( sum_i e_ij(x) ) where each g_i is a scalar function, and each e_ij is a function of a few of the components of x (or more generally is a low-dimensional range transformation of x). The structure tells us how f_j is composed, but says nothing about the e_ij (other than which variables are involved). CUTEst assumes that derivatives of e_ij involve all of its involved components, but structurally knows no more (only a compiler or some other more complicated description could know more). Thus if e_ij = x_1 x_2 (for example), the 1,1 and 2,2 second derivatives of e_ij are always zero, but structurally we cannot know this. Indeed if e_ij = x^a y^b for any a, b /= 0, the element Hessian would be dense. This is beyond what CUTEst was designed to deal with, and of course what you see for HATFLDH in which there are two elements x_1 x_3 and x_2 x_4. I will not be able to do more here without a complete rewrite of the SIF decoder, and I'll leave this for someone else if they really think it is important.

I'll let you know when I have dealt with 1. above

[nimgould]

@jfowkes I have updated the fortran cish function to record a better representation of individual Hessians for constrained problems (issue 1 above). Could you see if this does the job with pycutest, please?

I've also added a new fortran function, cishp, that records just the patterns that correspond to the matrices found by cish; this is the analog of ushp and cshp that do the same for unconstrained Hessians, and constrained Lagrangian Hessians found by ush and csh. I don't think you have python versions of these, do you? These are possibly more useful for @johannahaffner

[jfowkes]

Thanks @nimgould, I will test this with PyCUTEst's isphess function and report back.

No we don't currently have interfaces to ushp and cshp, the primary design goal of PyCUTEst was to replicate the CUTEst Matlab interface as closely as possible and these aren't exposed in the Matlab interface. I am happy to add interfaces to these together with cishp if there is demand, however this will take some effort as it requires writing Python C API wrappers.

[nimgould]

OK, @jfowkes of course determining demand depends on users knowing that this is a possibility ...

I should add that I suspect that some of the older SIF examples (including some that I wrote) could have been better written to expose sparsity. Oh well

[amontoison]

@nimgould @jfowkes
Be careful when you add new files to not forget to also update the Meson build system.
The new file was not compiled. I will fix it.
Do you need a new release?

[nimgould]

Ooops, sorry, @amontoison , I did this in a bit of a rush between other projects, I need to slow down! Let's wait for @jfowkes to test the
Python version before we issue a new release. Many thanks for the updates.

[jfowkes]

@nimgould @amontoison our automated PyCUTEst tests ran on Sunday 18th and everything passed:
https://github.com/jfowkes/pycutest/actions
Our CI tests clone the latest CUTEst release from GitHub and I believe the changes were committed before then. I don't yet interface to cishp so this was not an issue for me. Feel free to do a new release.

[amontoison]

@nimgould
I have a few random errors with cish in a PR of CUTEst.jl.
The platforms are Windows and linux ARM.
From what I understand it seems that the memory randomly contain junk values.
Maybe it is similar to this issue? #63

[nimgould]

Do you have a specific cutest example that it fails on that I can test. I don't think that it is the cigr issue, as the cish code has the correct call. But cish was "improved" fairly recently, so probably I added a bug.

[amontoison]

I can reproduce the error on my laptop with the problem "HS14".

I use x = [10, -10], and the content of the vector containing the nonzeros H_val is not filled when calling cish with the first constraint.

I always get [-0.5, -2.0] with the second constraint.

Maybe the first constraint is linear, but I would expect cish to fill the vector with 0.

If I call cish on the first constraint after the second constraint, the vector containing the nonzeros of the j-th Hessian is still [-0.5, -2.0].

My value of lh is the number of nonzeros for the full Hessian of the Lagrangian.
The vectors H_val, H_row and H_col have length lh and I expect cish to just fill correctly the vector H_val.

[nimgould]

I tested the fortran version using three compilers (gfortran, nagfor, ifort) with debug testing on this example, calling for the Hessian of each constraint in turn, twice, and get

CALL CUTEST_cish for constraint 1, H_ne = 0
CALL CUTEST_cish for constraint 2, H_ne = 2

• H(sparse)

• row    col     val        row    col     val    
  

•   1      1 -5.0000E-01      2      2 -2.0000E+00
  

CALL CUTEST_cish for constraint 1, H_ne = 0
CALL CUTEST_cish for constraint 2, H_ne = 2

• H(sparse)

• row    col     val        row    col     val    
  

•   1      1 -5.0000E-01      2      2 -2.0000E+00
  

which is what I would expect (the second constraint is quadratic, so its Hessian is (diagonal and) constant, while the first is linear, and its Hessian is zero).

Thus I suspect that the issue may lie either in the C interface to cish or in the Julia interface to the C

[nimgould]

Sorry, I just read the above more carefully. Since for the first constraint, nnzh = 0, I am not sure I follow what "fill the vector with 0" when the vector of nonzeros is of length zero. I presume you are checking nnzh before you try to access its nonzeros?

[jfowkes]

@amontoison I'm afraid that this is a bug in the Julia interface, the Python interface (which also calls the C interface) works as expected:

$ python test_bug.py 
H1 <COOrdinate sparse matrix of dtype 'float64'
	with 2 stored elements and shape (2, 2)>
  Coords	Values
  (0, 0)	-0.5
  (1, 1)	-2.0
H2 <COOrdinate sparse matrix of dtype 'float64'
	with 0 stored elements and shape (2, 2)>
H1 <COOrdinate sparse matrix of dtype 'float64'
	with 2 stored elements and shape (2, 2)>
  Coords	Values
  (0, 0)	-0.5
  (1, 1)	-2.0
H2 <COOrdinate sparse matrix of dtype 'float64'
	with 0 stored elements and shape (2, 2)>

Here is the test script that I used:

import numpy as np
import pycutest

p = pycutest.import_problem('HS14')

x = np.array([10,-10])

H1 = p.isphess(x, 0)
print('H1', H1)

H2 = p.isphess(x, 1)
print('H2', H2)

H1 = p.isphess(x, 0)
print('H1', H1)

H2 = p.isphess(x, 1)
print('H2', H2)

note that isphess calls CUTEst_cish.