Getting Faster Results

Warning

The solution described below is useful when you mathematically know a problem is DCP-compliant and none of your data inputs will change the nature of the problem. We recommend that users check the DCP-compliance of a problem (via a call to is_dcp(prob) for example) at least once to ensure this is the case. Not verifying DCP-compliance may result in garbage!

Introduction

As was remarked in the introduction to CVXR, its chief advantage is flexibility: you can specify a problem in close to mathematical form and CVXR solves it for you, if it can. Behind the scenes, CVXR compiles the domain specific language and verifies the convexity of the problem before sending it off to solvers. If the problem violates the rules of Disciplined Convex Programming it is rejected.

Therefore, it is generally slower than tailor-made solutions to a given problem.

An Example

To understand the speed issues, let us consider the global warming data from the Carbon Dioxide Information Analysis Center (CDIAC) again. The data points are the annual temperature anomalies relative to the 1961–1990 mean. We will fit the nearly-isotonic approximation \(\beta \in {\mathbf R}^m\) by solving

\[ \begin{array}{ll} \underset{\beta}{\mbox{Minimize}} & \frac{1}{2}\sum_{i=1}^m (y_i - \beta_i)^2 + \lambda \sum_{i=1}^{m-1}(\beta_i - \beta_{i+1})_+, \end{array} \] where \(\lambda \geq 0\) is a penalty parameter and \(x_+ =\max(x,0)\).

This can be solved as follows.

suppressMessages(suppressWarnings(library(CVXR)))
data(cdiac)
y <- cdiac$annual
m <- length(y)
lambda <- 0.44
beta <- Variable(m)
obj <- 0.5 * sum((y - beta)^2) + lambda * sum(pos(diff(beta)))
prob <- Problem(Minimize(obj))
soln <- solve(prob)
betaHat <- soln$getValue(beta)

This is the recommended way to solve a problem.

However, suppose we wished to construct bootstrap confidence intervals for the estimate using 100 resamples. It is clear that this computation time can quickly become limiting .

Below, we show how one can get at the problem data and directly call a solver to get faster results.

Profile the code

Profiling a single fit to the model is useful to figure out where most of the time is spent.

library(profvis)
suppressMessages(suppressWarnings(library(CVXR)))
data(cdiac)
y <- cdiac$annual
profvis({
    beta <- Variable(m)
    obj <- Minimize(0.5 * sum((y - beta)^2) + lambda * sum(pos(diff(beta))))
    prob <- Problem(obj)
    soln <- solve(prob)
    betaHat <- soln$getValue(beta)
})