I'm trying to reproduce the example "robust least-squares" from the cvxopt documentation:
import numpy as np
import cvxopt
def robls(A, b, rho):
m, n = A.size
def F(x=None, z=None):
if x is None: return 0, cvxopt.matrix(0.0, (n,1))
y = A*x-b
w = cvxopt.sqrt(rho + y**2)
f = sum(w)
Df = cvxopt.div(y, w).T * A
assert not np.isnan(f)
assert not np.isnan(np.array(Df)).any()
if z is None: return f, Df
H = A.T * cvxopt.spdiag(z[0]*rho*(w**-3)) * A
assert not np.isnan(np.array(H)).any()
return f, Df, H
return cvxopt.solvers.cp(F)['x']
The example is supposed to solve min sqrt{ρ + (Ax - b)²} for x. If the matrix A is identity, then the solution will obviously be x=b. This also works nicely for
b = np.ones(2)
robls(cvxopt.matrix(np.eye(2)), cvxopt.matrix(b), 1e-2)
and also when I use b = 3 * np.ones(2) or b = 4 * np.ones(2) instead.
However, for b = 2 * np.ones(2) I get an error:
pcost dcost gap pres dres
0: 0.0000e+00 4.0050e+00 1e+00 1e+00 1e+00
1: -1.5820e+03 1.5840e+03 1e-02 6e+02 1e+00
2: -9.8358e+10 9.8358e+10 1e-04 4e+10 1e+00
3: -2.3551e+34 2.3551e+34 1e-06 9e+33 1e+00
4: -3.2329e+104 3.2329e+104 1e-08 1e+104 1e+00
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-17-5e542f22beda> in <module>
16
17 b = 2 * np.ones(2)
---> 18 robls(cvxopt.matrix(np.eye(2)), cvxopt.matrix(b), 1e-2)
<ipython-input-17-5e542f22beda> in robls(A, b, rho)
13 assert not np.isnan(np.array(H)).any()
14 return f, Df, H
---> 15 return cvxopt.solvers.cp(F)['x']
16
17 b = 2 * np.ones(2)
~/.anaconda3/envs/myenv2020/lib/python3.7/site-packages/cvxopt/cvxprog.py in cp(F, G, h, dims, A, b, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal, **kwargs)
1964
1965 sol = cpl(c, F_e, G_e, h, dims, A_e, b, kktsolver_e, xnewcopy_e,
-> 1966 xdot_e, xaxpy_e, xscal_e, options = options)
1967
1968 sol['x'] = sol['x'][0]
~/.anaconda3/envs/myenv2020/lib/python3.7/site-packages/cvxopt/cvxprog.py in cpl(c, F, G, h, dims, A, b, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal, **kwargs)
1054 while backtrack:
1055 xcopy(x, newx); xaxpy(dx, newx, alpha = step)
-> 1056 t = F(newx)
1057 if t is None: newf = None
1058 else: newf, newDf = t[0], t[1]
~/.anaconda3/envs/myenv2020/lib/python3.7/site-packages/cvxopt/cvxprog.py in F_e(x, z)
1778 else:
1779 if z is None:
-> 1780 v = F(x[0])
1781 if v is None or v[0] is None: return None, None
1782 val = matrix(v[0], tc = 'd')
<ipython-input-7-b500673e8f31> in F(x, z)
4 if x is None: return 0, cvxopt.matrix(0.0, (n,1))
5 y = A*x-b
----> 6 w = cvxopt.sqrt(rho + y**2)
7 f = sum(w)
8 Df = cvxopt.div(y, w).T * A
ValueError: domain error
By debugging into the code, I see that the vector x is filled with nan:
ipdb> p np.array(x)
array([[nan],
[nan]])
Why is that? Also, how do I fix it?
