`CVXR`

?`CVXR`

is an R package that provides an object-oriented modeling
language for convex optimization, similar to `CVX`

, `CVXPY`

, `YALMIP`

,
and `Convex.jl`

. It allows the user to formulate convex optimization
problems in a natural mathematical syntax rather than the restrictive
standard form required by most solvers. The user specifies an
objective and set of constraints by combining constants, variables,
and parameters using a library of functions with known mathematical
properties. `CVXR`

then applies
signed disciplined convex programming (DCP)
to verify the problemâ€™s convexity. Once verified, the problem is
converted into standard conic form using graph implementations and
passed to a cone solver such
as ECOS
or SCS.

The paper by @fu:naras:boyd:2017 is the main reference. Further documentation, along with a number of tutorial examples, is also available on the CVXR website.

Below we provide a simple example to get you started.

Consider a simple linear regression problem where it is desired to estimate a set of parameters using a least squares criterion.

We generate some synthetic data where we know the model completely, that is

\[ Y = X\beta + \epsilon \]

where \(Y\) is a \(100\times 1\) vector, \(X\) is a \(100\times 10\) matrix, \(\beta = [-4,\ldots ,-1, 0, 1, \ldots, 5]\) is a \(10\times 1\) vector, and \(\epsilon \sim N(0, 1)\).

```
set.seed(123)
n <- 100
p <- 10
beta <- -4:5 # beta is just -4 through 5.
X <- matrix(rnorm(n * p), nrow=n)
colnames(X) <- paste0("beta_", beta)
Y <- X %*% beta + rnorm(n)
```

Given the data \(X\) and \(Y\), we can estimate the \(\beta\) vector using
`lm`

function in R that fits a standard regression model.

```
ls.model <- lm(Y ~ 0 + X) # There is no intercept in our model above
m <- data.frame(ls.est = coef(ls.model))
rownames(m) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(m)
```

ls.est | |
---|---|

\(\beta_{1}\) | -3.9196886 |

\(\beta_{2}\) | -3.0117048 |

\(\beta_{3}\) | -2.1248242 |

\(\beta_{4}\) | -0.8666048 |

\(\beta_{5}\) | 0.0914658 |

\(\beta_{6}\) | 0.9490454 |

\(\beta_{7}\) | 2.0764700 |

\(\beta_{8}\) | 3.1272275 |

\(\beta_{9}\) | 3.9609565 |

\(\beta_{10}\) | 5.1348845 |

These are the least-squares estimates and can be seen to be reasonably close to the original \(\beta\) values -4 through 5.

`CVXR`

formulationThe `CVXR`

formulation states the above as an optimization problem:

\[
\begin{array}{ll}
\underset{\beta}{\mbox{minimize}} & \|y - X\beta\|_2^2,
\end{array}
\]
which directly translates into a problem that `CVXR`

can solve as shown
in the steps below.

- Step 0. Load the
`CVXR`

library

```
suppressWarnings(library(CVXR, warn.conflicts=FALSE))
```

- Step 1. Define the variable to be estimated

```
betaHat <- Variable(p)
```

- Step 2. Define the objective to be optimized

```
objective <- Minimize(sum((Y - X %*% betaHat)^2))
```

Notice how the objective is specified using functions such as `sum`

,
`*%*`

and `^`

, that are familiar to R users despite that fact that
`betaHat`

is no ordinary R expression but a `CVXR`

expression.

- Step 3. Create a problem to solve

```
problem <- Problem(objective)
```

- Step 4. Solve it!

```
result <- solve(problem)
```

- Step 5. Extract solution and objective value

```
## Objective value: 97.847586
```

We can indeed satisfy ourselves that the results we get matches that
from `lm`

.

```
m <- cbind(coef(ls.model), result$getValue(betaHat))
colnames(m) <- c("lm est.", "CVXR est.")
rownames(m) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(m)
```

lm est. | CVXR est. | |
---|---|---|

\(\beta_{1}\) | -3.9196886 | -3.9196886 |

\(\beta_{2}\) | -3.0117048 | -3.0117048 |

\(\beta_{3}\) | -2.1248242 | -2.1248242 |

\(\beta_{4}\) | -0.8666048 | -0.8666048 |

\(\beta_{5}\) | 0.0914658 | 0.0914658 |

\(\beta_{6}\) | 0.9490454 | 0.9490454 |

\(\beta_{7}\) | 2.0764700 | 2.0764700 |

\(\beta_{8}\) | 3.1272275 | 3.1272275 |

\(\beta_{9}\) | 3.9609565 | 3.9609565 |

\(\beta_{10}\) | 5.1348845 | 5.1348845 |

On the surface, it appears that we have replaced one call to `lm`

with
at least five or six lines of new R code. On top of that, the code
actually runs slower, and so it is not clear what was really achieved.

So suppose we knew for a fact that the $\beta$s were nonnegative and
we wish to take this fact into account. This
is
nonnegative least squares regression and
`lm`

would no longer do the job.

In `CVXR`

, the modified problem merely requires the addition of a constraint to the
problem definition.

```
problem <- Problem(objective, constraints = list(betaHat >= 0))
result <- solve(problem)
m <- data.frame(CVXR.est = result$getValue(betaHat))
rownames(m) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(m)
```

CVXR.est | |
---|---|

\(\beta_{1}\) | 0.0000000 |

\(\beta_{2}\) | 0.0000000 |

\(\beta_{3}\) | 0.0000000 |

\(\beta_{4}\) | 0.0000000 |

\(\beta_{5}\) | 1.2374488 |

\(\beta_{6}\) | 0.6234665 |

\(\beta_{7}\) | 2.1230663 |

\(\beta_{8}\) | 2.8035640 |

\(\beta_{9}\) | 4.4448016 |

\(\beta_{10}\) | 5.2073521 |

We can verify once again that these values are comparable to those obtained from another R package, say nnls.

```
if (requireNamespace("nnls", quietly = TRUE)) {
nnls.fit <- nnls::nnls(X, Y)$x
} else {
nnls.fit <- rep(NA, p)
}
```

```
m <- cbind(result$getValue(betaHat), nnls.fit)
colnames(m) <- c("CVXR est.", "nnls est.")
rownames(m) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(m)
```

CVXR est. | nnls est. | |
---|---|---|

\(\beta_{1}\) | 0.0000000 | 0.0000000 |

\(\beta_{2}\) | 0.0000000 | 0.0000000 |

\(\beta_{3}\) | 0.0000000 | 0.0000000 |

\(\beta_{4}\) | 0.0000000 | 0.0000000 |

\(\beta_{5}\) | 1.2374488 | 1.2374488 |

\(\beta_{6}\) | 0.6234665 | 0.6234665 |

\(\beta_{7}\) | 2.1230663 | 2.1230663 |

\(\beta_{8}\) | 2.8035640 | 2.8035640 |

\(\beta_{9}\) | 4.4448016 | 4.4448016 |

\(\beta_{10}\) | 5.2073521 | 5.2073521 |

As you no doubt noticed, we have done nothing that other R packages could not do.

So now suppose that we know, for some extraneous reason, that the sum of \(\beta_2\) and \(\beta_3\) is nonpositive and but all other $\beta$s are nonnegative.

It is clear that this problem would not fit into any standard
package. But in `CVXR`

, this is easily done by adding a few
constraints.

To express the fact that \(\beta_2 + \beta_3\) is nonpositive, we construct a row matrix with zeros everywhere, except in positions 2 and 3 (for \(\beta_2\) and \(\beta_3\) respectively).

```
A <- matrix(c(0, 1, 1, rep(0, 7)), nrow = 1)
colnames(A) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(A)
```

\(\beta_{1}\) | \(\beta_{2}\) | \(\beta_{3}\) | \(\beta_{4}\) | \(\beta_{5}\) | \(\beta_{6}\) | \(\beta_{7}\) | \(\beta_{8}\) | \(\beta_{9}\) | \(\beta_{10}\) |
---|---|---|---|---|---|---|---|---|---|

0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

The sum constraint is nothing but \[ A\beta <= 0 \]

which we express in R as

```
constraint1 <- A %*% betaHat <= 0
```

*NOTE*: The above constraint can also be expressed simply as

```
constraint1 <- betaHat[2] + betaHat[3] <= 0
```

but it is easier working with matrices in general with `CVXR`

.

For the nonnegativity for rest of the variables, we construct a \(10\times 10\) matrix \(A\) to have 1's along the diagonal everywhere except rows 2 and 3 and zeros everywhere.

```
B <- diag(c(1, 0, 0, rep(1, 7)))
colnames(B) <- rownames(B) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(B)
```

\(\beta_{1}\) | \(\beta_{2}\) | \(\beta_{3}\) | \(\beta_{4}\) | \(\beta_{5}\) | \(\beta_{6}\) | \(\beta_{7}\) | \(\beta_{8}\) | \(\beta_{9}\) | \(\beta_{10}\) | |
---|---|---|---|---|---|---|---|---|---|---|

\(\beta_{1}\) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

\(\beta_{2}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

\(\beta_{3}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

\(\beta_{4}\) | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |

\(\beta_{5}\) | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |

\(\beta_{6}\) | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |

\(\beta_{7}\) | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |

\(\beta_{8}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |

\(\beta_{9}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |

\(\beta_{10}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |

The constraint for positivity is \[ B\beta > 0 \]

which we express in R as

```
constraint2 <- B %*% betaHat >= 0
```

Now we are ready to solve the problem just as before.

```
problem <- Problem(objective, constraints = list(constraint1, constraint2))
result <- solve(problem)
```

And we can get the estimates of \(\beta\).

```
m <- data.frame(CVXR.soln = result$getValue(betaHat))
rownames(m) <- paste0("$\\beta_{", 1:p, "}$")
knitr::kable(m)
```

CVXR.soln | |
---|---|

\(\beta_{1}\) | 0.0000000 |

\(\beta_{2}\) | -2.8446952 |

\(\beta_{3}\) | -1.7109771 |

\(\beta_{4}\) | 0.0000000 |

\(\beta_{5}\) | 0.6641308 |

\(\beta_{6}\) | 1.1781109 |

\(\beta_{7}\) | 2.3286139 |

\(\beta_{8}\) | 2.4144893 |

\(\beta_{9}\) | 4.2119052 |

\(\beta_{10}\) | 4.9483245 |

This demonstrates the chief advantage of `CVXR`

: flexibility. Users
can quickly modify and re-solve a problem, making our package ideal
for prototyping new statistical methods. Its syntax is simple and
mathematically intuitive. Furthermore, `CVXR`

combines seamlessly with
native R code as well as several popular packages, allowing it to be
incorporated easily into a larger analytical framework. The user is
free to construct statistical estimators that are solutions to a
convex optimization problem where there may not be a closed form
solution or even an implementation. Such solutions can then be
combined with resampling techniques like the bootstrap to estimate
variability.