How to solve quadratically constrained linear program in R with clarabel

Question

I'm looking for some help converting a quadratic constraint into the form needed for the clarabel optimization package (second-order cone programming solver).

Firstly, here is just the linear part of the problem:

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)

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

And here is test data for the quadratic constraint:

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')

The quadratic constraint I'm trying to impose is sqrt((x-y) %*% Q %*% (x-y)) <= 0.1, which is equivalent to x %*% Q %*% x - 2 * x %*% Q %*% y - (0.1^2 - y %*% Q %*% y) <= 0.

By putting the quadratic term in the objective, I did a simple search to find the approximate solution to have a guide of what to expect:

lambda = 13.8703
x = clarabel(A, b, -a + lambda * -2 * drop(Q %*% y), 2 * lambda * Q, cones=list(z=1, l=20))$x

So far so good, but my attempts at following various guides on convert QCQP to SOCP (eg here, or here) either produce infeasible errors or seem to ignore the constraint entirely.

This is what I've tried (both 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]

Any help would be much appreciated. Thanks!

Answer 1

It turns out my original attempt above was both incorrect in how it passed the cone parameters to clarabel, and also incorrect in that it requires a symmetric matrix square root (as per many of the QCQP to SOCP transformations). The below actually does successfully enforce the quadratic constraint:

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

Unfortunately in my real problem, I've got a very large sparse matrix, so I was looking for a transformation that just uses the Cholesky factorization. Using the above doesn't work due to the linear term, but that can be avoided by adding additional variables with equality constraints for the linear term (note change in A matrix and size of zero cone):

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 was a bit unsatisfactory too though, as doubling the number of variables (and constraints) in my problem seems undesirable from a performance perspective. In the end I found a transformation here on page 9 that can work with a Cholesky factorization:

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