Hi, I’m trying to solve the following problem, which is the formulation of the multivariate Denton method, if anyone knows about it. I’m using the R package CVXR.
min_y (y-x)‘Ω(y-x)
s.t. Hy = (vecz’ vecY’)', y >= 0
When the solver is finished I get a solution y which does not satisfy Hy = (vecz’ vecY’)'. The discrepancies are as big as -0.14.
I can’t seem to figure out why this is happening. I’ve tried all the free solvers and also I’ve tried setting a lower tolerance, but to no avail.
This is my code, given that matrices Omega and H are built.
y <- Variable(rows = nrow(x), cols = ncol(x))
vecy <- vec(y)
problem <- Problem(objective = Minimize(quad_form(vecy, Omega)),
constraints = list(H%*%vecy == c(vecz, vecY),
vecy >= 0))
result <- solve(problem, solver = solver, feastol = 1e-5)
To check the output I use the very same formula for the constraint:
round(H%*%stack(as.data.frame(solution))$values - c(vecz, vecY),5)
and it can be seen that there are multiple entries which don’t satisfy the constraints. The solution status nontheless is “optimal”.
[,1]
[1,] 0.00000
[2,] 0.00000
[3,] 0.00000
[4,] 0.00000
[5,] 0.00000
[6,] 0.00000
[7,] 0.00000
[8,] 0.00000
[9,] 0.00000
[10,] 0.00000
[11,] 0.00000
[12,] 0.00000
[13,] 0.00000
[14,] 0.00000
[15,] 0.00000
[16,] 0.00000
[17,] 0.00000
[18,] 0.00000
[19,] 0.00000
[20,] 0.00000
[21,] 0.00000
[22,] 0.00000
[23,] 0.00000
[24,] 0.00000
[25,] 0.00000
[26,] 0.00000
[27,] 0.00000
[28,] 0.00000
[29,] 0.00000
[30,] 0.00000
[31,] 0.00000
[32,] 0.00000
[33,] 0.00000
[34,] 0.00000
[35,] 0.00000
[36,] 0.00000
[37,] 0.00000
[38,] 0.00000
[39,] 0.00000
[40,] 0.00000
[41,] 0.00000
[42,] 0.00000
[43,] 0.00000
[44,] 0.00000
[45,] -0.14286
[46,] -0.14286
[47,] -0.14286
[48,] -0.14285
[49,] 0.00001
[50,] 0.00000
[51,] 0.00000
[52,] 0.00000
[53,] 0.00000
[54,] 0.00000
[55,] 0.00000
[56,] 0.00000
[57,] 0.00000
[58,] 0.00000
[59,] 0.00000
[60,] 0.00000
[61,] 0.14286
[62,] 0.00000
[63,] 0.00000
[64,] 0.00000
[65,] 0.00000
[66,] 0.00000
[67,] 0.00000
[68,] 0.00000
[69,] 0.00000
[70,] 0.00000
[71,] 0.00000
[72,] 0.00000
[73,] 0.00001
[74,] 0.14285
[75,] 0.00000
[76,] 0.00000
[77,] 0.00000
[78,] 0.00000
[79,] 0.00000
[80,] 0.00000
[81,] 0.00000
[82,] 0.00000
[83,] 0.00000
[84,] 0.00000
[85,] 0.00000
[86,] 0.00000
[87,] 0.14286
[88,] 0.00000
The solution to this problem has a nice expression in matrix form since it’s quadratic with linear constraints (except for the non-negativity constraint). When computed in that way, it also exhibits this problem for this data.
Edit: I’ve just seen the verbose option. Here’s the console log:
problem: variables n = 294, constraints m = 382
nnz(P) + nnz(A) = 1023
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
rho = 1.00e-01 (adaptive),
sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off
iter objective pri res dua res rho time
1 0.0000e+00 2.41e+05 3.57e+07 1.00e-01 5.98e-04s
200 1.3044e+00 1.43e-01 9.16e-05 1.60e-06 3.83e-03s
400 2.7572e-01 1.43e-01 2.97e-05 1.60e-06 6.75e-03s
600 1.0134e-01 1.43e-01 1.18e-05 1.60e-06 1.09e-02s
650 8.6386e-02 1.43e-01 9.44e-06 1.60e-06 1.18e-02s
status: solved
solution polish: unsuccessful
number of iterations: 650
optimal objective: 0.0864
run time: 1.24e-02s
optimal rho estimate: 1.00e-06