Converting QCQP into SOCP for clarabel

[aquasync]

Hi there,

This is probably more my lack of understanding of SOCP, rather than an actual issue with clarabel (which works really well!), though may be a good addition to the docs: I don't seem to be able to get clarabel to enforce my quadratic constraints.

Here is the test code I'm using:

library(clarabel)
library(Matrix)

a = c(1.425, 0.27, -0.085, -0.733, -0.534, 1.402, -0.199, -0.875, -1.586, -2.126)
xl = rep(0, 10)
xu = rep(1, 10)
A = rbind(rep(1, 10), Diagonal(10), -Diagonal(10))
b = c(1, xu, -xl)

# purely lp
x = clarabel(A, b, -a, cones=list(z=1, l=20))$x

y = c(0.242, 0.058, 0.204, 0.112, 0.123, 0.143, 0, 0.065, 0, 0.053)
Q = structure(c(1.291, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.018, 0.525, 0,
                0, 0, 0, 0, 0, 0, 0, 0.034, 0.02, 0.543, 0, 0, 0, 0, 0, 0, 0,
                0.012, 0.011, 0.015, 0.536, 0, 0, 0, 0, 0, 0, 0.039, 0.021, 0.036,
                0.012, 0.874, 0, 0, 0, 0, 0, 0.064, 0.02, 0.046, 0.009, 0.054,
                0.735, 0, 0, 0, 0, 0.02, 0.007, 0.019, 0.001, 0.025, 0.035, 1.791,
                0, 0, 0, 0.031, 0.018, 0.029, 0.011, 0.032, 0.043, 0.03, 1.242,
                0, 0, 0.008, 0.008, 0.012, 0.045, 0.006, 0.004, 0.008, 0.007,
                0.55, 0, 0.005, 0.007, 0.009, 0.041, 0.003, 0, 0.01, 0.006, 0.062,
                0.479), .Dim = c(10L, 10L))
Q = as(as(Q + lower.tri(Q) * t(Q), 'CsparseMatrix'), 'symmetricMatrix')

# hoisting quadratic term into objective and searching for appropriate weight to satisfy constraint:
# sqrt((x-y) %*% Q %*% (x-y)) <= 0.1
lambda = 13.8703
x = clarabel(A, b, -a + lambda * -2 * drop(Q %*% y), 2 * lambda * Q, cones=list(z=1, l=20))$x

By rewriting the quadratic constraint into standard form (x %*% Q %*% x - 2 * x %*% Q %*% y - (0.1^2 - y %*% Q %*% y) <= 0), I can add it as a second-order cone, but my attempts so far are infeasible:

S = chol(Q)
lin = drop(-2 * Q %*% y)
rhs = drop(0.1^2 - y %*% Q %*% y)
bb = -drop(solve(S, lin))
gamma = -sqrt(drop(crossprod(bb)) + rhs)

res = clarabel(rbind(A, 0, S), c(b, rhs, lin), -a, cones=list(z=1, l=20, q=11))
solver_status_descriptions()[res$status]

res = clarabel(rbind(A, 0, S), c(b, gamma, bb), -a, cones=list(z=1, l=20, q=11))
solver_status_descriptions()[res$status]

I take it I'm just missing something basic with the QCLP to SOCP conversion here?

[goulart-paul]

I am not an R user so can't directly suggest a fix for you, but perhaps it is just a minor problem in your reformulation. I think you want to use the hyperbolic constraint transformation that appears in the bottom half of slide 7 in this talk: https://inst.eecs.berkeley.edu/~ee127/fa19/Lectures/12_socp.pdf. I think that's not quite the same as what you have.

[aquasync]

Thanks for the pointer Paul - I was following a different transformation similar to that seen in Matlab here.

In the end it seems my problem was in using the Cholesky factorization, whereas that particular transformation (and the one from your link) requires a symmetric matrix square root.

Eg this works:

S = as(expm::sqrtm(Q), 'sparseMatrix')
lin = -2 * drop(Q %*% y)
bb = 0.5*drop(solve(S, -lin))
target = 0.1^2 - drop(y %*% Q %*% y)
x = clarabel(rbind(A, -0.1, S), c(b, 0, bb), -a, cones=list(z=1, l=20, q=11))$x
sqrt((x - y) %*% Q %*% (x - y)) # => 0.1

I was trying to avoid a matrix square root though, as I have a huge sparse matrix. I had some luck by using adding additional variables with equality constraints to avoid the need for a linear part on the quadratic:

S = chol(Q)
A2 = rbind(cbind(Diagonal(10), -Diagonal(10)), cbind(A, Matrix(0, nrow(A), 10)))
b2 = c(y, b)
a2 = c(a, rep(0, 10))
x = clarabel(rbind(A2, c(rep(-0.1, 20)), cbind(Matrix(0, 10, 10), S)), c(b2, rep(0, 11)), -a2, cones=list(z=11, l=20, q=11))$x[1:10]
sqrt((x - y) %*% Q %*% (x - y)) # => 0.1

That is however not particularly appealing as it means having to almost double the number of variables in the problem.

I finally had some luck following the transformation from page 9 here, which doesn't require the matrix square root:

S = chol(Q)
lin = -2 * drop(Q %*% y)
target = 0.1^2 - drop(y %*% Q %*% y)
x = clarabel(rbind(A, lin/2, S, -lin/2), c(b, (1+target)/2, rep(0, 10), (1-target)/2), -a, cones=list(z=1, l=20, q=12))$x
sqrt((x - y) %*% Q %*% (x - y)) # => 0.1

Thanks again.