ValueError: domain error in misc.py

[Yihao Lu]

Hi there,
I am currently using v1.1.3 to solve conic programming. When running
with around 20 vars, I see ValueError in misc.py line 449: s[:m] =
base.sqrt( s[:m] ) .
Some people proposed using s[:m] = base.sqrt( abs(s[:m]) ) , as it is
claimed s is only to control number of iterations. Is this solution
correct or is there a better one? thanks

[Lieven Vandenberghe]

s should remain positive, so taking the absolute value will not help.

Can you send us an example?

Lieven

[Yihao Lu]

Hi Lieven,
thanks for reply.
I am not sure I still have the example, but I think it happens when I
set abstol/reltol/feastol to around 10^-14 range. I do see 's' becomes
negative. It works when I change them to 10e-11 range. Do you still
need me find an example?

On Feb 20, 7:34 pm, Lieven Vandenberghe

[Martin]

We will need an example in order to investigate this further.

On Thursday, April 4, 2013 10:02:31 PM UTC+2, Paul Vecchio wrote:

Has there been any update on solving this issues?  I am also running into something similar.  My S continuously grows until it reaches nan. 

On Friday, January 18, 2013 2:55:41 PM UTC-8, Ely Spears wrote:

I am also seeing this error. For me, the entries don't become negative but they can become a NaN value. I see them increasing dramatically as the number of iterations grows, and for a 1000x1000 problem, when the iteration hits around 50, that's when I see the value error. I added the "abs" trick just to be sure, and I do still see the ValueError.

[Paul Vecchio]

My code is based off the following process of fitting Rational functions with a quadratic solver: http://hal.inria.fr/docs/00/67/88/85/PDF/main_tvcg.pdf

The paper uses different notation than CVXOPT so I'll try to make sure that all my notation is correct in this example.  Let's use two examples.  1 where i'm trying to fit a simple constant set of data and another where I'm trying to fit a simple diagonal line.

The core of the problem is solving for the following. 

arg min ||x||_2 (Which is minimizing the euclidean norm of the x coefficients )

subject to

Px  - 1/cn ||P||_2 >= 0, which can be rewritten Px >= 1/cn ||P||_2

To get this  I set:

P  = to an nxn identity matrix so that the general quadratic .5 * x' * Px + qT*x reduces to the euclidean norm
q  = is a vector of zeros equal n
G  = n x 2m my matrix of inequalities (n is the number of coefficients and m is the number of data points * 2)
cn = cond(P)
h  = 1/cn * rowNorm(P)

Since CVXOPT's equality has Gx <= h I have to multiply G and h by -1 to flip the inequality to the form of Gx >= h. 

At this point for the simple constant value data fit with 2 coefficients.

P = [ 1.00e+00  0.00e+00]
      [ 0.00e+00  1.00e+00]

q = [ 0.00e+00]
      [ 0.00e+00]

G = [-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[-1.00e+00  5.00e-01]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]
[ 1.00e+00 -1.50e+00]

h =  [-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-2.62e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]
[-4.23e-01]


At this point CVX returns the correct answer. 'coeffs': array([[ 0.70710637],  [ 0.70710673]]), 'cn': 2000.0005000002723, 'point': (1, 1)}
Basically meaning that the rational function is 0.70710637/70710637 * 1 which matches the input array of 1 at all x values. This essentially works for any number of coefficients for a constant value. 

The second case is when i try to fit anything besides a constant value.  In this case lets assume that we are fitting a diagonal line. from x = 0..1 and y = 0..1

in this case we'll try to fit with 4 coefficients.

P = [ 1.00e+00  0.00e+00  0.00e+00  0.00e+00]
[ 0.00e+00  1.00e+00  0.00e+00  0.00e+00]
[ 0.00e+00  0.00e+00  1.00e+00  0.00e+00]
[ 0.00e+00  0.00e+00  0.00e+00  1.00e+00]


q = [ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]

G = [-1.00e+00 -0.00e+00 -0.00e+00  0.00e+00]
[-1.00e+00 -1.11e-01 -1.23e-02  1.11e-01]
[-1.00e+00 -2.22e-01 -4.94e-02  2.22e-01]
[-1.00e+00 -3.33e-01 -1.11e-01  3.32e-01]
[-1.00e+00 -4.44e-01 -1.98e-01  4.43e-01]
[-1.00e+00 -5.56e-01 -3.09e-01  5.54e-01]
[-1.00e+00 -6.67e-01 -4.44e-01  6.65e-01]
[-1.00e+00 -7.78e-01 -6.05e-01  7.75e-01]
[-1.00e+00 -8.89e-01 -7.90e-01  8.86e-01]
[-1.00e+00 -1.00e+00 -1.00e+00  9.97e-01]
[ 1.00e+00  0.00e+00  0.00e+00 -0.00e+00]
[ 1.00e+00  1.11e-01  1.23e-02 -1.11e-01]
[ 1.00e+00  2.22e-01  4.94e-02 -2.23e-01]
[ 1.00e+00  3.33e-01  1.11e-01 -3.34e-01]
[ 1.00e+00  4.44e-01  1.98e-01 -4.46e-01]
[ 1.00e+00  5.56e-01  3.09e-01 -5.57e-01]
[ 1.00e+00  6.67e-01  4.44e-01 -6.69e-01]
[ 1.00e+00  7.78e-01  6.05e-01 -7.80e-01]
[ 1.00e+00  8.89e-01  7.90e-01 -8.92e-01]
[ 1.00e+00  1.00e+00  1.00e+00 -1.00e+00]

h = [-9.50e-04]
[-9.62e-04]
[-9.97e-04]
[-1.06e-03]
[-1.14e-03]
[-1.24e-03]
[-1.37e-03]
[-1.52e-03]
[-1.70e-03]
[-1.90e-03]
[-9.50e-04]
[-9.62e-04]
[-9.97e-04]
[-1.06e-03]
[-1.14e-03]
[-1.24e-03]
[-1.37e-03]
[-1.53e-03]
[-1.70e-03]
[-1.90e-03]

This is the point where i get a domain error where s and z both basically grow to large number until getting to NaN. When i wrote this code in matlab it seemed to work fine, so i'm not exactly sure what i'm doing wrong.  If you need any more information let me know.

[Martin]

I'm not quite sure that I follow your reformulation of the problem. Are you solving the problem as a QP with solvers.qp? Have you tried formulating the problem as a second-order cone program and solving it using solvers.conelp? 

Martin

[Paul Vecchio]

I am using solvers.qp for the solution.  I have not tried formulating the problem as a second-order cone problem though.   I don’t think that is something I really want to do as I am not familiar with that process and I do not know what issues that will cause by deviating from the papers process.  When I wrote this in MatLab I used the quadprog function to solve this problem and did not seem to run into any trouble.

You seem to have formed the problem correctly assuming that b = [b1,b2, bn] includes the delta (inverse of the conditional number of G) * Euclidian norm of G,   although you might have missed the part where I multiply the constraints of the problem to match the form that CVXOPT accepts.

 CVXOPT only accepts the constraints in the form  Gx <= h,  which is not the form that the paper takes.   Instead it is Gx >= h.   

The above constraint is of that form but the signs are flipped.  The inequality is multiplied by -1 to achieve the desired form   -Gx <=  -h which is equal to Gx >= h,  as shown in the CVX tutorial from: http://courses.csail.mit.edu/6.867/wiki/images/a/a7/Qp-cvxopt.pdf

I understand you concern about the rows being shared but doesn’t the multiplication by -1 flip your sign, making the inequalities feasible?

On Saturday, April 6, 2013 3:41:45 AM UTC-7, Martin wrote:

I'm not quite sure that I follow your reformulation of the problem. Are you solving the problem as a QP with solvers.qp? Have you tried formulating the problem as a second-order cone program and solving it using solvers.conelp? I don't know if I understand the problem that you are trying to solve correctly.  It seems that the problem in the paper that you refer to can be expressed as. 

minimize ||c||_2
subject to
B*c1 - Fl*c2 >= b1
-B*c1 + Fu*c2 >= b2

where c =[c1;c2] and B, Fl, and Fu are matrices.  Equivalently, the constraints can be expressed as

b1 + Fl*c2 <= B*c1 <= Fu*c2 - b2

Notice that if Fu and Fl share a row (this seems to be the case in your example), then you have an inequality of the form

pos. const. + r' *c2 <= bi * c1 <= r'*c2 - pos. const.

where r' is the common row. Notice that the lower bound is larger than the upper bound, so that the set of inequalities is infeasible.

Am I missing something?

 

b1 + Fl*c2 <= B*c1 <=

[Martin]

In your second example, you have the constraints

   9.50e-4 <= x1 <= -9.50e-4 

which follow from rows 1 and 11 of G*x <= h, and hence the problem is infeasible. 

Multiplying both sides of an infeasible inequality by -1 does not make it feasible.

[Paul Vecchio]

I understand you objection to the example dataset that I posted.  I went back the problem and code and I’m pretty sure that in its current form the inequalities are correct. I’ll post an example dataset if you want me to later.

I implemented a try, catch statement around the quadratic problem and instead of crashing I would just ignore the error and move to my next iteration.  However, I still receive the error Terminated (singular KKT matrix), which I believe is produced when there is no solution to the set of basis functions used to solve the problem.  Eventually if there is a solution it will converge to a feasible set of coefficients that product a rational function that fits the input data within the specified error bounds.  If you believe this is not a reasonable explanation let me know.  In any case I saw your suggestion as well as read in numerous places that the SOCP solver is more robust than the QP solver and tried to reformulate the problem.   From what I’ve read, LPs, QCQPs, and QPs are special cases of SOCPs.  However, inputting reformulations I’ve see into the conelq or socp solver does not seem straight forward.   I’ve explained in my attached PDF, how I tried to reformulate the QP as an SOCP and then input it into CVXOPT.

[Martin]

Hi Paul,

The problem

minimize   t 

subject to  Ax >= b,  ||x|| <= t

can be expressed as the following "cone LP":

minimize   1*t + 0'*x

subject to

  [  0  -A  ] [ t ]     [ s1 ]    [ -b ]

  [ -1   0  ] [ x ]  +  [ s2 ]  = [  0 ]

  [  0  -I  ]           [ s3 ]    [  0 ]

   s1 >= 0,  ||s3|| <= s2

Now suppose A is m-by-n, and b is m-by-1. Then you can solve this in CVXOPT with the following code:

from cvxopt import matrix, spmatrix, sparse, solvers

A = ...

b = ...

c = matrix(0.0, (n+1,1)); c[0] = 1.0

G = sparse([[spmatrix(-1.0,[m],[0],(m+1+n,1))],[-A, spmatrix([],[],[],(1,n)), -spmatrix(1.0,range(n),range(n))]])

h = matrix([-b, matrix(0.0, (n+1,1))])

sol = solvers.conelp(c, G, h, {'l' : m, 'q' : [n+1], 's' : [] })

x = sol['x'][1:]

Note that this does not address the issue with infeasible constraints.

Martin

[Paul Vecchio]

Thanks Martin.  I'll try it out.

[Paul Vecchio]

So i tried out your code and it seems like there might be an error somewhere.

Let a be a matrix A of size m x n where m is the number of rows and n is the number of columns. 

An A of (200,3) generates the following shapes:
c = (4,1)
G = (204,4)
h = (204,1)

When putting these values into conelp I get the error

cvxopt\coneprog.py", line 522, in conelp raise TypeError("'h' must be a 'd' matrix of size (%d,1)" %cdim) TypeError: 'h' must be a 'd' matrix of size (4,1)

Though it is TypeError it seems that also the dimensions of the h matrix are also incorrect.  I'm going to continue looking into this myself but maybe you might have some insight?

[Martin]

I am not sure what the problem is. You can generate a strictly feasible problem instance as follows:

from cvxopt import matrix, spmatrix, sparse, solvers, normal

m,n = 200,3

A = normal(m,n)

x0 = normal(n,1)

b = A*x0 - 0.01

c = matrix(0.0, (n+1,1)); c[0] = 1.0

G = sparse([[spmatrix(-1.0,[m],[0],(m+1+n,1))],[-A, spmatrix([],[],[],(1,n)), -spmatrix(1.0,range(n),range(n))]])

h = matrix([-b, matrix(0.0, (n+1,1))])

sol = solvers.conelp(c, G, h, {'l' : m, 'q' : [n+1], 's' : [] })

x = sol['x'][1:]

[Paul Vecchio]

Thanks for all the help Martin.  Your simple case helped me determine where my error was.  The conelp solver does not display the "singular KKT" error that the qp solver would get when the problem has no solution. 

[Yury Rojek]

I think I have an example that cvxopt fails to find a solution with exactly the same error. I might be wrong, but I think that zero is a feasible point. Here is my data:
===================================================
c=
[-1.95e+10]
[-2.39e+11]
[-1.09e+10]
[ 1.39e+10]
[ 5.00e-01]
[ 5.00e-01]
[ 0.00e+00]

G0=
[-1.00e+00     0         0         0         0         0         0    ]
[    0     -1.00e+00     0         0         0         0         0    ]
[    0         0     -1.00e+00     0         0         0         0    ]
[    0         0         0     -1.00e+00     0         0         0    ]
[    0         0         0         0     -1.00e+00     0         0    ]
[    0         0         0         0         0     -1.00e+00     0    ]
[    0         0         0         0         0         0     -1.00e+00]
[ 1.00e+00     0         0         0         0         0         0    ]
[    0      1.00e+00     0         0         0         0         0    ]
[    0         0      1.00e+00     0         0         0         0    ]
[    0         0         0      1.00e+00     0         0         0    ]
[-1.95e+10 -2.39e+11 -1.09e+10  1.39e+10     0         0      1.00e+00]

h0=
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 5.78e-03]
[ 3.76e-03]
[ 1.04e-02]
[ 7.23e-03]
[ 0.00e+00]

G_1=
[ 0.00e+00  0.00e+00  0.00e+00  0.00e+00 -5.00e-01  0.00e+00  0.00e+00]
[ 5.00e+04  3.34e+05  3.64e+04 -3.51e+04  0.00e+00  0.00e+00  0.00e+00]
[-4.49e+04  1.45e+05 -1.39e+04 -6.06e+03  0.00e+00  0.00e+00  0.00e+00]
[-8.51e+03 -1.50e+05  1.84e+04  1.87e+03  0.00e+00  0.00e+00  0.00e+00]
[-3.46e+04 -5.45e+04 -6.60e+03  1.39e+04  0.00e+00  0.00e+00  0.00e+00]
[ 3.55e+02  1.70e+05  2.02e+03  5.33e+03  0.00e+00  0.00e+00  0.00e+00]
[-8.42e+03 -3.93e+04  4.66e+03 -1.17e+04  0.00e+00  0.00e+00  0.00e+00]
[ 6.04e+03 -4.51e+04 -2.11e+03 -1.17e+04  0.00e+00  0.00e+00  0.00e+00]
[ 1.13e+03 -4.95e+04  4.95e+03 -2.28e+03  0.00e+00  0.00e+00  0.00e+00]
[ 1.08e+03 -1.24e+05 -1.92e+03  9.56e+03  0.00e+00  0.00e+00  0.00e+00]
[-4.59e+03  3.88e+04 -3.47e+03  7.00e+03  0.00e+00  0.00e+00  0.00e+00]
[-5.03e+03  5.87e+04 -1.35e+03  8.01e+02  0.00e+00  0.00e+00  0.00e+00]
[ 7.05e+02  4.25e+04 -1.07e+03 -3.79e+03  0.00e+00  0.00e+00  0.00e+00]
[ 1.79e+03 -4.11e+04  2.05e+03  1.12e+03  0.00e+00  0.00e+00  0.00e+00]
[ 1.45e+02  4.86e+04 -8.33e+02  6.40e+01  0.00e+00  0.00e+00  0.00e+00]
[-1.16e+03  1.73e+04  3.76e+02 -6.63e+02  0.00e+00  0.00e+00  0.00e+00]
[ 5.08e+03 -7.76e+03  3.07e+03  1.83e+03  0.00e+00  0.00e+00  0.00e+00]
[-3.98e+03  1.45e+04 -2.18e+03  1.54e+02  0.00e+00  0.00e+00  0.00e+00]
[ 0.00e+00  0.00e+00  0.00e+00  0.00e+00 -5.00e-01  0.00e+00  0.00e+00]

h_1=
[ 5.00e-01]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[-5.00e-01]

G_2=
[    0         0         0         0         0     -5.00e-01     0    ]
[-1.70e+04     0         0         0         0         0         0    ]
[    0     -1.08e+05     0         0         0         0         0    ]
[    0         0     -1.24e+04     0         0         0         0    ]
[    0         0         0     -1.47e+04     0         0         0    ]
[    0         0         0         0         0     -5.00e-01     0    ]

h_2=
[ 5.00e-01]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[-5.00e-01]

G_3=
[    0         0         0         0         0         0     -5.00e-01]
[-1.16e+07     0         0         0         0         0         0    ]
[    0     -1.92e+08     0         0         0         0         0    ]
[    0         0     -1.89e+07     0         0         0         0    ]
[    0         0         0     -2.08e+07     0         0         0    ]
[    0         0         0         0         0         0     -5.00e-01]

h_3=
[ 5.00e-01]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[ 0.00e+00]
[-5.00e-01]

A=
[ 1.55e+06  1.52e+07  7.69e+05 -1.75e+07  0.00e+00  0.00e+00  0.00e+00]

b=
[ 0.00e+00]

===================================================

I apologize in advance if it is my bug and the task is indeed infeasible.

Thanks

Yury

[Martin]

The problem that you're trying to solve is very poorly scaled, so you cannot expect much in terms of accuracy when using double precision arithmetic. 

[Yury Rojek]

Martin,

You are right. It is poorly scaled. But does it have to give ValueError on negative "s" in this case?

Thanks, I will re-scale everything and try again.

Yury