xuy/pyipopt

Same return values (of eval_jac_g) in python yield different results in IPOPT

goett2win opened this issue · 3 comments

Hi,
I noticed that I get different results from IPOPT when returning the very same arrays (in terms of numpy.array_equal(a,b)) as sparsity structure for the constraint Jacobian.

This is the code:

matrix = np.array([0.,0.,0.,0.,0.,0.,1.,1.,1.,0.],
                 [0.,0.,0.,1.,1.,0.,0.,1.,0.,0.])

matrix_sparse = matrix[np.nonzero(matrix)]

def eval_jac_g(x, flag, user_data=None):
    if flag:
        rows1, cols1 = numpy.nonzero(matrix)
        cols2 = numpy.array([6,7,8,3,4,7])
        rows2 = numpy.array([0,0,0,1,1,1])
        print numpy.array_equal(rows1, rows2) # returns True
        print numpy.array_equal(cols1, cols2) # returns True
        return (rows1, cols1)                   # or (rows2, cols2)
    else:
        return matrix_sparse

I checked that np.array_equal(rows1,rows2) == True, same for cols. However, if I run my code, then I get different results: Assigning the arrays manually (rows1, cols1) gives me a successful run of IPOPT whereas (rows2, cols2) gives an "The index of a matrix is out of range." error from MA27.

I am quite clueless about how this can happen. Any ideas?
I can provide a full example if necessary.

xuy commented

Hey @goett2win , I think there is a misunderstanding on the return values of eval_jac_g.
The flag used for eval_jac_g, as the document says, is:

           if the flag is true, it supposed to return a tuple (row, col)
                   to indicate the sparse Jacobi matrix's structure.
           if the flag is false if returns the values of the Jacobi matrix 
                   with length nnzj.

In other words, you should not return matrix_sparse when flag is false. Instead, you should return the actual values of the Jacobi matrix, as a one dimension array.

Thanks for the reply, @xuy . However, matrix_sparse does indeed contain "the actual values of the Jacobi matrix, as a one dimension array"; it is equal to numpy.array([1., 1., 1., 1., 1., 1.]). So that return value should be fine, right?

Hi, for completeness, I added the IPOPT output on highest verbosity level. There you can see that there are nonzeros in the inequality constraint Jacobian, which should not be the case, and the "matrix index out of range" error (both emphasized).

[PyIPOPT] Ipopt will use Hessian approximation.

[PyIPOPT] Problem created
x0 = [ 0.12122838 0.11242599 0.17030618 -0.17910858 -0.05394987 -0.05788019
0.16637586 0.1251587 0.29153457 0.00393032]
nvar, ncon, nnzj = 10 2 6
g_L, g_U = [ 0. 0.] [ 0. 0.]
f(x0) = 0.111517094123
g(x0) = [ 0. 0.]

List of options:

                                Name   Value                # times used
               hessian_approximation = limited-memory            6
                         print_level = 11                        2

This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt


This is Ipopt version 3.11, running with linear solver ma27.

[Callback:E] eval_jac_g

/home/kheidem1/git/phd/NSN/directed_graph.py(78)eval_jac_g()
-> return (Jrows, Jcols)
(Pdb) [Callback:R] eval_jac_g(1)
[Callback:R] eval_jac_g
Number of nonzeros in equality constraint Jacobian...: 6
Number of nonzeros in inequality constraint Jacobian.: 3
Number of nonzeros in Lagrangian Hessian.............: 0

Hessian approximation will be done in the space of all 10 x variables.

[Callback:E] eval_grad_f
[Callback:R] eval_grad_f
Scaling parameter for objective function = 1.000000e+00
[Callback:E] eval_jac_g
[Callback:R] eval_jac_g(2)
[Callback:R] eval_jac_g
objective scaling factor = 1
No x scaling provided
No c scaling provided
No d scaling provided
DenseVector "original x_L unscaled" with 0 elements:
DenseVector "original x_U unscaled" with 0 elements:
DenseVector "original d_L unscaled" with 0 elements:
DenseVector "original d_U unscaled" with 0 elements:
DenseVector "modified x_L scaled" with 0 elements:
DenseVector "modified x_U scaled" with 0 elements:
DenseVector "modified d_L scaled" with 0 elements:
DenseVector "modified d_U scaled" with 0 elements:
DenseVector "initial x unscaled" with 10 elements:
initial x unscaled[ 1]= 1.2122838365799216e-01
initial x unscaled[ 2]= 1.1242599030077560e-01
initial x unscaled[ 3]= 1.7030618337552261e-01
initial x unscaled[ 4]=-1.7910857673273917e-01
initial x unscaled[ 5]=-5.3949872144319122e-02
initial x unscaled[ 6]=-5.7880193074747011e-02
initial x unscaled[ 7]= 1.6637586244509472e-01
initial x unscaled[ 8]= 1.2515870458842004e-01
initial x unscaled[ 9]= 2.9153456703351477e-01
initial x unscaled[ 10]= 3.9303209304278885e-03
Initial values of x sufficiently inside the bounds.
Initial values of s sufficiently inside the bounds.
[Callback:E] eval_jac_g
[Callback:R] eval_jac_g(2)
[Callback:R] eval_jac_g
[Callback:E] eval_grad_f
[Callback:R] eval_grad_f

CompoundVector "RHS[ 0]" with 4 components:

Component 1:
DenseVector "RHS[ 0][ 0]" with 10 elements:
RHS[ 0][ 0][ 1]=-1.2122838365799216e-01
RHS[ 0][ 0][ 2]=-1.1242599030077560e-01
RHS[ 0][ 0][ 3]=-1.7030618337552261e-01
RHS[ 0][ 0][ 4]= 1.7910857673273917e-01
RHS[ 0][ 0][ 5]= 5.3949872144319122e-02
RHS[ 0][ 0][ 6]= 5.7880193074747011e-02
RHS[ 0][ 0][ 7]=-1.6637586244509472e-01
RHS[ 0][ 0][ 8]=-1.2515870458842004e-01
RHS[ 0][ 0][ 9]=-2.9153456703351477e-01
RHS[ 0][ 0][ 10]=-3.9303209304278885e-03

Component 2:
DenseVector "RHS[ 0][ 1]" with 0 elements:

Component 3:
DenseVector "RHS[ 0][ 2]" with 2 elements:
Homogeneous vector, all elements have value 0.0000000000000000e+00

Component 4:
DenseVector "RHS[ 0][ 3]" with 0 elements:
Homogeneous vector, all elements have value 0.0000000000000000e+00

CompoundSymMatrix "KKT" with 4 rows and columns components:
Component for row 0 and column 0:

SumSymMatrix "KKT[0][0]" of dimension 10 with 2 terms:
Term 0 with factor 0.0000000000000000e+00 and the following matrix:

DiagMatrix "Term: 0" with 10 rows and columns, and with diagonal elements:
  DenseVector "Term: 0" with 10 elements:
  Homogeneous vector, all elements have value  0.0000000000000000e+00

Term 1 with factor 1.0000000000000000e+00 and the following matrix:

DiagMatrix "Term: 1" with 10 rows and columns, and with diagonal elements:
  DenseVector "Term: 1" with 10 elements:
  Homogeneous vector, all elements have value  1.0000000000000000e+00

Component for row 1 and column 0:
This component has not been set.
Component for row 1 and column 1:

DiagMatrix "KKT[1][1]" with 0 rows and columns, and with diagonal elements:
DenseVector "KKT[1][1]" with 0 elements:
Homogeneous vector, all elements have value 1.0000000000000000e+00
Component for row 2 and column 0:

GenTMatrix "KKT[2][0]" of dimension 2 by 10 with 6 nonzero elements:
KKT[2][0][ 1, 7]= 1.0000000000000000e+00 (0)
KKT[2][0][ 1410, 1]= 1.0000000000000000e+00 (1)
KKT[2][0][ 1, 8]=-1.0000000000000000e+00 (2)
KKT[2][0][ 1, 1]= 1.0000000000000000e+00 (3)
KKT[2][0][ 1, 9]=-1.0000000000000000e+00 (4)
KKT[2][0][ 2, 2]= 1.0000000000000000e+00 (5)
Component for row 2 and column 1:
This component has not been set.
Component for row 2 and column 2:

DiagMatrix "KKT[2][2]" with 2 rows and columns, and with diagonal elements:
DenseVector "KKT[2][2]" with 2 elements:
Homogeneous vector, all elements have value -0.0000000000000000e+00
Component for row 3 and column 0:

GenTMatrix "KKT[3][0]" of dimension 0 by 10 with 3 nonzero elements:
Uninitialized!
Component for row 3 and column 1:

IdentityMatrix "KKT[3][1]" with 0 rows and columns and the factor -1.0000000000000000e+00.
Component for row 3 and column 2:
This component has not been set.
Component for row 3 and column 3:

DiagMatrix "KKT[3][3]" with 0 rows and columns, and with diagonal elements:
DenseVector "KKT[3][3]" with 0 elements:
Homogeneous vector, all elements have value -0.0000000000000000e+00
******* KKT SYSTEM *******
(0) KKT[1][1] = 0.000000000000000e+00
(1) KKT[2][2] = 0.000000000000000e+00
(2) KKT[3][3] = 0.000000000000000e+00
(3) KKT[4][4] = 0.000000000000000e+00
(4) KKT[5][5] = 0.000000000000000e+00
(5) KKT[6][6] = 0.000000000000000e+00
(6) KKT[7][7] = 0.000000000000000e+00
(7) KKT[8][8] = 0.000000000000000e+00
(8) KKT[9][9] = 0.000000000000000e+00
(9) KKT[10][10] = 0.000000000000000e+00
(10) KKT[1][1] = 1.000000000000000e+00
(11) KKT[2][2] = 1.000000000000000e+00
(12) KKT[3][3] = 1.000000000000000e+00
(13) KKT[4][4] = 1.000000000000000e+00
(14) KKT[5][5] = 1.000000000000000e+00
(15) KKT[6][6] = 1.000000000000000e+00
(16) KKT[7][7] = 1.000000000000000e+00
(17) KKT[8][8] = 1.000000000000000e+00
(18) KKT[9][9] = 1.000000000000000e+00
(19) KKT[10][10] = 1.000000000000000e+00
(20) KKT[11][7] = 1.000000000000000e+00
(21) KKT[1420][1] = 1.000000000000000e+00
(22) KKT[11][8] = -1.000000000000000e+00
(23) KKT[11][1] = 1.000000000000000e+00
(24) KKT[11][9] = -1.000000000000000e+00
(25) KKT[12][2] = 1.000000000000000e+00
(26) KKT[11][11] = -0.000000000000000e+00
(27) KKT[12][12] = -0.000000000000000e+00
(28) KKT[302][1] = 0.000000000000000e+00
(29) KKT[13][1] = 3.952525166729972e-323
(30) KKT[13][2] = 1.909796211918152e-313
In Ma27TSolverInterface::InitializeStructure: Using overestimation factor LiwFact = 2.000000e+00

Matrix structure given to MA27 with dimension 12 and 31 nonzero entries:
A[ 1, 1]
A[ 2, 2]
A[ 3, 3]
A[ 4, 4]
A[ 5, 5]
A[ 6, 6]
A[ 7, 7]
A[ 8, 8]
A[ 9, 9]
A[ 10, 10]
A[ 1, 1]
A[ 2, 2]
A[ 3, 3]
A[ 4, 4]
A[ 5, 5]
A[ 6, 6]
A[ 7, 7]
A[ 8, 8]
A[ 9, 9]
A[ 10, 10]
A[ 11, 7]
A[ 1420, 1]
A[ 11, 8]
A[ 11, 1]
A[ 11, 9]
A[ 12, 2]
A[ 11, 11]
A[ 12, 12]
A[ 302, 1]
A[ 13, 1]
A[ 13, 2]
**Return values from MA27AD: IFLAG = 1, IERROR = 4
*** Error from MA27AD * IFLAG = 1 IERROR = 4
The index of a matrix is out of range.
Please check your implementation of the Jacobian and Hessian matrices.
Factorization failed with retval = 4

LowRankAugSystemSolver: AugSystemSolver returned retval = 4 for right hand side.
Total number of variables............................: 10
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 2
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0

[Callback:E] eval_g
[ 0. 0.]
[Callback:R] eval_g
[Callback:E] eval_f
[Callback:R] eval_f

...