CVXPY DCP errors

CVXPY is a popular tool to model and solve convex optimization problems. Sometimes, it throws a "DCP" error, even for convex problems. DCP stands for Disciplined Convex Programming, the underlying framework for working with guaranteed convex models. The error says: I can not verify this is convex.  

Here are some small examples of convex objectives (under minimization) one would expect to work.

Objective	CVXPY code	Result	Notes
\[\color{darkred}x^T \color{darkred}x\]	x.T@x	DCP error	print shows minimize x@x, i.e. transpose is dropped
\[\color{darkred}x^T \color{darkred}x\]	x@x	DCP error
\[\color{darkred}x^T \color{darkred}x\]	cp.sum_squares(x)		transformed into quad_over_lin(x, 1.0)
\[\color{darkred}x^T \color{darkblue}Q \color{darkred}x\]	x.T@Q@x		transformed into QuadForm(x,Q)
\[\color{darkred}y:=\color{darkred}x-\color{darkblue}p\]\[\color{darkred}x^T \color{darkblue}Q \color{darkred}y\]	y=x-p x.T@Q@y	DCP error
\[\color{darkred}x^T \color{darkblue}Q \color{darkred}x - \color{darkred}x^T \color{darkblue}Q \color{darkblue}p\]	x.T@Q@x - x.T@Q@p		first term transformed into QuadForm(x,Q)

Not everything makes sense to me. I am not sure why x.T@x is not properly recognized, but x.T@Q@x is.

Fix non positive semi-definite covariance matrix in R

In [1], the following R code is proposed to solve a Mean-Variance portfolio optimization problem:

set.seed(123)

nO <- 10L  ## number of observations
nA <- 100L  ## number of assets
mData <- array(rnorm(nO * nA, sd = 0.05), 
               dim = c(nO, nA)) #Creating sample stock observations

library("quadprog")
aMat <- array(1, dim = c(1,nA))
bVec <- 1
zeros <- array(0, dim = c(nA,1))

solQP <- solve.QP(cov(mData), zeros, t(aMat), bVec, meq = 1) #Minimize optimization
solQP$solution

I added the line that explicitly sets the seed to make the runs reproducible. 

The model is a bit difficult to decipher in this R code, but it is a simplified traditional MV model:

Mean-Variance Portfolio Optimization Model
\[\begin{align} \min\> & \color{darkred}x'\color{darkblue}Q\color{darkred}x - \color{darkblue}\lambda \bar{\color{darkblue}r}'\color{darkred}x\\ & \sum_i \color{darkred}x_i = 1 \\ & \color{darkred}x_i \ge 0\end{align}\]

Instead of having the expected returns in the objective, we also see versions of this model with a linear constraint that puts a minimum value on the expected return of the portfolio. 

In the example, we have \(\lambda=0\), so the objective only has a quadratic term to minimize the risk. We don't care about returns. 

Note that the data only has 10 observations but 100 instruments. This situation can often lead to some problems. Indeed when we run this, we see:

> solQP <- solve.QP(cov(mData), zeros, t(aMat), bVec, meq = 1) #Minimize optimization
Error in solve.QP(cov(mData), zeros, t(aMat), bVec, meq = 1) : 
  matrix D in quadratic function is not positive definite!
> 

As a side note: the data is just completely random, so we cannot expect useful solutions. However, at least, we should be able to produce some solutions.

 

Solution #1

The covariance matrix is in theory positive-semi definite. We can see this from using the definition of the covariance [2]. However, due to floating-point inaccuracies, we can end up with a covariance matrix that is just slightly non-positive semi-definite. We can see this by inspecting the eigenvalues:

> Q <- cov(mData)
> eigen(Q)$values
  [1]  3.977635e-02  3.789523e-02  3.258164e-02  3.093782e-02  2.605389e-02  2.466463e-02  2.127899e-02
  [8]  1.897506e-02  1.627900e-02  1.871945e-17  1.454315e-17  8.114384e-18  5.109196e-18  4.392427e-18
 [15]  3.680766e-18  2.728485e-18  2.649302e-18  2.192073e-18  1.937634e-18  1.914487e-18  1.856811e-18
 [22]  1.449321e-18  1.391434e-18  1.211216e-18  1.196721e-18  1.127119e-18  9.776146e-19  9.670197e-19
 [29]  8.257462e-19  7.572554e-19  6.283522e-19  6.030820e-19  4.056541e-19  3.610128e-19  3.604800e-19
 [36]  2.770106e-19  2.530797e-19  2.263314e-19  2.133854e-19  2.089700e-19  2.024130e-19  2.006560e-19
 [43]  1.688435e-19  9.412675e-20  7.033666e-20  7.012531e-20  5.942729e-20  5.673501e-20  4.444257e-20
 [50]  4.259787e-20  3.124233e-21 -2.017011e-20 -2.269857e-20 -3.034668e-20 -3.075493e-20 -3.262642e-20
 [57] -4.345710e-20 -6.351674e-20 -7.347496e-20 -7.739535e-20 -8.065131e-20 -9.318415e-20 -1.004645e-19
 [64] -1.048992e-19 -1.185271e-19 -1.315539e-19 -1.362383e-19 -1.376705e-19 -1.456396e-19 -1.697860e-19
 [71] -1.700317e-19 -1.745068e-19 -1.794480e-19 -1.835296e-19 -1.960933e-19 -2.017457e-19 -2.332259e-19
 [78] -2.377082e-19 -3.232830e-19 -3.468225e-19 -3.844707e-19 -4.310383e-19 -4.467493e-19 -4.540707e-19
 [85] -4.977900e-19 -6.159103e-19 -6.338311e-19 -6.484161e-19 -6.489508e-19 -6.551306e-19 -6.664641e-19
 [92] -7.276688e-19 -7.422977e-19 -7.927766e-19 -1.106209e-18 -1.707077e-18 -1.742870e-18 -2.465438e-18
 [99] -1.859346e-17 -3.473638e-17

Most of the eigenvalues are basically zero with quite a few just a bit negative. A positive-definite matrix would have only positive values. A positive semi-definite (PSD) matrix would allow positive and zero eigenvalues but no negative ones.

There are algorithms to find a close-by positive-definite matrix. In R this is available as nearPD() in the Matrix package. Let's try this out:

> Q <- nearPD(cov(mData))$mat
> solQP <- solve.QP(Q, zeros, t(aMat), bVec, meq = 1) 
> solQP$solution
  [1]  0.007861808  0.004649581  0.012677143  0.009135964  0.011926439  0.013624409  0.007916747  0.004307772
  [9]  0.011808486  0.011085302  0.010229029  0.013080240  0.010085504  0.010998879  0.011099784  0.008520520
 [17]  0.007451730  0.012424371  0.012677330  0.007013515  0.009070013  0.010422999  0.011041407  0.008768456
 [25]  0.007594051  0.010845828  0.008091472  0.010219700  0.017100764  0.005583947  0.011981802  0.007366402
 [33]  0.013770522  0.009063828  0.010744898  0.007836608  0.009580356  0.007065588  0.008219370  0.010592597
 [41]  0.011775598  0.011112051  0.005555992  0.014596684  0.016019481  0.008870919  0.010940604  0.005662200
 [49]  0.007056997  0.012274299  0.006510091  0.012764160  0.008258651  0.009945869  0.009891680  0.007318763
 [57]  0.014960558  0.008255708  0.009655097  0.010054760  0.007309219  0.008761360  0.011616297  0.010300254
 [65]  0.012043733  0.011963495  0.009991699  0.010353212  0.013896332  0.011619403  0.007160011  0.007202327
 [73]  0.008530819  0.010918658  0.004754215  0.008877062  0.008526937  0.010251250  0.008190446  0.006754629
 [81]  0.013209816  0.014849748  0.012213860  0.008888433  0.018188014  0.009756690  0.014786487  0.014893706
 [89]  0.009737739  0.013256532  0.009990737 -0.001155889  0.009264992  0.007733719  0.014505131  0.011040198
 [97]  0.012859850  0.008987104  0.008458713  0.004497735

Much better. When printing the eigenvalues of the new Q matrix, we see:

> eigen(Q)$values
  [1] 3.977635e-02 3.789522e-02 3.258164e-02 3.093782e-02 2.605389e-02 2.466462e-02 2.127899e-02 1.897506e-02
  [9] 1.627900e-02 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [17] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [25] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [33] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [41] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [49] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [57] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [65] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10
 [73] 3.977635e-10 3.977635e-10 3.977635e-10 3.977635e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10
 [81] 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10
 [89] 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10 3.977634e-10
 [97] 3.977634e-10 3.977634e-10 3.977633e-10 3.977633e-10

All small (and negative) eigenvalues are now mapped to some tolerance. This means the Q matrix is now numerically positive definite.

Solution #2

A simple DIY approach would be to fix the Q matrix by adding a small number to the diagonal. The following code shows a small perturbation is enough:

> Q <- cov(mData) + diag(1e-6,nA)

Some solvers do something like this automatically.

Fixing the model

Let's first fix our model a bit:

• We want to pay attention to returns. So let's add a constraint that the expected portfolio returns are at least 0.035. This also has the effect that just a few assets will be in our solution. Which makes them easier to compare. The constraint can be stated as \[\sum_i \bar{r}_i \ge 0.035\] where \[\bar{r}_i = \frac{\sum_t r_{t,i}}{T}\] are the mean returns. Here \(T\) indicates the number of observations, and \(r_{t,i}\) are the returns.
• We also want to explicitly state that variables are non-negative: \(x_i \ge 0\). This is called "no short-selling",  

The portfolio model becomes:

Portfolio Model algebraic formulation	matrix notation
\[\begin{align} \min & \sum_{i,j} \color{darkblue}q_{i,j} \color{darkred}x_i \color{darkred}x_j \\ & \sum_i \color{darkred}x_i = 1 \\ & \sum_i \bar{\color{darkblue}r}_i \color{darkred}x_i \ge 0.035\\ & \color{darkred}x_i \ge 0\end{align}\]	\[\begin{align} \min \>& \color{darkred}x' \color{darkblue}Q \color{darkred}x \\ & \color{darkblue}e'\color{darkred}x = 1 \\ & \bar{\color{darkblue}r}' \color{darkred}x \ge 0.035\\ & \color{darkred}x \ge 0\end{align}\]

QuadProg does not have facilities for bounds on the variables so we need to form explicit constraints.  Its input-model looks like: \[\begin{align} \min \> & -d'b + \frac{1}{2} b'Db \\ & A'b\ge b_0\end{align}\] where the first meq constraints are equalities. The factor 0.5 is a relic from quadratic optimization (it makes the gradient a bit easier). If the model includes linear objective coefficients, it is important to take this 0.5 factor into account. Our model has only quadratic objective coefficients, so we don't care in this case. All in all, this is often a bit of an inconvenient input format. 

R trick. In the examples below we want to display a sparse vector. A simple trick using names can help here. 

> v <- c(0,0,3,0,0,0,1,0,0,0,2,0,0,0,0)
> v
 [1] 0 0 3 0 0 0 1 0 0 0 2 0 0 0 0
> names(v) <- paste("v",1:length(v),sep="")
> v[v>0]
 v3  v7 v11 
  3   1   2

The R code for our optimization model can look like:

> Q <- nearPD(cov(mData))$mat
> aMat <- rbind(
+     c(rep(1,nA)),     # sum(x)=1
+     colMeans(mData),  # sum(mu(i)*x(i))>=0.035   
+     diag(1,nA))       # for bounds x(i)>=0
> bVec <- c(1,
+          0.035,
+          rep(0,nA))
> solQP <- solve.QP(Q, zeros, t(aMat), bVec, meq = 1) 
> xsol <- solQP$solution
> names(xsol) = paste("x",1:100,sep="")
> xsol[xsol>1e-4]
      x30       x75       x88 
0.1465990 0.5314489 0.3219522 

Another issue with QuadProg is that requires a strictly positive definite matrix. This also means it does not allow linear variables: all variables need appear to quadratically in the objective (and without all zero quadratic coefficients). After all, if the quadratic coefficient matrix would look like \[D = \left[\begin{array}{c|c}Q & 0 \\ \hline 0 & 0 \end{array} \right]\] we would end up with a matrix that is not strictly positive-definite.

Often, for these types of models, CVXR is a more agreeable tool.  Here is the same model in CVXR, and solved with the OSQP QP solver:

> library(CVXR)
> x <- Variable(nA)
> w <- Variable(nO)
> prob <- Problem(
+     Minimize(quad_form(x,cov(mData))),
+     list(sum(x) == 1,
+          t(colMeans(mData)) %*% x >= 0.035,
+          x >= 0)
+     )
> result <- solve(prob, verbose=T)
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 100, constraints m = 102
          nnz(P) + nnz(A) = 5350
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter  objective    pri res    dua res    rho        time
   1   0.0000e+00   1.00e+00   1.00e+02   1.00e-01   2.88e-03s
 200   8.7535e-04   5.96e-04   1.36e-06   1.75e-02   9.38e-03s
 400   9.2707e-04   5.09e-04   1.16e-06   1.75e-02   1.41e-02s
 600   9.7703e-04   4.36e-04   9.94e-07   1.75e-02   1.99e-02s
 800   1.0239e-03   3.72e-04   8.50e-07   1.75e-02   2.58e-02s
1000   1.0671e-03   3.18e-04   7.27e-07   1.75e-02   3.11e-02s
1200   1.1062e-03   2.72e-04   6.21e-07   1.75e-02   3.76e-02s
1400   1.1413e-03   2.33e-04   5.31e-07   1.75e-02   4.35e-02s
1600   1.1726e-03   1.99e-04   4.54e-07   1.75e-02   4.90e-02s
1800   1.2001e-03   1.70e-04   3.88e-07   1.75e-02   5.51e-02s
2000   1.2243e-03   1.45e-04   3.32e-07   1.75e-02   6.20e-02s
2200   1.2455e-03   1.24e-04   2.84e-07   1.75e-02   6.81e-02s
2400   1.2639e-03   1.06e-04   2.43e-07   1.75e-02   7.33e-02s
2600   1.2799e-03   9.09e-05   2.07e-07   1.75e-02   7.83e-02s
2800   1.2938e-03   7.77e-05   1.77e-07   1.75e-02   8.51e-02s
3000   1.3058e-03   6.65e-05   1.52e-07   1.75e-02   9.20e-02s
3200   1.3161e-03   5.68e-05   1.30e-07   1.75e-02   1.00e-01s
3400   1.3251e-03   4.86e-05   1.11e-07   1.75e-02   1.08e-01s
3600   1.3327e-03   4.15e-05   9.48e-08   1.75e-02   1.21e-01s
3800   1.3394e-03   3.55e-05   8.11e-08   1.75e-02   1.28e-01s
4000   1.3450e-03   3.04e-05   6.93e-08   1.75e-02   1.35e-01s
4200   1.3499e-03   2.60e-05   5.92e-08   1.75e-02   1.43e-01s
4400   1.3541e-03   2.22e-05   5.07e-08   1.75e-02   1.50e-01s
4550   1.3568e-03   1.97e-05   4.50e-08   1.75e-02   1.56e-01s
plsh   1.3790e-03   3.33e-16   1.21e-15  ---------   1.58e-01s

status:               solved
solution polish:      successful
number of iterations: 4550
optimal objective:    0.0014
run time:             1.58e-01s
optimal rho estimate: 2.58e-02

> xsol <- result$getValue(x)
> names(xsol) = paste("x",1:100,sep="")
> xsol[xsol>1e-4]
      x30       x75       x88 
0.1465989 0.5314489 0.3219522 
>

Here we use our covariance matrix directly. Both CVXR and OSQP accept it, even though it has tiny negative numbers. However, the OSQP documentation warns that this may be a bit dangerous [3]:

We recommend the user to check the convexity of their problem before passing it to OSQP! If the user passes a non-convex problem we do not assure the solver will be able to detect it.

OSQP will try to detect non-convex problems by checking if the residuals diverge or if there are any issues in the initial factorization (if a direct method is used). It will detect non-convex problems when one or more of the eigenvalues of P are “clearly” negative, i.e., when P + sigma * I is not positive semidefinite. However, it might fail to detect non-convexity when P has slightly negative eigenvalues, i.e., when P + sigma * I is positive semidefinite and P is not.

Notice that CVXR will complain if the Q matrix is really non-positive semi-definite.

Solution #3

  

In [2] a different model is proposed for these situations where we have more instruments than observations.

Means-adjusted-returns Portfolio Optimization Model
\[\begin{align} \min\> & \sum_t \color{darkred}w_t^2 \\ &\color{darkred}w_t = \sum_i (\color{darkblue}r_{t,i} - \bar{\color{darkblue}r}_i)\color{darkred}x_i \\ & \sum_i \color{darkred}x_i = 1 \\ & \sum_i \bar{\color{darkblue}r}_i  \color{darkred}x_i  \ge 0.035\\ & \color{darkred}x_i \ge 0\end{align}\]

Some of the advantages of this model are:

• The quadratic objective is benign. It is sparse and safely positive definite.
• The data is smaller than a full-blown covariance matrix given that the number of observations is smaller than the number of instruments.
• The number of quadratic variables is smaller (again, because \(T<N\)).

 

This model is easily implemented in CVXR:

> library(CVXR)
> x <- Variable(nA)
> w <- Variable(nO)
> adjRet <- mData - matrix(colMeans(mData),nrow=nO,ncol=nA,byrow=T)
> prob <- Problem(
+     Minimize(sum_squares(w)),
+     list(w == adjRet %*% x,
+          sum(x)==1,
+          t(colMeans(mData)) %*% x >= 0.035,
+          x >= 0)
+     )
> result <- solve(prob, verbose=T)
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 110, constraints m = 112
          nnz(P) + nnz(A) = 1320
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter  objective    pri res    dua res    rho        time
   1   0.0000e+00   1.00e+00   1.00e+02   1.00e-01   1.15e-03s
 200   5.3515e-03   1.49e-03   4.04e-05   1.00e-01   4.33e-03s
 400   7.0787e-03   7.87e-04   6.50e-06   1.00e-01   8.41e-03s
 600   7.3521e-03   7.13e-04   5.90e-06   1.00e-01   1.23e-02s
 800   7.6413e-03   6.45e-04   5.34e-06   1.00e-01   1.61e-02s
1000   7.9373e-03   5.84e-04   4.84e-06   1.00e-01   1.93e-02s
1200   8.2334e-03   5.29e-04   4.38e-06   1.00e-01   2.27e-02s
1400   8.5244e-03   4.79e-04   3.96e-06   1.00e-01   2.64e-02s
1600   8.8067e-03   4.33e-04   3.59e-06   1.00e-01   3.00e-02s
1800   9.0778e-03   3.92e-04   3.25e-06   1.00e-01   3.34e-02s
2000   9.3358e-03   3.55e-04   2.94e-06   1.00e-01   3.69e-02s
2200   9.5798e-03   3.22e-04   2.66e-06   1.00e-01   4.05e-02s
2400   9.8091e-03   2.91e-04   2.41e-06   1.00e-01   4.43e-02s
2600   1.0024e-02   2.64e-04   2.18e-06   1.00e-01   4.81e-02s
2800   1.0224e-02   2.39e-04   1.98e-06   1.00e-01   5.12e-02s
3000   1.0410e-02   2.16e-04   1.79e-06   1.00e-01   5.48e-02s
3200   1.0582e-02   1.96e-04   1.62e-06   1.00e-01   5.84e-02s
3400   1.0741e-02   1.77e-04   1.47e-06   1.00e-01   6.19e-02s
3600   1.0887e-02   1.60e-04   1.33e-06   1.00e-01   6.55e-02s
3800   1.1022e-02   1.45e-04   1.20e-06   1.00e-01   6.90e-02s
4000   1.1145e-02   1.31e-04   1.09e-06   1.00e-01   7.53e-02s
4200   1.1259e-02   1.19e-04   9.85e-07   1.00e-01   8.18e-02s
4400   1.1362e-02   1.08e-04   8.92e-07   1.00e-01   8.95e-02s
4600   1.1457e-02   9.75e-05   8.07e-07   1.00e-01   9.44e-02s
4800   1.1544e-02   8.83e-05   7.31e-07   1.00e-01   1.04e-01s
5000   1.1623e-02   7.99e-05   6.62e-07   1.00e-01   1.09e-01s
5200   1.1695e-02   7.24e-05   5.99e-07   1.00e-01   1.13e-01s
5400   1.1761e-02   6.55e-05   5.42e-07   1.00e-01   1.20e-01s
5600   1.1821e-02   5.93e-05   4.91e-07   1.00e-01   1.26e-01s
5800   1.1876e-02   5.37e-05   4.45e-07   1.00e-01   1.31e-01s
6000   1.1925e-02   4.86e-05   4.02e-07   1.00e-01   1.37e-01s
6200   1.1970e-02   4.40e-05   3.64e-07   1.00e-01   1.45e-01s
6400   1.2011e-02   3.98e-05   3.30e-07   1.00e-01   1.58e-01s
6600   1.2049e-02   3.61e-05   2.99e-07   1.00e-01   1.72e-01s
6800   1.2082e-02   3.27e-05   2.70e-07   1.00e-01   1.81e-01s
7000   1.2113e-02   2.96e-05   2.45e-07   1.00e-01   1.90e-01s
7200   1.2141e-02   2.68e-05   2.22e-07   1.00e-01   1.98e-01s
7400   1.2166e-02   2.42e-05   2.01e-07   1.00e-01   2.01e-01s
7600   1.2189e-02   2.19e-05   1.82e-07   1.00e-01   2.05e-01s
7800   1.2210e-02   1.99e-05   1.64e-07   1.00e-01   2.08e-01s
plsh   1.2411e-02   6.66e-16   1.28e-15  ---------   2.09e-01s

status:               solved
solution polish:      successful
number of iterations: 7800
optimal objective:    0.0124
run time:             2.09e-01s
optimal rho estimate: 3.27e-01

> xsol <- result$getValue(x)
> names(xsol) = paste("x",1:100,sep="")
> xsol[xsol>1e-4]
      x30       x75       x88 
0.1465989 0.5314489 0.3219522 

Using QuadProg this is a much more difficult exercise. The quadratic coefficient matrix \(D\) will look like:  \[D = \left[\begin{array}{c|c}Q & 0 \\ \hline 0 & \varepsilon I \end{array} \right]\] The R code can look like:

> Q <- nearPD(cov(mData))$mat
> D <- rbind(
+   cbind(Q,matrix(0,nrow=nA,ncol=nO)),
+   cbind(matrix(0,nrow=nO,ncol=nA),diag(1e-6,nO))
+ )
> aMat <- rbind(
+   cbind(adjRet,diag(-1,nO)),     # sum(i,adjRet(t,i)*x(i))-w(t)=0
+   c(rep(1,nA),rep(0,nO)),        # sum(x)=1
+   c(colMeans(mData),rep(0,nO)),  # sum(mu(i)*x(i))>=0.035   
+   cbind(diag(1,nA),matrix(0,nrow=nA,ncol=nO)) # for bounds x(i)>=0
+ )
> bVec <- c(rep(0,nO),
+           1,
+           0.035,
+           rep(0,nA))
> zeros <- rep(0,nA+nO)
> solQP <- solve.QP(D, zeros, t(aMat), bVec, meq = nO+1) 
> xsol <- solQP$solution[1:nA]
> names(xsol) = paste("x",1:100,sep="")
> xsol[xsol>1e-4]
      x30       x75       x88 
0.1465990 0.5314489 0.3219522 

Conclusion 

QuadProg is a well-known QP solver for R. Unfortunately it has some drawbacks:

• It is not very suited for more complex models: the matrix oriented input format is just too cumbersome. 
• The quadratic coefficient matrix must be strictly positive definite.
• All variables are quadratic, linear variables are not available. 

For quick experimentation with these types of models, CVXR is likely a better environment. This exercise showed: 

• The models are more readable, assuming you are familiar with matrix notation.
• We can formulate models with linear variables without a problem.
• Models with a quadratic coefficient matrix that are borderline positive semi-definite are accepted by both CVXR and the solver I used: OSQP.

Two useful R tricks are discussed:

• nearPD(A) to find a positive-definite matrix close to A.
• Printing a sparse vector using names. 

References

1.  R minimum variance portfolio: solve non invertible matrix, https://stackoverflow.com/questions/64067591/r-minimum-variance-portfolio-solve-non-invertible-matrix
2. Covariance Matrix not positive definite in portfolio models, https://yetanothermathprogrammingconsultant.blogspot.com/2018/04/covariance-matrix-not-positive-definite.html
3. OSQP, Status values and errors, https://osqp.org/docs/interfaces/status_values.html

CVXPY large memory allocation

In [1] a simple regression problem was stated and solved with CVXPY. The number of observations is very large (\(200,000\)), while the number of coefficients to estimate is moderate (\(100\)). The first formulation is simply:

Regression I
\[\begin{align}\min_{\color{darkred}w}\>& ||\color{darkblue}y-\color{darkblue}X\color{darkred}w||^2_2 \end{align}\]

In Python code, this can look like:

import cvxpy as cp
import numpy as np

N = 200000
M = 100

X = np.random.normal(0, 1, size=(N, M))
y = np.random.normal(0, 1, size=N)

w = cp.Variable(M)
prob = cp.Problem(cp.Minimize(cp.sum_squares(y - X @ w)))
prob.solve()
print("status:",prob.status)
print("obj:",prob.value)

Unfortunately, this is giving the following runtime error:

[ec2-user@ip-172-30-0-79 etc]$ python3 ls0.py 
Traceback (most recent call last):
  File "ls0.py", line 12, in <module>
    prob.solve()
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 289, in solve
    return solve_func(self, *args, **kwargs)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 567, in _solve
    self._construct_chains(solver=solver, gp=gp)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 510, in _construct_chains
    raise e
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 501, in _construct_chains
    self._intermediate_chain.apply(self)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/chain.py", line 65, in apply
    problem, inv = r.apply(problem)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/qp2quad_form/qp2symbolic_qp.py", line 60, in apply
    return super(Qp2SymbolicQp, self).apply(problem)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/canonicalization.py", line 58, in apply
    problem.objective)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/canonicalization.py", line 96, in canonicalize_tree
    canon_arg, c = self.canonicalize_tree(arg)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/canonicalization.py", line 99, in canonicalize_tree
    canon_expr, c = self.canonicalize_expr(expr, canon_args)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/canonicalization.py", line 108, in canonicalize_expr
    return self.canon_methods[type(expr)](expr, args)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/qp2quad_form/atom_canonicalizers/quad_over_lin_canon.py", line 29, in quad_over_lin_canon
    return SymbolicQuadForm(t, eye(affine_expr.size)/y, expr), [affine_expr == t]
  File "/home/ec2-user/.local/lib/python3.7/site-packages/numpy/lib/twodim_base.py", line 201, in eye
    m = zeros((N, M), dtype=dtype, order=order)
MemoryError: Unable to allocate array with shape (200000, 200000) and data type float64
[ec2-user@ip-172-30-0-79 etc]$ 

This is interesting: CVXPY seems to allocate a \(200,000\times 200,000\) matrix here. Actually it is an identity matrix! This can be seen from the name eye. This seems to be related to the reformulation into a QP model.

To get this a bit more under control, lets make the size a bit smaller and use a memory profiler. At least we get some information about the memory usage.

[ec2-user@ip-172-30-0-79 etc]$ python3 -m memory_profiler ls1.py 
-----------------------------------------------------------------
           OSQP v0.6.0  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2019
-----------------------------------------------------------------
problem:  variables n = 20100, constraints m = 20000
          nnz(P) + nnz(A) = 2040000
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1   0.0000e+00   4.38e+00   2.98e+04   1.00e-01   1.03e+00s
  50   1.9804e+04   2.13e-09   4.01e-07   1.00e-01   1.65e+00s
plsh   1.9804e+04   4.01e-15   1.03e-12   --------   2.45e+00s

status:               solved
solution polish:      successful
number of iterations: 50
optimal objective:    19804.0226
run time:             2.45e+00s
optimal rho estimate: 1.35e-02

status: optimal
obj: 19804.022648294507
Filename: ls1.py

Line #    Mem usage    Increment   Line Contents
================================================
      4   50.188 MiB   50.188 MiB   @profile
     5                             def f():
     6                             #    N = 200000
     7   50.188 MiB    0.000 MiB       N = 20000
     8   50.188 MiB    0.000 MiB       M = 100
     9                             
    10   50.188 MiB    0.000 MiB       np.random.seed(123)
    11   65.477 MiB   15.289 MiB       X = np.random.normal(0, 1, size=(N, M))
    12   65.707 MiB    0.230 MiB       y = np.random.normal(0, 1, size=N)
    13                             
    14   65.707 MiB    0.000 MiB       w = cp.Variable(M)
    15   65.707 MiB    0.000 MiB       prob = cp.Problem(cp.Minimize(cp.sum_squares(y - X @ w)))
    16  412.086 MiB  346.379 MiB       prob.solve(verbose=True)
    17  412.086 MiB    0.000 MiB       print("status:",prob.status)
    18  412.086 MiB    0.000 MiB       print("obj:",prob.value)

We see indeed this is solved as QP model (it is solved by the QP solver OSQP), and in the solve statement we have this spike in memory usage.

The allocation is measured in MiBs or Mebibytes. A mebibyte is equal to \(2^{20}=1,048,576\) bytes.

It is probably a bad idea to allocate this identity matrix like this as a fully allocated matrix. The better approach would be not to use this matrix at all. A second best approach would be to make this identity matrix a sparse matrix.

Approach 1: use a conic programming solver

This is an easy fix. Just use a solver like ECOS:

[ec2-user@ip-172-30-0-79 etc]$ python3 -m memory_profiler ls1.py 

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -2.601e-05  +2e+04  5e-01  7e-06  1e+00  1e+04    ---    ---    1  2  - |  -  - 
 1  -9.312e-01  +2.050e-02  +3e+02  2e-02  1e-07  1e+00  2e+02  0.9829  1e-04   1  2  2 |  0  0
 2  -5.192e+00  +2.947e+00  +1e+02  7e-03  5e-08  8e+00  7e+01  0.6186  7e-02   2  3  3 |  0  0
 3  +6.518e+02  +8.866e+02  +3e+02  6e-02  4e-07  2e+02  1e+02  0.1800  9e-01   2  2  2 |  0  0
 4  +9.646e+02  +9.849e+02  +4e+00  1e-03  8e-09  2e+01  2e+00  0.9826  1e-04   7  7  6 |  0  0
 5  +1.130e+03  +1.144e+03  +7e-01  2e-03  1e-09  1e+01  4e-01  0.8836  2e-02   2  3  2 |  0  0
 6  +1.493e+03  +1.521e+03  +3e-01  1e-02  4e-10  3e+01  2e-01  0.9890  1e-01   1  1  1 |  0  0
 7  +2.057e+03  +2.080e+03  +5e-02  8e-03  9e-11  2e+01  3e-02  0.9271  5e-02   1  2  1 |  0  0
 8  +2.134e+03  +2.152e+03  +4e-02  4e-03  4e-11  2e+01  2e-02  0.8136  2e-01   2  2  2 |  0  0
 9  +2.244e+03  +2.256e+03  +3e-02  4e-02  2e-11  1e+01  1e-02  0.7107  1e-01   1  1  1 |  0  0
10  +1.995e+03  +2.004e+03  +4e-03  6e-02  2e-11  9e+00  4e-03  0.0005  8e-01   1  1  1 |  0  0
11  +6.609e+02  +6.615e+02  +8e+00  8e-03  2e-12  6e-01  4e+00  0.9890  8e-01   0  0  0 |  0  0
12  +6.293e+02  +6.423e+02  +1e+01  5e-03  2e-12  1e+01  8e+00  0.5325  3e-01   1  1  1 |  0  0
13  +5.873e+02  +7.017e+02  +3e+01  2e-02  3e-12  1e+02  2e+01  0.4675  5e-01   2  3  3 |  0  0
14  +7.549e+02  +1.687e+03  +4e+00  7e-03  4e-12  9e+02  3e+00  0.9890  1e-01   1  1  1 |  0  0
15  +1.292e+03  +1.300e+03  +2e+00  2e-03  3e-12  8e+00  9e-01  0.8813  2e-01   3  2  2 |  0  0
16  +2.141e+03  +2.149e+03  +2e-01  2e-02  2e-12  8e+00  9e-02  0.9890  5e-03   1  1  1 |  0  0
17  +2.819e+03  +2.827e+03  +2e-02  4e-02  4e-12  8e+00  1e-02  0.9052  3e-02   1  1  1 |  0  0
18  +2.757e+03  +2.774e+03  +5e-02  2e-02  2e-12  2e+01  2e-02  0.9746  4e-01   1  2  2 |  0  0
19  +2.615e+03  +2.631e+03  +1e-02  8e-02  2e-12  2e+01  8e-03  0.1350  2e-01   1  1  1 |  0  0
20  +2.002e+03  +2.009e+03  +2e-03  6e-02  3e-13  7e+00  3e-03  0.0029  9e-01   1  1  1 |  0  0
21  +2.884e+03  +2.959e+03  +2e-03  8e-02  3e-12  8e+01  9e-03  0.7080  7e-01   0  0  0 |  0  0
22  +2.861e+03  +2.934e+03  +2e-03  9e-02  3e-12  7e+01  9e-03  0.0007  9e-01   1  1  1 |  0  0
23  +1.847e+03  +1.857e+03  +3e-02  5e-02  2e-13  1e+01  2e-02  0.0026  1e+00   0  0  0 |  0  0
24  +6.545e+02  +6.548e+02  +1e+01  9e-03  2e-12  3e-01  6e+00  0.9890  6e-01   0  0  0 |  0  0
25  +7.003e+02  +7.142e+02  +2e+01  4e-03  2e-12  1e+01  1e+01  0.6078  2e-01   1  1  0 |  0  0
26  +6.232e+02  +7.771e+02  +5e+01  2e-02  3e-12  2e+02  3e+01  0.5706  4e-01   2  3  3 |  0  0
27  +5.916e+02  +1.396e+03  +8e+00  3e-03  2e-12  8e+02  6e+00  0.9890  1e-01   1  1  1 |  0  0
28  +1.468e+03  +1.658e+03  +1e+00  2e-03  3e-12  2e+02  7e-01  0.9890  9e-02   1  1  1 |  0  0
29  +2.046e+03  +2.089e+03  +4e-01  1e-02  4e-13  4e+01  2e-01  0.9890  3e-02   1  1  1 |  0  0
30  +1.720e+03  +1.744e+03  +7e-02  5e-02  3e-14  2e+01  6e-02  0.1873  2e-01   1  1  1 |  0  0
31  +6.708e+02  +6.722e+02  +2e+01  7e-03  2e-12  1e+00  1e+01  0.9890  6e-01   0  0  0 |  0  0
32  +5.985e+02  +6.158e+02  +5e+01  3e-03  9e-13  2e+01  3e+01  0.9260  3e-01   1  1  1 |  0  0
33  +7.045e+02  +9.047e+02  +1e+02  2e-02  2e-12  2e+02  6e+01  0.6591  5e-01   1  2  3 |  0  0
34  +4.246e+02  +8.139e+02  +3e+01  2e-03  1e-13  4e+02  3e+01  0.9890  2e-01   1  1  1 |  0  0
35  +1.140e+03  +1.436e+03  +7e-01  3e-03  7e-13  3e+02  3e+00  0.9600  1e-02   1  1  1 |  0  0
36  +1.597e+03  +1.707e+03  +1e+00  4e-03  3e-14  1e+02  9e-01  0.9890  2e-01   2  2  2 |  0  0
37  +2.184e+03  +2.226e+03  +1e-01  5e-03  1e-12  4e+01  1e-01  0.9890  6e-02   1  2  2 |  0  0
38  +2.540e+03  +2.556e+03  +2e-01  4e-03  2e-12  2e+01  1e-01  0.9890  2e-01   1  1  2 |  0  0
39  +3.055e+03  +3.065e+03  +5e-02  6e-03  2e-12  1e+01  3e-02  0.9890  8e-02   1  2  2 |  0  0
40  +3.450e+03  +3.458e+03  +4e-02  5e-03  2e-12  8e+00  2e-02  0.9890  1e-01   1  1  2 |  0  0
41  +3.932e+03  +3.939e+03  +1e-02  7e-03  2e-12  6e+00  8e-03  0.9890  1e-01   1  1  2 |  0  0
42  +4.347e+03  +4.352e+03  +1e-02  6e-03  2e-12  6e+00  5e-03  0.9890  1e-01   1  1  2 |  0  0
43  +4.796e+03  +4.801e+03  +5e-03  8e-03  2e-12  5e+00  3e-03  0.9890  1e-01   1  1  2 |  0  0
44  +5.206e+03  +5.211e+03  +3e-03  7e-03  2e-12  4e+00  2e-03  0.9890  1e-01   1  1  2 |  0  0
45  +5.633e+03  +5.637e+03  +2e-03  1e-02  2e-12  4e+00  1e-03  0.9890  1e-01   1  1  2 |  0  0
46  +6.028e+03  +6.032e+03  +1e-03  9e-03  2e-12  4e+00  8e-04  0.9890  1e-01   1  1  2 |  0  0
47  +6.430e+03  +6.433e+03  +9e-04  1e-02  2e-12  3e+00  5e-04  0.9890  1e-01   1  1  2 |  0  0
48  +6.811e+03  +6.814e+03  +7e-04  1e-02  2e-12  3e+00  4e-04  0.9890  1e-01   1  1  2 |  0  0
49  +7.188e+03  +7.191e+03  +5e-04  1e-02  2e-12  3e+00  3e-04  0.9890  1e-01   1  1  2 |  0  0
50  +7.551e+03  +7.553e+03  +3e-04  1e-02  2e-12  2e+00  2e-04  0.9890  1e-01   1  1  2 |  0  0
51  +7.908e+03  +7.910e+03  +2e-04  2e-02  2e-12  2e+00  1e-04  0.9890  2e-01   1  1  2 |  0  0
52  +8.251e+03  +8.253e+03  +2e-04  2e-02  2e-12  2e+00  1e-04  0.9890  2e-01   1  1  2 |  0  0
53  +8.586e+03  +8.588e+03  +1e-04  2e-02  2e-12  2e+00  8e-05  0.9890  2e-01   1  1  2 |  0  0
54  +8.911e+03  +8.913e+03  +1e-04  2e-02  2e-12  2e+00  6e-05  0.9890  2e-01   1  1  2 |  0  0
55  +9.228e+03  +9.230e+03  +8e-05  2e-02  2e-12  2e+00  4e-05  0.9890  2e-01   1  1  2 |  0  0
56  +9.534e+03  +9.535e+03  +6e-05  2e-02  2e-12  1e+00  4e-05  0.9890  2e-01   1  1  2 |  0  0
57  +9.779e+03  +9.780e+03  +5e-05  1e-01  2e-12  1e+00  3e-05  0.8357  2e-01   1  1  1 |  0  0
58  +9.766e+03  +9.768e+03  +2e-05  3e-01  2e-12  1e+00  1e-05  0.0001  9e-01   1  1  1 |  0  0
59  +3.518e+03  +3.517e+03  -6e+00  7e+00  9e-13  -5e-01  -3e+00  0.0000  1e+00   1  1  1 |  0  0
Unreliable search direction detected, recovering best iterate (58) and stopping.

Close to PRIMAL INFEASIBLE (within feastol=5.7e-06).
Runtime: 20.262543 seconds.

status: infeasible_inaccurate
obj: None
Filename: ls1.py

Line #    Mem usage    Increment   Line Contents
================================================
     4   50.312 MiB   50.312 MiB   @profile
     5                             def f():
     6                             #    N = 200000
     7   50.312 MiB    0.000 MiB       N = 20000
     8   50.312 MiB    0.000 MiB       M = 100
     9                             
    10   50.312 MiB    0.000 MiB       np.random.seed(123)
    11   65.578 MiB   15.266 MiB       X = np.random.normal(0, 1, size=(N, M))
    12   65.836 MiB    0.258 MiB       y = np.random.normal(0, 1, size=N)
    13                             
    14   65.836 MiB    0.000 MiB       w = cp.Variable(M)
    15   65.836 MiB    0.000 MiB       prob = cp.Problem(cp.Minimize(cp.sum_squares(y - X @ w)))
    16   76.824 MiB   10.988 MiB       prob.solve(solver=cp.ECOS,verbose=True)
    17   76.824 MiB    0.000 MiB       print("status:",prob.status)
    18   76.824 MiB    0.000 MiB       print("obj:",prob.value)

The results are mixed. We certainly don't see this crazy memory allocation anymore. It is now a very small 10 MiB. However the solver is not numerically stable enough to solve this problem.

There is another solver that comes with CVXPY that we can try:

[ec2-user@ip-172-30-0-79 etc]$ python3 -m memory_profiler ls1.py 
----------------------------------------------------------------------------
SCS v2.1.1 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 2000002
eps = 1.00e-04, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 0, rho_x = 1.00e-03
Variables n = 101, constraints m = 20002
Cones:soc vars: 20002, soc blks: 1
WARN: aa_init returned NULL, no acceleration applied.
Setup time: 7.44e-01s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 1.22e+18  7.87e+18  1.00e+00 -9.64e+18  4.27e+21  4.18e+21  1.92e-02 
   100| 8.68e+17  7.02e+17  5.09e-01  1.71e+21  5.25e+21  3.54e+21  6.67e-01 
   200| 8.04e+17  5.33e+17  4.23e-01  2.16e+21  5.33e+21  3.16e+21  1.31e+00 
   300| 7.37e+17  4.39e+17  3.72e-01  2.39e+21  5.22e+21  2.83e+21  1.96e+00 
   400| 6.73e+17  3.74e+17  3.35e-01  2.51e+21  5.04e+21  2.53e+21  2.59e+00 
   500| 6.14e+17  3.23e+17  3.05e-01  2.57e+21  4.83e+21  2.26e+21  3.22e+00 
   600| 5.60e+17  2.82e+17  2.80e-01  2.59e+21  4.60e+21  2.01e+21  3.85e+00 
   700| 5.11e+17  2.48e+17  2.57e-01  2.58e+21  4.37e+21  1.79e+21  4.50e+00 
   800| 4.67e+17  2.20e+17  2.36e-01  2.56e+21  4.14e+21  1.58e+21  5.15e+00 
   900| 4.27e+17  1.96e+17  2.17e-01  2.52e+21  3.92e+21  1.40e+21  5.79e+00 
  1000| 3.91e+17  1.74e+17  1.99e-01  2.47e+21  3.70e+21  1.23e+21  6.46e+00 
  1100| 3.58e+17  1.56e+17  1.81e-01  2.42e+21  3.49e+21  1.07e+21  7.09e+00 
  1200| 3.28e+17  1.40e+17  1.64e-01  2.36e+21  3.29e+21  9.29e+20  7.72e+00 
  1300| 3.01e+17  1.26e+17  1.48e-01  2.30e+21  3.10e+21  7.98e+20  8.36e+00 
  1400| 2.77e+17  1.14e+17  1.31e-01  2.24e+21  2.92e+21  6.78e+20  9.00e+00 
  1500| 2.55e+17  1.03e+17  1.15e-01  2.18e+21  2.75e+21  5.67e+20  9.65e+00 
  1600| 2.35e+17  9.31e+16  9.88e-02  2.12e+21  2.59e+21  4.65e+20  1.03e+01 
  1700| 2.17e+17  8.44e+16  8.25e-02  2.06e+21  2.43e+21  3.71e+20  1.09e+01 
  1800| 2.00e+17  7.67e+16  6.62e-02  2.01e+21  2.29e+21  2.84e+20  1.16e+01 
  1900| 1.85e+17  6.98e+16  4.97e-02  1.95e+21  2.15e+21  2.03e+20  1.23e+01 
  2000| 1.72e+17  6.37e+16  3.30e-02  1.89e+21  2.02e+21  1.29e+20  1.29e+01 
  2100| 1.59e+17  5.81e+16  1.60e-02  1.84e+21  1.90e+21  5.98e+19  1.36e+01 
  2200| 8.80e-03  1.43e-01  1.60e-05  1.20e+04  1.20e+04  1.94e-16  1.42e+01 
  2300| 6.45e-03  1.37e-01  1.62e-05  1.23e+04  1.23e+04  6.83e-17  1.49e+01 
  2400| 6.25e-03  1.34e-01  1.40e-05  1.26e+04  1.26e+04  2.15e-16  1.55e+01 
  2500| 6.06e-03  1.30e-01  1.21e-05  1.28e+04  1.28e+04  2.25e-16  1.62e+01 
  2600| 5.89e-03  1.27e-01  1.04e-05  1.30e+04  1.30e+04  7.78e-17  1.68e+01 
  2700| 5.71e-03  1.24e-01  8.98e-06  1.33e+04  1.33e+04  8.14e-17  1.75e+01 
  2800| 5.55e-03  1.21e-01  7.70e-06  1.35e+04  1.35e+04  8.46e-17  1.82e+01 
  2900| 5.39e-03  1.18e-01  6.56e-06  1.37e+04  1.37e+04  8.78e-17  1.88e+01 
  3000| 5.24e-03  1.15e-01  5.57e-06  1.39e+04  1.39e+04  9.04e-17  1.95e+01 
  3100| 5.09e-03  1.12e-01  4.68e-06  1.41e+04  1.41e+04  9.40e-17  2.01e+01 
  3200| 4.94e-03  1.09e-01  3.90e-06  1.42e+04  1.42e+04  2.91e-16  2.08e+01 
  3300| 4.80e-03  1.07e-01  3.21e-06  1.44e+04  1.44e+04  9.96e-17  2.14e+01 
  3400| 4.67e-03  1.04e-01  2.60e-06  1.46e+04  1.46e+04  1.03e-16  2.21e+01 
  3500| 4.54e-03  1.01e-01  2.05e-06  1.48e+04  1.48e+04  3.18e-16  2.27e+01 
  3600| 4.41e-03  9.87e-02  1.57e-06  1.49e+04  1.49e+04  1.09e-16  2.33e+01 
  3700| 4.29e-03  9.62e-02  1.14e-06  1.51e+04  1.51e+04  1.12e-16  2.40e+01 
  3800| 4.16e-03  9.37e-02  7.64e-07  1.52e+04  1.52e+04  1.14e-16  2.47e+01 
  3900| 4.05e-03  9.12e-02  4.28e-07  1.54e+04  1.54e+04  1.17e-16  2.53e+01 
  4000| 3.93e-03  8.88e-02  1.30e-07  1.55e+04  1.55e+04  1.20e-16  2.60e+01 
  4100| 3.82e-03  8.65e-02  1.34e-07  1.57e+04  1.57e+04  1.23e-16  2.66e+01 
  4200| 3.72e-03  8.42e-02  3.67e-07  1.58e+04  1.58e+04  1.25e-16  2.73e+01 
  4300| 3.61e-03  8.20e-02  5.72e-07  1.59e+04  1.59e+04  1.28e-16  2.80e+01 
  4400| 3.51e-03  7.98e-02  7.52e-07  1.60e+04  1.60e+04  3.92e-16  2.86e+01 
  4500| 3.41e-03  7.77e-02  9.11e-07  1.62e+04  1.62e+04  1.33e-16  2.93e+01 
  4600| 3.32e-03  7.56e-02  1.05e-06  1.63e+04  1.63e+04  1.36e-16  3.00e+01 
  4700| 3.22e-03  7.35e-02  1.17e-06  1.64e+04  1.64e+04  4.15e-16  3.07e+01 
  4800| 3.13e-03  7.15e-02  1.27e-06  1.65e+04  1.65e+04  1.41e-16  3.13e+01 
  4900| 3.04e-03  6.96e-02  1.36e-06  1.66e+04  1.66e+04  4.30e-16  3.20e+01 
  5000| 2.96e-03  6.77e-02  1.44e-06  1.67e+04  1.67e+04  4.37e-16  3.27e+01 
----------------------------------------------------------------------------
Status: Solved/Inaccurate
Hit max_iters, solution may be inaccurate
Timing: Solve time: 3.27e+01s
Lin-sys: nnz in L factor: 2025055, avg solve time: 5.82e-03s
Cones: avg projection time: 3.02e-05s
Acceleration: avg step time: 7.51e-07s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 0.0000e+00, dist(y, K*) = 3.6380e-12, s'y/|s||y| = -6.2180e-16
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 2.9568e-03
dual res:   |A'y + c|_2 / (1 + |c|_2) = 6.7689e-02
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 1.4408e-06
----------------------------------------------------------------------------
c'x = 16701.6612, -b'y = 16701.6131
============================================================================
status: optimal_inaccurate
obj: 16701.661202598487
Filename: ls1.py

Line #    Mem usage    Increment   Line Contents
================================================
     4   50.438 MiB   50.438 MiB   @profile
     5                             def f():
     6                             #    N = 200000
     7   50.438 MiB    0.000 MiB       N = 20000
     8   50.438 MiB    0.000 MiB       M = 100
     9                             
    10   50.438 MiB    0.000 MiB       np.random.seed(123)
    11   65.699 MiB   15.262 MiB       X = np.random.normal(0, 1, size=(N, M))
    12   65.957 MiB    0.258 MiB       y = np.random.normal(0, 1, size=N)
    13                             
    14   65.957 MiB    0.000 MiB       w = cp.Variable(M)
    15   65.957 MiB    0.000 MiB       prob = cp.Problem(cp.Minimize(cp.sum_squares(y - X @ w)))
    16   76.953 MiB   10.996 MiB       prob.solve(solver=cp.SCS,verbose=True)
    17   76.953 MiB    0.000 MiB       print("status:",prob.status)
    18   76.953 MiB    0.000 MiB       print("obj:",prob.value)

That solver also has problems. In this case we ran out of iterations. We can allow for more iterations by adding the argument max_iters=nnnn to the solve statement. SCS is using a first order algorithm. These algorithms typically use a very large number of iterations.

Approach 2: use a different formulation

We can use a different formulation:

Regression II
\[\begin{align}\min_{\color{darkred}w,\color{darkred}r}\>& ||\color{darkred}r||^2_2 \\ & \color{darkred}r = \color{darkblue}y - \color{darkblue}X\color{darkred}w \end{align}\]

Here we add a bunch of (free) variables \(r\) for the residuals, and also a bunch of linear constraints. The payback is that we have a much simpler quadratic objective. In addition it is 100% positive (semi-) definite. A reasonable implementation can look like:

import cvxpy as cp
import numpy as np

# N = 200000
N = 20000
M = 100

X = np.random.normal(0, 1, size=(N,M))
y = np.random.normal(0, 1, size=(N,1))

w = cp.Variable((M,1))
r = cp.Variable((N,1))
prob = cp.Problem(cp.Minimize(r.T @ r), [r == y - X @ w])
prob.solve(solver=cp.OSQP,verbose=True)
print("status:",prob.status)
print("obj:",prob.value)

This should work, but it doesn't. We see:

[ec2-user@ip-172-30-0-79 etc]$ python3 ls2a.py 
Traceback (most recent call last):
  File "ls2a.py", line 14, in <module>
    prob.solve(solver=cp.OSQP,verbose=True)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 289, in solve
    return solve_func(self, *args, **kwargs)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 567, in _solve
    self._construct_chains(solver=solver, gp=gp)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 510, in _construct_chains
    raise e
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/problems/problem.py", line 499, in _construct_chains
    construct_intermediate_chain(self, candidate_solvers, gp=gp)
  File "/home/ec2-user/.local/lib/python3.7/site-packages/cvxpy/reductions/solvers/intermediate_chain.py", line 70, in construct_intermediate_chain
    raise DCPError("Problem does not follow DCP rules. Specifically:\n" + append)
cvxpy.error.DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP. Its following subexpressions are not:
var1.T * var1

Bummer. This is a perfectly convex problem!

Some alternative formulations to generate a QP do not help:

• cp.sum_squares(r) yields a large allocation when forming a QP  
• sum([r[i]**2 for i in range(N)] is very slow and very memory hungry
• quad_form with an identity matrix is not efficient 

So I had no luck in being able to generate a QP problem for OSQP without very large memory allocations.

This did not work out so well. I am not really able to solve relatively large, but easy least squares problems using standard formulations.

Conclusion

• CVXPY thinks \(r^Tr = \sum_i r_i^2\) is not convex.
• We cannot always generate large problems for the OSQP QP solver due to large dense memory allocations in CVXPY. I probably would consider this a bug: it is never a good idea to physically create and allocate a large dense identity matrix in any somewhat serious code.
• CVXPY can generate large instances for the conic solvers ECOS and SCS. But they have some troubles solving this problem.
• A commercial solver like Cplex has no problem with this; 2 iterations, 0.2 seconds. 
• We can conclude that CVXPY (and its collection of standard solvers) is best suited for smaller and medium-sized problems. When tackling large problems you may hit some issues. Of course, many interesting and relevant problems are not that big and can be solved conveniently using CVXPY.

References

1. Least squares problem run out of memory, https://stackoverflow.com/questions/59315300/least-squares-problem-run-out-of-memory-in-cvxpy