The first two examples are from the original Clarabel documentation and the third is from SCS.
Suppose that we want to solve the following 2-dimensional quadratic programming problem:
\[ \begin{array}{ll} \text{minimize} & 3x_1^2 + 2x_2^2 - x_1 - 4x_2\\ \text{subject to} & -1 \leq x \leq 1, ~ x_1 = 2x_2 \end{array} \]
We will show how to solve this problem using Clarabel in R.
The first step is to put the problem data into the standard form expected by the solver.
The Clarabel solver’s default configuration expects problem data in
the form \(\frac{1}{2}x^\top P x + q^\top
x\).
We therefore define the objective function data as
\[ P = 2 \cdot \begin{bmatrix} 3 & 0 \\ 0 & 2\end{bmatrix} \mbox{ and } q = \begin{bmatrix} -1 \\ -4\end{bmatrix}. \]
The solver’s default configuration expects constraints in the form \(Ax + s = b\), where \(s \in \mathcal{K}\) for some composite cone \(\mathcal{K}\). We have 1 equality constraint and 4 inequalities, so we require the first element of \(s\) to be zero (i.e. the first constraint will correspond to the equality) and all other elements \(s_i \ge 0\). Our cone constraint on \(s\) is therefore
\[ s \in \mathcal K = \{0\}^1 \times \mathbb{R}^4_{\ge 0}. \]
Define the constraint data as
\[ A = \begin{bmatrix} 1 & -2 \\ 1 & 0 \\ 0 & 1 \\ -1 & 0 \\ 0 & -1\end{bmatrix} \mbox{ and } b=\begin{bmatrix} 0 \\ 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}. \]
Note that Clarabel expects inputs in Compressed Sparse Column (CSC) format for both \(P\) and \(A\) and will try to convert them if not so.
P <- Matrix::Matrix(2 * c(3, 0, 0, 2), nrow = 2, ncol = 2, sparse = TRUE)
P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix
q <- c(-1, -4)
A <- Matrix::Matrix(c(1, 1, 0, -1, 0, -2, 0, 1, 0, -1), ncol = 2, sparse = TRUE)
b <- c(0, 1, 1, 1, 1)
cones <- list(z = 1L, l = 4L) ## 1 equality and 4 inequalities, in order
s <- clarabel(A = A, b = b, q = q, P = P, cones = cones)
cat(sprintf("Solution status, description: = (%d, %s)\n",
s$status, solver_status_descriptions()[s$status]))
#> Solution status, description: = (2, Solver terminated with a solution.)
cat(sprintf("Solution: (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
#> Solution: (x1, x2) = (0.428571, 0.214286)We want to solve the following 2-dimensional optimization problem:
\[ \begin{array}{ll} \text{minimize} & x_2^2\\[2ex] \text{subject to} & \left\|\begin{pmatrix} 2x_1 \\ x_2 \end{pmatrix} - \begin{pmatrix} 2 \\ 2 \end{pmatrix}\right\|_2 \le 1 \end{array} \]
The Clarabel solver’s default configuration expects problem data in
the form \(\frac{1}{2}x^\top P x + q^\top
x\).
We therefore define the objective function data as
\[ P = 2 \cdot \begin{bmatrix} 0 & 0 \\ 0 & 1\end{bmatrix} \mbox{ and } q = \begin{bmatrix} 0 \\ 0\end{bmatrix}. \]
The solver’s default configuration expects constraints in the form \(Ax + s = b\), where \(s \in \mathcal{K}\) for some composite cone \(\mathcal{K}\). We have a single constraint on the 2-norm of a vector, so we rewrite
\[ \left\|\begin{pmatrix} 2x_1 \\ x_2 \end{pmatrix} - \begin{pmatrix} 2 \\ 2 \end{pmatrix}\right\|_2 \le 1 \quad \Longleftrightarrow \quad \begin{pmatrix} 1 \\ 2x_1 - 2\\ x_2 - 2 \end{pmatrix} \in \mathcal{K}_{SOC} \] which puts our constraint in the form \(b - Ax \in \mathcal{K}_{SOC}\).
P <- Matrix::Matrix(2 * c(0, 0, 0, 1), nrow = 2, ncol = 2, sparse = TRUE)
P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix
q <- c(0, 0)
A <- Matrix::Matrix(c(0, -2.0, 0, 0, 0, 1.0), nrow = 3, ncol = 2, sparse = TRUE)
b <- c(1, -2, -2)
cones <- list(q = 3L)
s <- clarabel(A = A, b = b, q = q, P = P, cones = cones)
cat(sprintf("Solution status, description: = (%d, %s)\n",
s$status, solver_status_descriptions()[s$status]))
#> Solution status, description: = (2, Solver terminated with a solution.)
cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
#> Solution (x1, x2) = (1.000000, -1.000000)Semidefinite cones have to be specified in a particular form. We borrow from the documentation for the SCS solver which has similar calling conventions.
The symmetric positive semidefinite cone of matrices is the set
\[ \{S \in \mathbf{R}^{k \times k} \mid S = S^\top, x^\top S x \geq 0 \ \forall x \in \mathbb{R}^k \} \]
and for short, we use \(S \succeq 0\) to denote membership. Clarabel vectorizes this cone in a special way which we detail here.
Clarabel assumes that the input data corresponding to semidefinite cones have been vectorized by scaling the off-diagonal entries by \(\sqrt{2}\) and stacking the upper triangular elements column-wise. (SCS uses the lower triangular elements.) For a \(k \times k\) matrix variable (or data matrix) this operation would create a vector of length \(k(k+1)/2\). Scaling by \(\sqrt{2}\) is required to preserve the inner-product.
This must be done for the rows of both \(A\) and \(b\) that correspond to semidefinite cones and must be done independently for each semidefinite cone.
More explicitly, we want to express \(\text{Trace}(Y S)\) as \(\text{vec}(Y)^\top \text{vec}(S)\), where the \(\text{vec}\) operation takes the (assumed to be symmetric) \(k \times k\) matrix
\[ \begin{aligned} S = \begin{bmatrix} S_{11} & S_{12} & \ldots & S_{1k} \\ S_{21} & S_{22} & \ldots & S_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ S_{k1} & S_{k2} & \ldots & S_{kk} \\ \end{bmatrix} \end{aligned} \]
and produces a vector consisting of the upper triangular elements scaled and arranged as
\[ \text{vec}(S) = (S_{11}, \sqrt{2} S_{12}, S_{22}, \sqrt{2}S_{13}, \ldots, \sqrt{2} S_{1k}, \sqrt{2}S_{2k}, \sqrt{2}S_{3k}, \dots, \sqrt{2}S_{k-1,k}, S_{kk}) \in \mathbb{R}^{k(k+1)/2}. \]
To recover the matrix solution this operation must be inverted on the components of the vectors returned by Clarabel corresponding to each semidefinite cone. That is, the off-diagonal entries must be scaled by \(1/\sqrt{2}\) and the upper triangular entries are filled in by copying the values of lower triangular entries. Explicitly, the inverse operation takes vector \(s \in \mathbb{R}^{k(k+1)/2}\) and produces the matrix
\[ \begin{aligned} \text{mat}(s) = \begin{bmatrix} s_{1} & s_{2} / \sqrt{2} & \ldots & s_{k(k-1)/2-1} / \sqrt{2} \\ s_{2} / \sqrt{2} & s_{3} & \ldots & s_{k(k-1)/2} / \sqrt{2} \\ \vdots & \vdots & \ddots & \vdots \\ s_{k(k-1)/2-1} / \sqrt{2} & s_{k(k-1)/2} / \sqrt{2} & \ldots & s_{k(k+1) / 2} \\ \end{bmatrix} \in \mathbb{R}^{k \times k}. \end{aligned} \]
So the cone definition that Clarabel uses is
\[ \mathcal{S}_+^k = \{ \text{vec}(S) \mid S \succeq 0\} = \{s \in \mathbb{R}^{k(k+1)/2} \mid \text{mat}(s) \succeq 0 \}. \]
Below are two functions to implement both \(\text{vec}\) and \(\text{mat}\).
#' Return an vectorization of symmetric matrix using the upper triangular part,
#' still in column order.
#' @param S a symmetric matrix
#' @return vector of values
vec <- function(S) {
n <- nrow(S)
sqrt2 <- sqrt(2.0)
upper_tri <- upper.tri(S, diag = FALSE)
S[upper_tri] <- S[upper_tri] * sqrt2
S[upper.tri(S, diag = TRUE)]
}
#' Return the symmetric matrix from the [vec] vectorization
#' @param v a vector
#' @return a symmetric matrix
mat <- function(v) {
n <- (sqrt(8 * length(v) + 1) - 1) / 2
sqrt2 <- sqrt(2.0)
S <- matrix(0, n, n)
upper_tri <- upper.tri(S, diag = TRUE)
S[upper_tri] <- v / sqrt2
S <- S + t(S)
diag(S) <- diag(S) / sqrt(2)
S
}Consider the problem:
\[ \begin{array}{ll} \text{minimize} & x_1 - x_2 + x_3\\ \text{subject to} & B_1 - A_{11}x_1 - A_{12}x_2 - A_{13}x_3 \succeq 0\\ & B_2 - A_{21}x_1 - A_{22}x_2 - A_{23}x_3 \succeq 0 \end{array} \] where \[ A_{11}=\begin{bmatrix} -7 & -11 \\ -11 & 3\end{bmatrix}\mbox{, } A_{12}=\begin{bmatrix} 7 & -18 \\ -18 & 8\end{bmatrix}\mbox{, } A_{13}=\begin{bmatrix} -2 & -8 \\ -8 & 1\end{bmatrix}, \]
\[ A_{21}=\begin{bmatrix} -21 & -11 & 0\\ -11 & 10 & 8\\ 0 & 8 & 5\end{bmatrix}\mbox{, } A_{22}=\begin{bmatrix} 0 & 10 & 16\\ 10 & -10 & -10\\ 16 & -10 & 3\end{bmatrix}\mbox{, } A_{23}=\begin{bmatrix} -5 & 2 & -17\\ 2 & -6 & 8\\ -17 & 8 & 6\end{bmatrix}, \]
and
\[ B_1=\begin{bmatrix} 33 & -9 \\ -9 & 26\end{bmatrix}\mbox{, } B_2=\begin{bmatrix} 14 & 9 & 40\\ 9 & 91 & 10\\ 40 & 10 & 15\end{bmatrix}. \]
The constraints involve symmetric positive semidefinite cones over variables \(x \in \mathbb{R}^n\) and \(S \in \mathbb{R}^{k \times k}\)
\[ B - \sum_{i=1}^n \mathcal{A}_i x_i = S \succeq 0 \]
where data \(B, \mathcal{A}_1, \ldots, \mathcal{A}_n \in \mathbb{R}^{k \times k}\) are symmetric. We can write this in the canonical form over a new variable \(s \in \mathcal{S}_+^k\):
\[ \begin{aligned} \begin{align} s &= \text{vec}(S)\\ &= \text{vec}(B - \sum_{i=1}^n \mathcal{A}_i x_i) \\ &= \text{vec}(B) - \sum_{i=1}^n \text{vec}(\mathcal{A}_i) x_i \\ &= b - Ax \end{align} \end{aligned} \]
using the fact that \(\text{vec}\) is linear, where \(b = \text{vec}(B)\) and
\[ A = \begin{bmatrix} \text{vec}(\mathcal{A}_1) & \text{vec}(\mathcal{A}_2) & \cdots & \text{vec}(\mathcal{A}_n) \end{bmatrix} \]
i.e., the vectors \(\text{vec}(\mathcal{A}_i)\) stacked columnwise. This is in a form that we can input into Clarabel. To recover the matrix solution from the optimal solution returned by Clarabel, we simply use \(S^\star = \text{mat}(s^\star)\).
We have two such constraints and will therefore construct the vectors for both cones in order and specify the appropriate dimensions, 2 and 3, respectively.
q <- c(1, -1, 1) # objective: x_1 - x2 + x_3
A11 <- matrix(c(-7, -11, -11, 3), nrow = 2)
A12 <- matrix(c(7, -18, -18, 8), nrow = 2)
A13 <- matrix(c(-2, -8, -8, 1), nrow = 2)
A21 <- matrix(c(-21, -11, 0, -11, 10, 8, 0, 8, 5), nrow = 3)
A22 <- matrix(c(0, 10, 16, 10, -10, -10, 16, -10, 3), nrow = 3)
A23 <- matrix(c(-5, 2, -17, 2, -6, 8, -17, 8, 6), nrow = 3)
B1 <- matrix(c(33, -9, -9, 26), nrow = 2)
B2 <- matrix(c(14, 9, 40, 9, 91, 10, 40, 10, 15), nrow = 3)
A <- rbind(
cbind(vec(A11), vec(A12), vec(A13)), # first psd constraint
cbind(vec(A21), vec(A22), vec(A23)) # second psd constraint
)
b <- c(vec(B1), vec(B2)) # stack both psd constraints
cones <- list(s = c(2, 3)) # cone dimensions
s <- clarabel(A = A, b = b, q = q, cone = cones)
cat(sprintf("Solution status, description: = (%d, %s)\n",
s$status, solver_status_descriptions()[s$status]))
#> Solution status, description: = (2, Solver terminated with a solution.)
cat(sprintf("Solution (x1, x2, x3) = (%f, %f, %f)\n", s$x[1], s$x[2], s$x[3]))
#> Solution (x1, x2, x3) = (-0.367752, 1.898333, -0.887460)Clarabel has a number of parameters that control its behavior,
including verbosity, time limits, and tolerances; see help on
clarabel_control(). As an example, in the last problem, we
can reduce the number of iterations.
P <- Matrix::Matrix(2 * c(0, 0, 0, 1), nrow = 2, ncol = 2, sparse = TRUE)
P <- as(P, "symmetricMatrix") # P needs to be a symmetric matrix
q <- c(0, 0)
A <- Matrix::Matrix(c(0, -2.0, 0, 0, 0, 1.0), nrow = 3, ncol = 2, sparse = TRUE)
b <- c(1, -2, -2)
cones <- list(q = 3L)
s <- clarabel(A = A, b = b, q = q, P = P, cones = cones,
control = list(max_iter = 3)) ## Reduced number of iterations
cat(sprintf("Solution status, description: = (%d, %s)\n",
s$status, solver_status_descriptions()[s$status]))
#> Solution status, description: = (5, Solver terminated with a solution (reduced accuracy))
cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
#> Solution (x1, x2) = (1.000000, -0.999998)Note the different status, which should always be checked in code.