---
title: "Clarabel Solver Examples"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Clarabel Solver Examples}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup, echo = FALSE}
library(clarabel)
```

## Introduction

The first two examples are from the original [Clarabel
documentation](https://clarabel.org) and the third is from
[SCS](https://www.cvxgrp.org/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

```{r}
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]))
cat(sprintf("Solution: (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
```

## 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

```{r}
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]))
cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
```

## 3. Semidefinite Cone Programming

Semidefinite cones have to be specified in a particular form. We
borrow from the [documentation for the SCS
solver](https://www.cvxgrp.org/scs/api/cones.html#sdcone) 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](https://www.cvxgrp.org/scs/index.html) 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}$. 

```{r}
#' 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.

```{r, echo = TRUE}
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, cones = cones)
cat(sprintf("Solution status, description: = (%d, %s)\n",
            s$status, solver_status_descriptions()[s$status]))
cat(sprintf("Solution (x1, x2, x3) = (%f, %f, %f)\n", s$x[1], s$x[2], s$x[3]))
```

## 4. Cone Specifications

The following cones can be specified in Clarabel. 

```{r, echo = FALSE}
parameter_df <- data.frame(
  Parameter = c("z", "l", "q", "s", "ep", "p", "gp"),
  Type = c("integer", "integer", "integer", "integer", "integer", "numeric", "list"),
  Length = c("1", "1", ">= 1", ">= 1", "1", ">= 1", ">= 1"),
  Description = c(
    "Number of primal zero cones (dual free cones), which corresponds to the primal equality constraints",
    "Number of linear cones (non-negative cones)",
    "Vector of second-order cone sizes",
    "Vector of positive semidefinite cone sizes",
    "Number of primal exponential cones",
    "Vector of primal power cone parameters",
    "List of named lists of two items, `a` : the numeric vector of at least 2 exponent terms, and `n` : an integer dimension of generalized power cone parameters"
  ),
  Definition = c("$\\{ 0 \\}^{z}$",
                 "$\\{ x \\in \\mathbb{R}^{l} : x_i \\ge 0, \\forall i=1,\\dots,l \\}$",
                 "$\\{ (t,x) \\in \\mathbb{R}^{q}  :  \\lVert x\\rVert_2  \\leq t \\}$",
                 "Upper triangular part of the positive semidefinite cone $S^s_+$. The elements $x$ of this cone represent the columnwise stacking of the upper triangular part of a positive semidefinite matrix $X \\in S^s_+$, so that $x \\in R^d$ with $d = s(s+1)/2$",
  "$\\{(x, y, z) : y > 0,~~ ye^{x/y} \\le z \\}$",
  "$\\{(x, y, z) : x^p y^{(1-p)} \\ge  \\lVert z\\rVert,~ (x,y) \\ge 0 \\}$ with $p \\in (0,1)$",
  "$\\{(x, y) \\in R^{len(a)} \\times R^n : \\prod\\limits_{a_i \\in a} x_i^{a_i} \\ge \\lVert y\\rVert_2,~ x \\ge 0 \\}$ with $a_i \\in (0,1)$ and $\\sum a_i = 1$"
  )
)
names(parameter_df)[5] <- "Definition (per parameter element)"
knitr::kable(parameter_df)
```
Generalized power cone parameters are specified as list of two-item
lists, with component named $a$ denoting the exponents and the named
component $n$ denoting the dimension. 

One can specify cones in any order if `strict_cone_order` is set to
`FALSE` in the call to `clarabel()` but one has to ensure that
parameter types are strictly specified for the values, e.g. `5L` for integers, `0.`
for reals etc.  

## 5. 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.

```{r}
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]))
cat(sprintf("Solution (x1, x2) = (%f, %f)\n", s$x[1], s$x[2]))
```

Note the different status, which should always be checked in code.