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

Nonlinear variant of a knapsack problem

In [1] a problem is posed:

Original Problem
\[\begin{align}\max & \sum_i \log_{100}(\color{darkblue}w_i) \cdot \color{darkred}x_i \\ & \frac{\sum_i \color{darkblue}w_i \color{darkred}x_i}{\sqrt{\color{darkred}k}} \le 10,000\\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \color{darkred}x_i \in \{0,1\} \end{align}\]

The data vector \(w_i\) is assumed to integer valued with \(w_i\ge 1\). Hence the logarithm can be evaluated without a problem. Also we can assume that the number of items can be around 1,000.

I don't think I have ever seen \(\log_{100}()\) being used in a model. Most programming (and modeling) languages only support natural logarithms \(\ln()\) and may be \(\log_{10}()\). We can convert things by: \[\log_{100}(x) = \frac{\ln(x)}{\ln(100)}\] This means the objective can be written as \[\max \frac{1}{\ln(100)} \sum_i \ln(w_i) x_i\] Essentially, the \(\log_{100}()\) function just adds a scaling factor. We can simplify the objective to \[\max \sum_i \ln(w_i)x_i\] (The objective value will be different, but the optimal solution will be the same). Note that this is a linear objective.

log100(x) vs ln(x)

In this post I'll discuss a few different formulations:

• As as MINLP (Mixed Integer Nonlinear Programming) model. This is an obvious choice considering the square root in the model.
• An MIQCP (Mixed Integer Quadratically Constrained Programming) model. This (convex) quadratic reformulation allows us to use alternative solvers.
• A linear MIP model. This formulation may not be immediately obvious. Of course with this model we can use widely available MIP solvers.

MINLP Formulation

The constraint can be rewritten as: \[\sum_i w_i x_i \le 10,000 \sqrt{k}\] If we ignore the all zero solution, we can assume \(k\ge 1\). This bound will make sure the square root function is always differentiable. With this, we have a standard MINLP model.

MINLP Model
\[\begin{align}\max & \sum_i \ln(\color{darkblue}w_i) \color{darkred}x_i \\ & \sum_i \color{darkblue}w_i \color{darkred}x_i \le 10,000\sqrt{\color{darkred}k}\\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \color{darkred}x_i \in \{0,1\} \\ & \color{darkred} k \ge 1 \end{align}\]

This model has only a single, well behaved, non-linearity. Literally, there is only one nonlinear nonzero element.  In addition, the model is convex. So we don't expect any problems. For \(w_i\),  I generated 1000 random integer values from the interval \([1,10000]\). All timings are in seconds.

MINLP Results
Solver	Obj	Time	Notes
Dicopt	1057.1355	0.6	3 NLP, 2 MIP subproblems
SBB	1056.9361	95	Node limit exceeded
Bonmin	1056.9472	600	Time limit exceeded
Bonmin	1057.1355	16	Option: bonmin.algorithm B-OA

The outer-approximation based algorithms (Dicopt, Bonmin with B-OA option) do much better than the branch & bound algorithms (SBB, default Bonmin). Even the global solvers do better:

Global Solver Results
Solver	Obj	Time
Baron	1057.1355	2
Antigone	1057.1355	1
Couenne	1057.1355	9

MIQCP Formulation

The problem can also be formulated as a convex MIQCP model:

MIQCP Model
\[\begin{align}\max & \sum_i \ln(\color{darkblue}w_i) \color{darkred}x_i \\ & \color{darkred}y^2 \le \color{darkred}k\\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \color{darkred}y = \frac{\sum_i \color{darkblue}w_i \color{darkred}x_i}{10,000} \\ & \color{darkred}x_i \in \{0,1\} \\ & \color{darkred} k \ge 1, \color{darkred}y \ge 0 \end{align}\]

Cplex has two ways to form relaxations:

• MIQCPstrat 1: SOCP-based branch & bound
• MIQCPstrat 2: OA-based branch & cut

The results are shown in the next section. We see again that the Outer Approximation approach works better.

MIP Formulation

Finally, we can also linearize this model by observing that \(k\in \{1,\dots,1000\}\). So we don't really have a continuous function \(f(k)=\sqrt{k}\), but rather only need function values at the integer points. We can exploit this by making this explicit. We can write: \[\begin{align} & k = \sum_i i\cdot \delta_i \\ & \sqrt{k} = \sum_i \sqrt{i}\cdot \delta_i \\ & \sum_i \delta_i = 1 \\ & \delta_i \in \{0,1\}\end{align}\] This is essentially a SOS1 (Special Ordered Set of Type 1) structure implementing a table lookup. The MIP model looks like:

MIP Model

\[\begin{align}\max & \sum_i \ln(\color{darkblue}w_i) \color{darkred}x_i \\ & \color{darkred} k = \sum_i \color{darkred}x_i \\ & \sum_i \color{darkblue}w_i \color{darkred} x_i \le 10,000 \color{darkred}q \\ & \color{darkred} k = \sum_i i \cdot \color{darkred}\delta_i \\ & \color{darkred} q = \sum_i \sqrt{i} \cdot \color{darkred}\delta_i \\ & \sum_i \color{darkred} \delta_i = 1 \\ & \color{darkred}x_i, \color{darkred}\delta_i \in \{0,1\} \\ & \color{darkred} k, \color{darkred} q \ge 1 \end{align}\]

When we solve this problem we see:

MIQCP and MIP Results
Model	Solver	Obj	Time	Notes
MIQCP	Cplex	1057.1355	600	SOCP-based branch & bound, time limit
MIQCP	Cplex	1057.1355	0.3	OA-based branch & cut
MIP	Cplex	1057.1355	0.6

Conclusion: the question in the original post was: how to solve this problem? Here we proposed three different models: an MINLP, MIQCP and MIP model. All these models can be solved quickly. It is noted that the quadratic and linear models are not approximations: they give the same solution as the original MINLP model. Pure non-linear branch & bound methods are having a bit of a problem with the MINLP model, but Outer-Approximation works very well.

References

1. How do we solve a variant of the knapsack problem in which the capacity of the knapsack keeps increasing as we add more items into the knapsack?, https://stackoverflow.com/questions/59242370/how-do-we-solve-a-variant-of-the-knapsack-problem-in-which-the-capacity-of-the-k
2. Pierre Bonami, Andrea Tramontani, Recent improvement to MISOCP in CPLEX, Informs 2015, http://www.optimizationdirect.com/data/informs2015_bonami.pdf

Appendix: automatic reformulation with CVXPY

CVXPY accepts the MINLP model, and can automatically reformulate this into a MISOCP model:

x = cp.Variable(n,integer=True)
obj = cp.Maximize(np.log(w)*x)
con = w*x <= 10000 * cp.sqrt(cp.sum(x)) 
prob = cp.Problem(obj,[con, x >= 0, x <= 1])
prob.solve(solver=cp.ECOS_BB, max_iters=10000)
print("status:",prob.status)
print("Obj:",prob.value)

Unfortunately, the MISOCP solver that comes with CVXPY is not the most reliable one. It shows:

status: optimal_inaccurate
Obj: 1057.1085712222387
X: [-1.52640229e-10 -4.44776117e-11 -6.08671797e-11 -9.71628737e-11
 -9.95050649e-11 -1.23087763e-10 -8.63183870e-11 -4.39955343e-11
  1.00000000e+00 -6.54372650e-11 -3.95051292e-11 -5.86051991e-11
...

The reported objective is close to what we have found earlier.

Elementwise vs matrix multiplication

Introduction

In this post I want to demonstrate some interesting issues in how CVXR and CVXPY deal with elementwise vs matrix multiplication and compare that to the host languages (R and Python, respectively). CVXR and CVXPY deviate a bit from what their host languages do. When implementing optimization models with these tools, it is important to be aware of these details.

Elementwise and matrix multiplication 

There are two often used methods to perform the multiplication of matrices (and vectors). The first is simply elementwise multiplication: \[c_{i,j} = a_{i,j} \cdot b_{i,j}\] In mathematics, this is sometimes referred to as the Hadamard product. The notation is typically: \[C = A \circ B\] but sometimes we see: \[C = A \odot B\] This product only works when \(A\) and \(B\) have the same shape. I.e. \[\begin{matrix} C&=&A&\circ&B\\ (m \times n) &&(m \times n)&&(m \times n)\end{matrix}\]

The standard matrix multiplication \(C=A\cdot B\) or \[c_{i,j} = \sum_k a_{i,k} b_{k,j}\] has a different rule for conformance:  \[\begin{matrix} C&=&A&\cdot&B\\ (m \times n) &&(m \times k)&&(k \times n)\end{matrix}\] Most of the time, the dot operator is dropped and we write \(C = A B\).

In optimization modeling both forms are used a lot.

R, CVXR

R has two operators for multiplication:

• * for elementwise multiplication
• %*% for matrix multiplication

A vector is considered as a column vector (i.e. an \(n\)-vector is like a \(n \times 1\) matrix).  There is one special thing in R, as shown here:

> m <- 2
> n <- 3
> A <- matrix(c(1,2,3,4,5,6),m,n)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> 
> x <- c(1,2,3)
>
 
#
# matrix multiplication
# 
> A %*% x
     [,1]
[1,]   22
[2,]   28
> 

#
# elementwise multiplication
#
> A * A
     [,1] [,2] [,3]
[1,]    1    9   25
[2,]    4   16   36

#
# but this also works
#
> A * x
     [,1] [,2] [,3]
[1,]    1    9   10
[2,]    4    4   18
>

The last multiplication is surprising as \(A\) is a \(2 \times 3\) matrix and \(x\) is different in shape. Well, R may extend and recycle vectors to make them as large as needed. In this case \(x\) is duplicated and then considered as a  \(2 \times 3\) matrix. More or less like:

> x
[1] 1 2 3
> matrix(x,2,3)
     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    2    1    3

The modeling tool CVXR follows R notation and implements both * and %*% (for elementwise and matrix multiplication). However, CVXR is not implementing the extending and recycling of vectors that are too small.

Concept: recycling

Just to emphasize the concept of recycling. If an operation requires two vectors of the same length, R may make the shorter vector longer by recycling (duplicating). Here is an example:

> a <- 1:10
> b <- 1:2
> c <- a + b
> c
 [1]  2  4  4  6  6  8  8 10 10 12

In this example the vector a has elements 1 through 10. The vector b is too short, so it is recycled. When added to a, b is functionally equal to rep(c(1,2),5). When multiples of b do not exactly fit a, we effectively have a fractional duplication number. E.g. when we use b <- 1:3, we get a message:
Warning message:
In a + b : longer object length is not a multiple of shorter object length

This recycling trick is somewhat unique to R (I don't know about other languages doing this).

Example

In [2] the product: \[b_{i,j} = a_{i,j} \cdot x_i \] is implemented in R with the elementwise product:

> m <- 2
> n <- 3
> A <- matrix(c(1,2,3,4,5,6),m,n)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> x <- c(1,2)
> B <- A*x
> B
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12

This again uses recycling of \(x\). CVXR is not doing this automatically. We can see this here:

> library(CVXR)
> x <- Variable(m)
> B <- Variable(m,n)
> e <- rep(1,n)
> problem <- Problem(Minimize(0),
+                    list(x == c(1,2), 
+                         B == A * x ))
Error in sum_shapes(lapply(object@args, function(arg) { : 
  Incompatible dimensions

Here \(x\) is now a CVXR variable. As the vector \(x\) is not recycled, we end up with two different shapes and elementwise multiplication is refused. So how do we do something like this in CVXR?

The recycling operation:

> x
[1] 1 2
> matrix(x,m,n)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2

can be expressed in matrix notation as: \[x \cdot e^T\] where \(e\) is a (column) vector of ones. This is sometimes called an outerproduct. I.e. we can write our assignment as \[B = A \circ (x \cdot e^T)\] In a CVXR model this can look like:

> library(CVXR)
> x <- Variable(m)
> B <- Variable(m,n)
> e <- rep(1,n)
> problem <- Problem(Minimize(0),
+                    list(x == c(1,2), 
+                         B == A * (x %*% t(e))))
> sol <- solve(problem)
> sol$status
[1] "optimal"
> sol$getValue(x)
     [,1]
[1,]    1
[2,]    2
> sol$getValue(B)
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12

Note: the constraint B == A * (x %*% t(e)) has both a matrix multiplication and a elementwise multiplication. This is rather funky.

Conclusion: if your matrix operations rely on recycling, you will need to rework things a bit to have this work correctly in CVXR. CVXR does not do recycling.

Python, CVXPY

In the previous section we saw that there are subtle differences between R's and CVXR's elementwise multiplication semantics.

Let's now look at Python and CVXPY.

Since Python 3.5 we have two multiplication operators:

• * for elementwise multiplication 
• @ for matrix multiplication 

CVXPY has different rules:

• *, @  and matmul for matrix multiplication
• multiply for elementwise multiplication

Example

import numpy as np
import cvxpy as cp

#
# In Python/numpy * indicates elementwise multiplication
#
A = np.array([[1,2],[3,4]])
B = np.array([[1,1],[2,2]])
C = A*B
print(A)
# output:
# [[ 1 2]
#  [ 3 4]]
print(C)
# output:
# [[ 1 2]
#  [ 6 8]]

#
# In CVXPY * indicates matrix multiplication
#
A = cp.Variable((2,2))
C = cp.Variable((2,2))
prob = cp.Problem(cp.Minimize(0),
                  [A == [[1,2],[3,4]],
                   C == A*B])
prob.solve(verbose=True)
print(A.value)
# output:
# [[1. 3.]
#  [2. 4.]]
print(C.value)
# [[ 7.  7.]
#  [10. 10.]]

Here we see some differences between Python/Numpy and CVXPY. First the interpretation of Python lists (the values for A) is different. And secondly: the semantics of * are different. This may cause some confusion.

References

1. CVXR, Convex Optimization in R, https://cvxr.rbind.io/
2. CVXR Elementwise multiplication of matrix with vector, https://stackoverflow.com/questions/59224555/cvxr-elementwise-multiplication-of-matrix-with-vector

Opt Art

New book by TSP and Domino art creator Robert Bosch [1].

Content:

1. Optimization and the Visual Arts?
2. Truchet Tiles
3. Linear Optimization and the Lego Problem
4. The Linear Assignment problem and Cartoon Mosaics
5. Domino Mosaics
6. From the TSP to Continuous Line Drawings
7. TSP Art with Side Constraints
8. Knight's Tours
9. Labyrinth Design with Tiling and Pattern Matching
10. Mosaics with Side Constraints
11. Game-of-life Mosaics

Yes, it contains color pictures.

TSP Art

An example of a TSP drawing (from [2]):

25,000-city TSPortrait of George Dantzig [2,3]

Original fotograph

I guess more cities are needed to prevent his teeth from disappearing.

Domino Art

You can submit your picture and create a domino mosaic through NEOS [4].  I tried this with the Dantzig picture. With a "high quality" output image this creates a large model:

Reduced MIP has 13915 rows, 1511565 columns, and 4534695 nonzeros.
Reduced MIP has 1511565 binaries, 0 generals, 0 SOSs, and 0 indicators.

This is not small at all. Although the model is large, it is easy to solve. It only took 11 nodes (most of the work was in the root node). The model is an assignment problem with side constraints [5]. Basically the costs are the difference in gray scale value between the picture and the domino piece for each cell. The constraints are: exactly one domino in each cell and place all dominos. NEOS solves this by processing the picture in Matlab and then solving the model with GAMS/CPLEX. The output PNG file is also notably large: 2.8 MB. Too large to include in this blog. A screen grab looks like:

Screen capture of hi-res image

After zooming in we see indeed the dominos:

Zooming in

With some effort we can recognize part of the frame of the eye glasses here.

References

1. Robert Bosch, Opt Art, From Mathematical Optimization to Visual Design, 2019,  https://press.princeton.edu/books/hardcover/9780691164069/opt-art
2. http://www2.oberlin.edu/math/faculty/bosch/tspart-page.html
3. https://en.wikipedia.org/wiki/George_Dantzig
4. https://neos-guide.org/content/domino-art
5. Robert Bosch, Opt Art, Math Horizons 14(3) (2006) 6-9, https://www.maa.org/sites/default/files/pdf/upload_library/22/Evans/Horizons-Feb06-p.06-09.pdf