---
title: "Introduction to CVXR"
author: "Anqi Fu and Balasubramanian Narasimhan"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to CVXR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

`CVXR` is an R package that provides an object-oriented modeling
language for convex optimization, similar to
[CVXPY](https://www.cvxpy.org/) for Python. It allows you to
formulate and solve convex optimization problems in a natural
mathematical syntax.

## A Simple Example: Least Squares

Consider a simple linear regression problem where we want to estimate
parameters using a least squares criterion.

We generate synthetic data where we know the true model:

$$ Y = X\beta + \epsilon $$

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

```{r}
set.seed(123)
n <- 100
p <- 10
beta <- -4:5

X <- matrix(rnorm(n * p), nrow = n)
Y <- X %*% beta + rnorm(n)
```

Using base R, we can estimate $\beta$ via `lm`:

```{r}
ls.model <- lm(Y ~ 0 + X)
```

### The CVXR formulation

The same problem can be expressed as:

$$
\underset{\beta}{\text{minimize}} \quad \|Y - X\beta\|_2^2
$$

In CVXR, this translates directly:

```{r}
library(CVXR)

betaHat <- Variable(p)
objective <- Minimize(sum((Y - X %*% betaHat)^2))
problem <- Problem(objective)
result <- psolve(problem) ## use default solver
```

The optimal value and estimated coefficients:

```{r}
cat("Optimal value:", result, "\n")
cbind(CVXR = round(value(betaHat), 3),
      lm   = round(coef(ls.model), 3))
```

## Adding Constraints

The real power of CVXR is the ability to add constraints easily.

### Nonnegative Least Squares

Suppose we know the $\beta$s should be nonnegative:

```{r}
problem <- Problem(objective, constraints = list(betaHat >= 0))
result <- psolve(problem, solver = "CLARABEL")
round(value(betaHat), 3)
```

### Custom Constraints

Now suppose $\beta_2 + \beta_3 \le 0$ and all other $\beta$s are
nonnegative:

```{r}
A <- matrix(c(0, 1, 1, rep(0, 7)), nrow = 1)
B <- diag(c(1, 0, 0, rep(1, 7)))

constraint1 <- A %*% betaHat <= 0
constraint2 <- B %*% betaHat >= 0

problem <- Problem(objective, constraints = list(constraint1, constraint2))
result <- psolve(problem, solver = "CLARABEL", verbose = TRUE) ## verbose = TRUE for details
round(value(betaHat), 3)
```

This demonstrates the chief advantage of CVXR: _flexibility_. Users
can quickly modify and re-solve a problem, making the package ideal
for prototyping new statistical methods. Its syntax is simple and
mathematically intuitive.

## Available Solvers

CVXR 1.8.2 supports 15 solvers, both open source and commercial:
Clarabel, SCS, OSQP, HiGHS, MOSEK, Gurobi, GLPK, GLPK_MI, ECOS,
ECOS_BB, CPLEX, CVXOPT, PIQP, SCIP, and XPRESS.

```{r}
installed_solvers()
```

You can specify a solver explicitly:

```{r, eval = FALSE}
psolve(problem, solver = "CLARABEL")
```

## What's New in 1.8.2

- **15 solvers**: SCIP and XPRESS added (up from 13 in 1.8.1).
- **Element-wise matrix indexing**: Constrain specific entries of a
  matrix variable using R's native indexing idioms:
  ```r
  ind <- which(!is.na(Rmiss), arr.ind = TRUE)
  prob <- Problem(Minimize(obj), list(X[ind] == Rmiss[ind]))
  ```
  Also supports logical matrix indexing (`X[mask]`) and linear
  integer indexing (`X[c(1, 5, 9)]`).
- **CVXPY 1.8.2 parity**: All applicable bug fixes ported.
- **Unified solver options** via `solver_opts()`.

See `NEWS.md` for full details.

## Further Reading

- The [CVXR website](https://cvxr.rbind.io) has many worked examples
- The [CVXPY documentation](https://www.cvxpy.org/) covers the
  underlying mathematical framework
- [The published paper](https://doi.org/10.18637/jss.v094.i14): Fu, Narasimhan, and Boyd (2020). "CVXR: An R Package for Disciplined
  Convex Optimization." _Journal of Statistical Software_, 94(14), DOI:10.18637/jss.v094.i14.

## Session Info

```{r}
sessionInfo()
```