# Clarabel Solver Examples

## Introduction

The first two examples are from the original Clarabel documentation and the third is from SCS.

## 1. Basic Quadratic Program Example

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.

### 1.1. Objective function

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}.$

### 1.2. Constraints

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.

### 1.3. Solution

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)

### 2. Basic Second-order Cone Programming Example

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}$

### 2.1. Objective function

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}.$

### 2.2. Constraints

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}$$.

### 2.3. Solution

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)

## 3. Semidefinite Cone Programming Example

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
}

### 3.1. Example

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) ## 4. Control parameters 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.