# A Gentle Introduction to `CVXR`

## Introduction

Welcome to `CVXR`

: a modeling language for describing and solving convex optimization problems that follows the natural, mathematical notation of convex optimization rather than the requirements of any particular solver. The purpose of this document is both to introduce the reader to `CVXR`

and to generate excitement for its possibilities in the field of statistics.

Convex optimization is a powerful and very general tool. As a practical matter, the set of convex optimization problems includes almost every optimization problem that can be solved exactly and efficiently (i.e. without requiring an exhaustive search). If an optimization problem can be solved, it is probably convex. This family of problems becomes even larger if you include those that can be solved *approximately* and efficiently. To learn more about the mathematics and applications of convex optimization, see Boyd and Vandenberghe 2009.

Convex optimization systems written in other languages are already widely used in practical applications. These include YALMIP and CVX (Matlab), CVXPY (Python), and Convex.jl (Julia). `CVXR`

Shares a lot of its code base with CVXcanon and CVXPY. As far as we know, this is the first full-featured general convex optimization package for R.

One of the great headaches of conventional numerical optimization is the process of deciding which algorithm to use and how to set its parameters. In convex optimization, the particular algorithm matters much less. So while a user of `CVXR`

is still free to choose from a number of different algorithms and to set algorithm parameters as they please, the vast majority of users will not need to do this. `CVXR`

will just work.

The uses for convex optimization in statistics are many and varied. A large number of parameter-fitting methods are convex, including least-squares, ridge, lasso, and isotonic regression, as well as many other kinds of problems such as maximum entropy or minimum Kullback-Leibler divergence over a finite set.

All of these examples, at least in their most basic forms, are established enough that they already have well-crafted R packages devoted to them. If you use `CVXR`

to solve these problems, it will work. It will probably be slower than a custom-built algorithm—for example, glmnet for fitting lasso or ridge regression models—but it will work. However, this is not the true purpose of `CVXR`

. If you want to build a well-established model, you should use one of the well-established packages for doing so. If you want to build your *own* model—one that is a refinement of an existing method, or perhaps even something that has never been tried before—then `CVXR`

is the package to do it. The advantage of `CVXR`

over glmnet and the like comes from its flexibility: A few lines of code can transform a problem from commonplace to state-of-the-art, and can often do the work of an entire package in the process.

This document is meant to familiarize you with the `CVXR`

package. It assumes basic knowledge of convex optimization and statistics as well as proficiency with R. A potential user of `CVXR`

should read this entire document and then refer to the tutorial examples.

Happy optimizing!

## Convex Optimization

A convex optimization problem has the following form:

\[ \begin{array}{ll} \mbox{minimize} & f_0(x)\\ \mbox{subject to} & f_i(x) \leq 0, \quad i=1,\ldots,m\\ & g_i(x) = 0, \quad i=1,\ldots,p \end{array} \]

where \(x\) is the variable, \(f_0\) and \(f_1,...,f_m\) are convex, and \(g_1,...,g_p\) are affine. \(f_0\) is called the objective function, \(f_i \leq 0\) are called the inequality constraints, and \(g_i = 0\) are called the equality constraints.

In `CVXR`

, you will specify convex optimization problems in a more convenient format than the one above.

A convex function is one that is upward curving. A concave function is downward curving. An affine function is flat, and is thus both convex and concave.

A convex optimization problem is one that attempts to minimize a convex function (or maximize a concave function) over a convex set of input points.

You can learn much more about convex optimization via Boyd and Vandenberghe (2004) as well as the CVX101 MOOC.

## ‘Hello World’

We begin with one of the simplest possible problems that presents all three of these features:

\[ \begin{array}{ll} \mbox{minimize} & x^2 + y^2 \\ \mbox{subject to} & x \geq 0, \quad 2x + y = 1 \end{array} \]

with scalar variables \(x\) and \(y\). This is a convex optimization problem with objective \(f_0(x,y) = x^2 + y^2\) and constraint functions \(f_1(x,y) = -x\) and \(g_1(x,y) = 2x - y - 1\).

Note that this problem is simple enough to be solved analytically, so we can confirm that `CVXR`

has produced the correct answer. Here’s how we formulate the problem in `CVXR`

.

```
# Variables minimized over
x <- Variable(1)
y <- Variable(1)
# Problem definition
objective <- Minimize(x^2 + y^2)
constraints <- list(x >= 0, 2*x + y == 1)
prob2.1 <- Problem(objective, constraints)
# Problem solution
solution2.1 <- solve(prob2.1)
solution2.1$status
```

`## [1] "optimal"`

`solution2.1$value`

`## [1] 0.2`

`solution2.1$getValue(x)`

`## [1] 0.3999978`

`solution2.1$getValue(y)`

`## [1] 0.2000044`

`# The world says 'hi' back.`

We now turn to a careful explanation of the code. The first lines create two Variable objects, `x`

and `y`

, both of length 1 (i.e. scalar variables).

```
x <- Variable(1)
y <- Variable(1)
```

`x`

and `y`

represent what we are allowed to adjust in our problem in order to obtain the optimal solution. They don’t have values yet, and they won’t until after we solve the problem. For now, they are just placeholders.

Next, we define the problem objective.

`objective <- Minimize(x^2 + y^2)`

This call to `Minimize()`

does *not* return the minimum value of the expression `x^2 + y^2`

the way a call to the native R function `min()`

would do (after all, `x`

and `y`

don’t have values yet). Instead, `Minimize()`

creates an Objective object, which defines the goal of the optimization we will perform, namely to find values for `x`

and `y`

which produce the smallest possible value of `x^2 + y^2`

.

The next line defines two constraints—an inequality constraint and an equality constraint, respectively.

`constraints <- list(x >= 0, 2*x + y == 1)`

Again, counter to what you might ordinarily expect, the expression `x >= 0`

does not return `TRUE`

or `FALSE`

the way `1.3 >= 0`

would. Instead, the `==`

and `>=`

operators have been overloaded to return Constraint objects, which will be used by the solver to enforce the problem’s constraints. (Without them, the solution to our problem would simply be \(x = y = 0\).)

Next, we define our Problem object, which takes our Objective object and our two Constraint objects as inputs.

`prob2.1 <- Problem(objective, constraints)`

Problem objects are very flexible in that they can have 0 or more constraints, and their objective can be to `Minimize()`

a convex expression (as shown above) *or* to `Maximize()`

a concave expression.

The call to `Problem()`

still does not actually *solve* our optimization problem. That only happens with the call to `solve()`

.

`solution2.1 <- solve(prob2.1)`

Behind the scenes, this call translates the problem into a format that a convex solver can understand, feeds the problem to the solver, and then returns the results in a list. For this problem, the list will contain among other things the optimal value of the objective function `x^2 + y^2`

, values for `x`

and `y`

that achieve that optimal objective value (along with a function `solution2.1$getValue`

to retrieve them conveniently), and some accompanying metadata such as `solution2.1$status`

, which confirms that the solution was indeed `"optimal"`

.

`solution2.1$status`

`## [1] "optimal"`

`solution2.1$value`

`## [1] 0.2`

`solution2.1$getValue(x)`

`## [1] 0.3999978`

`solution2.1$getValue(y)`

`## [1] 0.2000044`

In general, when you apply the `solve()`

method to a Problem, several things can happen:

`solution$status == "optimal"`

: The problem was solved. Values for the optimization variables were found, which satisfy all of the constraints and minimize/maximize the objective.`solution$status == "infeasible"`

: The problem was*not*solved because no combination of input variables exists that can satisfy all of the constraints. For a trivial example of when this might happen, consider a problem with optimization variable`x`

, and constraints`x >= 1`

and`x <= 0`

. Obviously, no value of`x`

exists that can satisfy both constraints. In this case,`solution$value`

is`+Inf`

for a minimization problem and`-Inf`

for a maximization problem, indicating infinite dissatisfaction with the result. The values for the input variables are`NA`

.`solution$status == "unbounded"`

: The problem was*not*solved because the objective can be made arbitrarily small for a minimization problem or arbitrarily large for a maximization problem. Hence there is no optimal solution because for any given solution, it is always possible to find something even more optimal. In this case,`solution$opt.val`

is`-Inf`

for a minimization problem and`+Inf`

for a maximization problem, indicating infinite satisfaction with the result. Again, the values of the the input variables are`NA`

.

### Modifying a CVXR Problem

Like any normal R object, the `Problem`

, `Minimize`

, `Maximize`

, and `Constraint`

objects can all be modified and computed upon after creation. Here is an example where we modify the problem we created above by changing its objective and adding a constraint, print the modified problem, check whether it is still convex, and then solve the modified problem:

```
# Modify the problem from example 1
prob2.2 <- prob2.1
objective(prob2.2) <- Minimize(x^2 + y^2 + abs(x-y))
constraints(prob2.2) <- c(constraints(prob2.1), y <= 1)
# Solve the modified problem
solution2.2 <- solve(prob2.2)
# Examine the solution
solution2.2$status
```

`## [1] "optimal"`

`solution2.2$value`

`## [1] 0.2222222`

`solution2.2$getValue(x)`

`## [1] 0.3333333`

`solution2.2$getValue(y)`

`## [1] 0.3333333`

### An Invalid Problem

Unfortunately, you can’t just type any arbitrary problem you like into `CVXR`

. There are restrictions on what kinds of problems can be handled. For example, if we tried to `Maximize()’ the objective from example 2.1, we get an error:

```
prob2.3 <- prob2.1
objective(prob2.3) <- Maximize(x^2 + y^2)
solve(prob2.3)
```

`## Error in CVXR::psolve(a, b, ...): Problem does not follow DCP rules.`

We would get a similar error if we tried to add the constraint `p_norm(x,2) == 1`

. This is because `CVXR`

uses a strict set of rules called Disciplined Convex Programming (DCP) to evaluate the convexity of any given problem. If you follow these rules, you are guaranteed that your problem is convex. If you don’t follow these rules, `CVXR`

will throw an exception. See the last section for further information on DCP.

## Simple Examples

We begin by showing what a standard linear regression problem looks like in `CVXR`

:

### Ordinary Least Squares

For illustration, we generate some synthetic data for use in this example.

```
set.seed(1)
s <- 1
n <- 10
m <- 300
mu <- rep(0, 9)
Sigma <- cbind(c(1.6484, -0.2096, -0.0771, -0.4088, 0.0678, -0.6337, 0.9720, -1.2158, -1.3219),
c(-0.2096, 1.9274, 0.7059, 1.3051, 0.4479, 0.7384, -0.6342, 1.4291, -0.4723),
c(-0.0771, 0.7059, 2.5503, 0.9047, 0.9280, 0.0566, -2.5292, 0.4776, -0.4552),
c(-0.4088, 1.3051, 0.9047, 2.7638, 0.7607, 1.2465, -1.8116, 2.0076, -0.3377),
c(0.0678, 0.4479, 0.9280, 0.7607, 3.8453, -0.2098, -2.0078, -0.1715, -0.3952),
c(-0.6337, 0.7384, 0.0566, 1.2465, -0.2098, 2.0432, -1.0666, 1.7536, -0.1845),
c(0.9720, -0.6342, -2.5292, -1.8116, -2.0078, -1.0666, 4.0882, -1.3587, 0.7287),
c(-1.2158, 1.4291, 0.4776, 2.0076, -0.1715, 1.7536, -1.3587, 2.8789, 0.4094),
c(-1.3219, -0.4723, -0.4552, -0.3377, -0.3952, -0.1845, 0.7287, 0.4094, 4.8406))
library(MASS)
```

```
##
## Attaching package: 'MASS'
```

```
## The following object is masked from 'package:CVXR':
##
## huber
```

```
X <- mvrnorm(m, mu, Sigma)
X <- cbind(rep(1, m), X)
trueBeta <- c(0, 0.8, 0, 1, 0.2, 0, 0.4, 1, 0, 0.7)
y <- X %*% trueBeta + rnorm(m, 0, s)
```

```
beta <- Variable(n)
objective <- Minimize(sum_squares(y - X %*% beta))
prob3.1 <- Problem(objective)
```

Here, `y`

is the response, `X`

is the matrix of predictors, `n`

is the number of predictors, and `beta`

is a vector of coefficients on the predictors. The Ordinary Least-Squares (OLS) solution for `beta`

minimizes the \(l_2\)-norm of the residuals (i.e. the root-mean-squared error). As we can see below, `CVXR`

’s solution matches the solution obtained by using `lm`

.

```
CVXR_solution3.1 <- solve(prob3.1)
lm_solution3.1 <- lm(y ~ 0 + X)
```

Obviously, if all you want to do is least-squares linear regression, you should simply use `lm`

. The chief advantage of `CVXR`

is its flexibility, as we will demonstrate in the rest of this section.

### Non-Negative Least Squares

Looking at example 3.1, you may notice that the OLS regression problem has an objective, but no constraints. In many contexts, we can greatly improve our model by constraining the solution to reflect our prior knowledge. For example, we may know that the coefficients `beta`

must be non-negative.

```
prob3.2 <- prob3.1
constraints(prob3.2) <- list(beta >= 0)
solution3.2 <- solve(prob3.2)
```

As we can see in the figure above, adding that one constraint produced a massive improvement in the accuracy of the estimates. Not only are the non-negative least-squares estimates much closer to the true coefficients than the OLS estimates, they have even managed to recover the correct sparsity structure in this case.

Like with OLS, there are already R packages available which implement non-negative least squares, such as `nnls`

. But that is actually an excellent demonstration of the power of `CVXR`

: A single line of code here, namely `prob3.2$constraints <- list(beta >= 0)`

, is doing the work of an entire package.

### Support Vector Classifiers

Another common statistical tool is the support vector classifier (SVC). The SVC is an affine function (hyperplane) that separates two sets of points by the widest margin. When the sets are not linearly separable, the SVC is determined by a trade-off between the width of the margin and the number of points that are misclassified.

For the binary case, where the response \(y_i \in \{-1,1\}\), the SVC is obtained by solving

\[ \begin{array}{ll} \mbox{minimize} & \frac{1}{2}\Vert\beta\Vert^2 + C\sum_{i=1}^m \xi_i \\ \mbox{subject to} & \xi_i \geq 0, \quad y_i(x_i^T\beta + \beta_0) \geq 1 - \xi_i, \quad i = 1,\ldots,m \end{array} \]

with variables \((\beta,\xi)\). Here, \(\xi\) is the amount by which a point can violate the separating hyperplane, and \(C > 0\) is a user-chosen penalty on the total violation. As \(C\) increases, fewer misclassifications will be allowed.

Below, we fit a SVC in `CVXR`

with \(C = 10\).

```
## Generate data
set.seed(10)
n <- 2
m <- 50
X <- matrix(rnorm(m*n), nrow = m, ncol = n)
y <- c(rep(-1, m/2), rep(1, m/2))
X[y == 1,] = X[y == 1,] + 1
```

```
## Define variables
cost <- 10
beta0 <- Variable()
beta <- Variable(n)
slack <- Variable(m)
# Form problem
objective <- (1/2) * sum_squares(vstack(beta, beta0)) + cost * sum(slack)
constraints <- list(y * (X %*% beta + beta0) >= 1 - slack, slack >= 0)
prob3.3 <- Problem(Minimize(objective), constraints)
solution3.3 <- solve(prob3.3)
```

## Disciplined Convex Programming (DCP)

Disciplined convex programming (DCP) is a system for constructing mathematical expressions with known sign and curvature from a given library of base functions. `CVXR`

uses DCP to ensure that the specified optimization problems are convex.

The user may find it helpful to read about how the DCP rules are applied in other languages such as Python, Matlab, and Julia.

`CVXR`

implements the same rules, and a short introduction is available here. The set of DCP functions are described in Function Reference.

## Session Info

`sessionInfo()`

```
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin17.7.0 (64-bit)
## Running under: macOS 10.14.1
##
## Matrix products: default
## BLAS/LAPACK: /usr/local/Cellar/openblas/0.3.3/lib/libopenblasp-r0.3.3.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices datasets utils methods base
##
## other attached packages:
## [1] MASS_7.3-51.1 CVXR_0.99-2 tidyr_0.8.2 ggplot2_3.1.0
##
## loaded via a namespace (and not attached):
## [1] gmp_0.5-13.2 Rcpp_0.12.19 pillar_1.3.0
## [4] compiler_3.5.1 plyr_1.8.4 bindr_0.1.1
## [7] R.methodsS3_1.7.1 R.utils_2.7.0 tools_3.5.1
## [10] bit_1.1-14 digest_0.6.18 lattice_0.20-35
## [13] evaluate_0.12 tibble_1.4.2 gtable_0.2.0
## [16] pkgconfig_2.0.2 rlang_0.3.0.1 Matrix_1.2-15
## [19] yaml_2.2.0 blogdown_0.9.2 xfun_0.4
## [22] bindrcpp_0.2.2 Rmpfr_0.7-1 ECOSolveR_0.4
## [25] withr_2.1.2 stringr_1.3.1 dplyr_0.7.7
## [28] knitr_1.20 bit64_0.9-7 rprojroot_1.3-2
## [31] grid_3.5.1 tidyselect_0.2.5 glue_1.3.0
## [34] R6_2.3.0 rmarkdown_1.10 bookdown_0.7
## [37] purrr_0.2.5 magrittr_1.5 scs_1.1-1
## [40] backports_1.1.2 scales_1.0.0 htmltools_0.3.6
## [43] assertthat_0.2.0 colorspace_1.3-2 labeling_0.3
## [46] stringi_1.2.4 lazyeval_0.2.1 munsell_0.5.0
## [49] crayon_1.3.4 R.oo_1.22.0
```

## Source

## References

Boyd, S., and L. Vandenberghe. 2004. *Convex Optimization*. Cambridge University Press.