0

I am using the CVXR modelling package to solve a convex optimization problem. I know for sure that the problem is convex and that it follows the DCP rules, but if I check the DCP rules using CVXR it returns False. However, if I take the exact same problem and check it using CVXPY it returns True (as expected)

What is happening here? I attach a minimal reproducible example of this behavior in R and Python:

R code using CVXR

library(splines2)
library(CVXR)
deriv_basis = splines2::dbs(seq(0, 1, length.out=100), degree=3, intercept=T, df=30, derivs=2)
R = t(deriv_basis) %*% deriv_basis
beta_var = CVXR::::Variable(nrow(R))
q = CVXR::quad_form(beta_var, R)
CVXR::is_dcp(q)

[1] FALSE

write.table(x=R, file='R.csv'), row.names=F, sep=';')

Python code using CVXPY

import cvxpy
import pandas as pd

R = pd.read_csv('R.csv', sep=';').values
beta_var = cvxpy.Variable(R.shape[1])
q = cvxpy.quad_form(beta_var, R)
q.is_dcp()

Out[1]: True

Can someone explain what is happening here and how to solve it so I can use CVXR?

Álvaro Méndez Civieta
  • 725
  • 3
  • 12
  • 33

2 Answers2

2

The problem is the negative eigenvalue in the R matrix. If you fix that by setting it to zero, say, then it satisfies the dcp condition. I have also fixed the syntax errors in the code in the question and removed the redundant :: . Another possibility (not shown) is to use nearest_spd in the pracma package to adjust the R matrix.

library(splines2)
library(CVXR)

deriv_basis <- dbs(seq(0, 1, length.out=100), degree = 3, 
  intercept = TRUE, df = 30, derivs = 2)
R <- t(deriv_basis) %*% deriv_basis
e <- eigen(R)

# check decomposition
all.equal(R, e$vectors %*% diag(e$values) %*% t(e$vectors), 
 check.attributes = FALSE)
## [1] TRUE

e$values  # note negative value
##  [1]  1.095213e+08  1.095213e+08  1.056490e+07  1.055430e+07  1.052481e+07
##  [6]  1.046063e+07  1.034247e+07  1.015017e+07  9.866358e+06  9.485145e+06
## [11]  8.643220e+06  8.280963e+06  7.549803e+06  6.731472e+06  5.853402e+06
## [16]  4.949804e+06  4.056714e+06  3.209045e+06  2.437320e+06  1.759963e+06
## [21]  1.214976e+06  7.785251e+05  4.590441e+05  2.428199e+05  1.107300e+05
## [26]  4.060476e+04  1.040537e+04  1.320942e+03  7.239578e-09 -5.019224e-09

# zap negative eigenvalues making them zero
R <- with(e, vectors %*% diag(pmax(values, 0)) %*% t(vectors))

beta_var <- Variable(nrow(R))
q <- quad_form(beta_var, R)
is_dcp(q)
## [1] TRUE
G. Grothendieck
  • 254,981
  • 17
  • 203
  • 341
  • Thank you for the answer! And do you know why is there a difference in the DCP rules between R and python? – Álvaro Méndez Civieta Jun 08 '21 at 07:35
  • 2
    It can't be assumed that the DCP rules behave differently. It is conceivable that the Python software produces a PSD R matrix. Check out the eigenvalues on the Python side as we did in the answer on the R side. – G. Grothendieck Jun 08 '21 at 08:58
0

cvxpy should give the same results as CVXR as implied by G. Grothendieck since the DCP rules are the same. Something seems to have broken in recent versions of cvxpy. I've opened an issue on cvxpy github.

Balasubramanian Narasimhan
  • 49
  • 3