Friday, December 20, 2019

A facility location problem

The following problem is from [1]:

I am trying to solve an uncapacitated facility location problem in which I would like to find the location to open at most \(k\) facilities, and assign \(n\) clients to those facilities, each with cost \(c_i\) (proportionate to the distance of the client from the facility it's assigned to) for assignment with a budget \(B\). So I would like to find the locations of the facilities, and assign clients to those facilities such that the maximum number of clients are serviced within the budget \(B\).  I haven't found a variant like this in which there is an overall budget; I've only seen versions where we want to minimize the overall cost while assigning all the clients to facilities.

Let's give this a shot.

Example Data


A good way to familiarize yourself with a problem, is to invent some example data. Here is some random test data I used for the models below:


----     35 SET i  clients

client1 ,    client2 ,    client3 ,    client4 ,    client5 ,    client6 ,    client7 ,    client8 ,    client9 
client10,    client11,    client12,    client13,    client14,    client15,    client16,    client17,    client18
client19,    client20,    client21,    client22,    client23,    client24,    client25


----     35 SET j  facilities

facility1,    facility2,    facility3,    facility4


----     35 SET k  xy-coordinates

x,    y


----     35 PARAMETER U                    =      100.000  box size for generating random data
            PARAMETER maxdist              =      141.421  maximum distance (big-M)
            PARAMETER budget               =      500.000  available budget

----     35 PARAMETER client  location and cost data

                   x           y        cost

client1       17.175      84.327       6.611
client2       55.038      30.114       7.558
client3       29.221      22.405       6.274
client4       34.983      85.627       2.839
client5        6.711      50.021       0.864
client6       99.812      57.873       1.025
client7       99.113      76.225       6.413
client8       13.069      63.972       5.453
client9       15.952      25.008       0.315
client10      66.893      43.536       7.924
client11      35.970      35.144       0.728
client12      13.149      15.010       1.757
client13      58.911      83.089       5.256
client14      23.082      66.573       7.502
client15      77.586      30.366       1.781
client16      11.049      50.238       0.341
client17      16.017      87.246       5.851
client18      26.511      28.581       6.212
client19      59.396      72.272       3.894
client20      62.825      46.380       3.587
client21      41.331      11.770       2.430
client22      31.421       4.655       2.464
client23      33.855      18.210       1.305
client24      64.573      56.075       9.334
client25      76.996      29.781       3.799

The clients are uniformly distributed over the \([0,100]\times[0,100]\) square.

Client locations

1. High-level MINLP Formulation


The first step is to translate the "word-problem" into a high-level Mixed-Integer Nonlinear Programming (MINLP) model. The goal is not to achieve the best performance possible, but rather to have a "reference model": a precise mathematical description of the problem.

We use the following naming:


Names
Indices
\(i\)Clients
\(j\)Facilities
\(k\)x or y coordinate
Data
\(\color{darkblue}{\mathit{Cost}}_i\)Unit cost
\(\color{darkblue}{\mathit{Client}}_{i,k}\)Client locations
\(\color{darkblue}{\mathit{Budget}}\)Available budget
\(\color{darkblue}U\)Size of box containing client locations
Variables
\(\color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\} \)Assignment of client to facility
\(\color{darkred}{\mathit{ClientUsed}}_{i} \in \{0,1\} \)Client is assigned to a facility
\(\color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\)Facility locations
\(\color{darkred}{\mathit{Dist}}_{i,j} \in [0,\sqrt{2}\color{darkblue}U] \)Distance between client and facility
Zero if not assigned

With this we can start developing our model:


MINLP Problem
\[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_j \color{darkred}{\mathit{Assign}}_{i,j} \\ &\color{darkred}{\mathit{Dist}}_{i,j} = \color{darkred}{\mathit{Assign}}_{i,j} \cdot \sqrt{\sum_k \left(\color{darkblue}{\mathit{Client}}_{i,k}-\color{darkred}{\mathit{Facility}}_{j,k} \right)^2} \\ & \sum_{i,j} \color{darkblue}{\mathit{Cost}}_i \cdot \color{darkred}{\mathit{Dist}}_{i,j} \le \color{darkblue}{\mathit{Budget}}\\ &\color{darkred}{\mathit{ClientUsed}}_{i}, \color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\} \\ & \color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\\ & \color{darkred}{\mathit{Dist}}_{i,j} \in [0,\sqrt{2}\color{darkblue}U]\end{align}\]


We don't need to add an assignment constraint like\[\sum_j {\mathit{Assign}}_{i,j}  \le 1 \>\>\forall i \] as this is already captured in the constraint \[\mathit{ClientUsed}_i = \sum_j {\mathit{Assign}}_{i,j}\]

The distance constraint can be interpreted as \[{\mathit{Dist}}_{i,j} = \begin{cases} \sqrt{\displaystyle\sum_k \left({\mathit{Client}}_{i,k}-{\mathit{Facility}}_{j,k} \right)^2} & \text{if $\mathit{Assign}_{i,j}=1$} \\ 0 & \text{if $\mathit{Assign}_{i,j}=0$}\end{cases}\]

We can even use a global MINLP solver to get optimal solutions for smaller data sets (or "good" solutions for slightly larger problems). This will give us more confidence that the model is correctly representing the problem.

2. Mixed Integer Quadratically Constrained Model


The second step is to get rid of some of the non-linearities so that a quadratic model remains. We need to work on the distance constraint. By making it a big-M constraint, we can make everything quadratic. We get rid of the square root and the multiplication by the Assign variable.


MIQCP Problem
\[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_j \color{darkred}{\mathit{Assign}}_{i,j} \\ &\color{darkred}{\mathit{Dist}}_{i,j}^2 \ge \sum_k \left(\color{darkblue}{\mathit{Client}}_{i,k}-\color{darkred}{\mathit{Facility}}_{j,k} \right)^2 - \color{darkblue}M (1-\color{darkred}{\mathit{Assign}}_{i,j}) \\ & \sum_{i,j} \color{darkblue}{\mathit{Cost}}_i \cdot \color{darkred}{\mathit{Dist}}_{i,j} \le \color{darkblue}{\mathit{Budget}} \\  &\color{darkred}{\mathit{ClientUsed}}_{i}, \color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\}\\ & \color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\\ & \color{darkred}{\mathit{Dist}}_{i,j} \in [0,\sqrt{2}\color{darkblue}U]\\&\color{darkblue} M = 2\color{darkblue}U^2 \end{align}\]

In this model we have only a bound on the variable \(\mathit{Dist}_{i,j}\). That means it can be a bit larger than needed, but only as long as it obeys the Budget constraint. If we want to report on the actual distances, we need to recalculate them.

If we run this with Cplex we see:

CPLEX Error  5002: 'QCP_row_for_distance2(client1.facility1)' is not convex.

The new version Gurobi 9 can solve this as a nonconvex MIQCP. Of course we can also solve this with a global MINLP solver.

The problem would be easier if the costs do not increase proportionally with the distance but with the squared distance. That would yield a simpler and convex MIQCP problem. [2]

3. Mixed Integer Second Order Cone Formulation


The constraint \[y^2 \ge \sum_i x_i^2 \>\>\text{where $y\ge 0$}\] is a Second Order Cone constraint. This type of constraint is supported by Cplex. But only if this is specified (almost) exactly like this. I.e. no additional linear terms in the constraints. The mechanics are not complicated, just a bit tedious. We need to split our distance constraint in parts to shoehorn it into a SOCP framework:

MISOCP Problem
\[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_j \color{darkred}{\mathit{Assign}}_{i,j} \\ & \color{darkred}{\mathit{Dist}}_{i,j} \ge \color{darkred}\Delta_{i,j} - \color{darkblue}M (1-\color{darkred}{\mathit{Assign}}_{i,j})  \\ &\color{darkred}{\mathit{\Delta}}_{i,j}^2 \ge  \sum_k  \color{darkred} D_{i,j,k}^2\\ & \color{darkred} D_{i,j,k} = \color{darkblue}{\mathit{Client}}_{i,k}-\color{darkred}{\mathit{Facility}}_{j,k}  \\  & \sum_{i,j} \color{darkblue}{\mathit{Cost}}_i  \cdot\color{darkred}{\mathit{Dist}}_{i,j} \le \color{darkblue}{\mathit{Budget}} \\ &\color{darkred}{\mathit{ClientUsed}}_{i}, \color{darkred}{\mathit{Assign}}_{i,j} \in \{0,1\}\\ & \color{darkred}{\mathit{Facility}}_{j,k} \in [0,\color{darkblue}U]\\ & \color{darkred}{\mathit{Dist}}_{i,j}, \color{darkred}\Delta_{i,j} \in [0,\sqrt{2}\color{darkblue}U]\\  & \color{darkred}D_{i,j,k} \>\text{free} \\ & \color{darkblue} M = \sqrt{2}\color{darkblue}U\end{align}\]


This is now a convex problem. It can be solved with MISOCP solvers like Cplex, Gurobi and Mosek.

This is just a reformulation of the previous MIQCP model. Solvers do not typically recognize that form as a second order cone problem. So often the task is left to the modeler to help the solver and produce a convex model.

There are two objections against this special way of writing SOCP constraints:

  1. We need to have fairly active knowledge of how SOCP constraints can look like. I.e. we need to remember a bag of reformulation tricks. If you don't, you may miss out on solver capabilities.
  2. Even if we know how to do this, the reformulation is making the model less clear and less intuitive. We need to add extra comments why additional, at first sight unnecessary, variables and constraints were added.

Some high level modeling systems can apply reformulations automatically [3]. I think that is a good idea in order to keep formulations as natural, and as close to the original problem as possible.

This is a small problem. But Cplex has no easy time in solving it. Interestingly, it found the optimal solution completely at the end.


 6939104 1303772       20.7395    24       19.0000       20.8047 2.53e+08    9.50%
 6971225 1302201       20.0000    20       19.0000       20.7915 2.54e+08    9.43%
 7005795 1298106    infeasible             19.0000       20.7779 2.55e+08    9.36%
 7039089 1292724       20.0000    24       19.0000       20.7573 2.56e+08    9.25%
 7072254 1283943    infeasible             19.0000       20.7430 2.58e+08    9.17%
 7107348 1281824    infeasible             19.0000       20.7274 2.58e+08    9.09%
 7143017 1262286    infeasible             19.0000       20.7002 2.60e+08    8.95%
 7179356 1248370       20.0000    29       19.0000       20.6153 2.61e+08    8.50%
 7218255 1222703    infeasible             19.0000       20.2398 2.62e+08    6.53%
*7227772 1220913      integral     0       20.0000       20.1541 2.62e+08    0.77%
Found incumbent of value 19.999999 after 1932.55 sec. (890057.09 ticks)
*7247324 926397      integral     0       20.0000       20.0058 2.63e+08    0.03%
                                                     Cone: 48                  
Found incumbent of value 20.000000 after 1940.95 sec. (892272.37 ticks)
 7252826 926392    infeasible             20.0000       20.0058 2.63e+08    0.03%
Elapsed time = 1945.70 sec. (892962.40 ticks, tree = 766.57 MB, solutions = 28)
 7280114 64647    infeasible             20.0000       20.0002 2.64e+08    0.00%
 7303087 40831        cutoff             20.0000       20.0000 2.64e+08    0.00%
 7323676 15673    infeasible             20.0000       20.0000 2.66e+08    0.00%

Cover cuts applied:  1
Mixed integer rounding cuts applied:  20
Gomory fractional cuts applied:  10
Cone linearizations applied:  3328

Root node processing (before b&c):
  Real time             =    0.25 sec. (50.77 ticks)
Parallel b&c, 16 threads:
  Real time             = 1995.13 sec. (905620.31 ticks)
  Sync time (average)   =  265.77 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) = 1995.38 sec. (905671.08 ticks)
MIQCP status(101): integer optimal solution


The solution looks like:


----    156 VARIABLE assign.L  client to facility

           facility1   facility2   facility3   facility4

client1        1.000
client3                                1.000
client4                    1.000
client5                    1.000
client6                                            1.000
client8                    1.000
client9                                1.000
client10                                           1.000
client11                               1.000
client12                               1.000
client14                   1.000
client15                                           1.000
client16                   1.000
client17       1.000
client18                               1.000
client20                                           1.000
client21                               1.000
client22                               1.000
client23                               1.000
client25                                           1.000


----    156 VARIABLE facility.L  coordinates

                    x           y

facility1      17.146      84.398
facility2      23.055      66.584
facility3      29.228      22.420
facility4      66.892      43.536


----    156 VARIABLE clientUsed.L  client is assigned to facility

client1  1.000,    client3  1.000,    client4  1.000,    client5  1.000
client6  1.000,    client8  1.000,    client9  1.000,    client10 1.000
client11 1.000,    client12 1.000,    client14 1.000,    client15 1.000
client16 1.000,    client17 1.000,    client18 1.000,    client20 1.000
client21 1.000,    client22 1.000,    client23 1.000,    client25 1.000


----    156 PARAMETER DistanceReport  distances for assignments

           facility1   facility2   facility3   facility4

client1        0.084
client3                                0.029
client4                   22.499
client5                   23.346
client6                                           35.973
client8                   10.337
client9                               13.634
client10                                           0.020
client11                              14.488
client12                              17.746
client14                   0.047
client15                                          17.008
client16                  20.391
client17       3.082
client18                               6.750
client20                                           4.991
client21                              16.154
client22                              17.930
client23                               6.298
client25                                          17.089


----    156 PARAMETER CostAllocation  cost = unit cost * distance for assignment
                                      s

           facility1   facility2   facility3   facility4

client1        0.557
client3                                0.180
client4                   63.867
client5                   20.177
client6                                           36.877
client8                   56.367
client9                                4.298
client10                                           0.155
client11                              10.542
client12                              31.173
client14                   0.352
client15                                          30.295
client16                   6.962
client17      18.036
client18                              41.934
client20                                          17.903
client21                              39.260
client22                              44.184
client23                               8.219
client25                                          64.929


----    156 PARAMETER CostReport  total cost vs budget

total 496.267,    limit 500.000

Basically: we can cover 20 clients with the given budget and 4 facilities. We recalculated the cost, as the optimization model may overestimate the distance and the cost (it only worries about the cost to be below the budget). In the picture below we see the assigned clients (and their corresponding facility) as colored clusters. Remaining unassigned clients are drawn as black points.

Solution

One thing we can see is that a small distance does not guarantee to be selected: it can be too costly because of high unit cost. This is visible in the light blue cluster: a far away point is included while a closer by point is not included. We should also note that this is not a min cost solution, just a solution that obeys the budget constraint. We are not minimizing the cost.

Symmetry
In the models above we can reduce symmetry by having the locations \(\mathit{facility}_{j,k}\) ordered by the \(x\)-coordinate. We can do this by adding the (linear) constraint \[\mathit{facility}_{j,x} \le \mathit{facility}_{j+1,x}\] This makes the solution easier to read, but more importantly will remove a lot of symmetry. This can reduce solution times.

4. Approximation: restrict locations to grid points


In the previous models, we calculated two quantities:

  1. \(x,y\) location of each facility (continuous variables)
  2. Assignment of clients to a facility (binary variable) 

This turned out to be a difficult model to solve. The usual approach is: if a model is too difficult to solve, solve a different problem. If we set up a grid, and allow a facility only to be placed on the grid points, we get a much simpler model. We can pre-compute the distance and cost between each client location and each possible facility location. The end result is just a linear, but large MIP model.

The model looks like:

MIP Problem
\[\begin{align}\max & \sum_i \color{darkred}{\mathit{ClientUsed}}_i\\ & \color{darkred}{\mathit{ClientUsed}}_i = \sum_p \color{darkred}{\mathit{Assign}}_{i,p} \\ & \sum_p \color{darkred}{\mathit{PointUsed}}_p = \color{darkblue}{\mathit{NumFacilities}}\\ & \color{darkred}{\mathit{Assign}}_{i,p} \le \color{darkblue}{\mathit{NumClients}} \cdot\color{darkred}{\mathit{PointUsed}}_p\\ & \sum_{i,p} \color{darkblue}{\mathit{Cost}}_{i,p} \cdot \color{darkred}{\mathit{Assign}}_{i,p} \le \color{darkblue}{\mathit{Budget}} \\ & \color{darkred}{\mathit{Assign}}_{i,p}, \color{darkred}{\mathit{PointUsed}}_p, \color{darkred}{\mathit{ClientUsed}}_i    \in \{0,1\} \end{align}\]

Here \(p\) are the grid points. As we have box of \([0,100]\times[0,100]\) for the client locations, I used \(101 \times 101 = 10201\) grid points. This is a fine grid, so we expect similar results as for the models where facilities can be located anywhere on the \(x\)-\(y\) plane. The parameter \(\mathit{Cost}_{i,p}\) is precalculated. The resulting model is very large (considering the original data set was small):


MODEL STATISTICS

BLOCKS OF EQUATIONS           5     SINGLE EQUATIONS      255,053
BLOCKS OF VARIABLES           4     SINGLE VARIABLES      265,252
NON ZERO ELEMENTS     1,030,352     DISCRETE VARIABLES    265,251

However it solves very fast (less than a minute). The results are:

----    121 VARIABLE assign.L  client to facility

              gp1803      gp2953      gp6043      gp6811

client1        1.000
client3                    1.000
client4        1.000
client5                    1.000
client6                                            1.000
client9                    1.000
client10                                           1.000
client11                   1.000
client12                   1.000
client13                               1.000
client15                                           1.000
client16                   1.000
client17       1.000
client18                   1.000
client19                               1.000
client20                                           1.000
client21                   1.000
client22                   1.000
client23                   1.000
client25                                           1.000


----    121 VARIABLE clientUsed.L  client is assigned to gridpoint

client1  1.000,    client3  1.000,    client4  1.000,    client5  1.000,    client6  1.000,    client9  1.000
client10 1.000,    client11 1.000,    client12 1.000,    client13 1.000,    client15 1.000,    client16 1.000
client17 1.000,    client18 1.000,    client19 1.000,    client20 1.000,    client21 1.000,    client22 1.000
client23 1.000,    client25 1.000


----    121 VARIABLE numClients.L          =       20.000  objective variable

----    121 PARAMETER DistanceReport  distances for assignments

              gp1803      gp2953      gp6043      gp6811

client1        0.696
client3                    0.635
client4       17.994
client5                   35.027
client6                                           36.025
client9                   13.202
client10                                           0.546
client11                  14.002
client12                  17.751
client13                               0.126
client15                                          16.483
client16                  32.622
client17       2.452
client18                   6.111
client19                              10.735
client20                                           5.372
client21                  16.678
client22                  18.504
client23                   6.820
client25                                          16.573


----    121 PARAMETER CostAllocation  cost = unit cost * distance for assignments

              gp1803      gp2953      gp6043      gp6811

client1        4.599
client3                    3.981
client4       51.078
client5                   30.272
client6                                           36.931
client9                    4.162
client10                                           4.328
client11                  10.189
client12                  31.181
client13                               0.661
client15                                          29.360
client16                  11.137
client17      14.346
client18                  37.964
client19                              41.799
client20                                          19.269
client21                  40.534
client22                  45.598
client23                   8.901
client25                                          62.968


----    121 PARAMETER CostReport  total cost vs budget

total 489.260,    limit 500.000


We see again that we can cover 20 clients with the given budget.

Obviously the problem becomes smaller and easier when we use a coarser grid. But we give up some precision. I tried with a \(19 \times 19 = 361\) grid, with coordinates along the axes: \(5,10,\dots,90,95\). This model had only 9,411 binary variables and solved instantaneously. But that optimal solution covered only 19 clients.

Of course this we can also use this approximation model as a way to get a good initial solution. We can solve the MISOCP model using the MIPSTART option. This can shave some time from a solve from scratch.

Conclusion


A variant of the facility location problem seeks to find how many clients we can serve for a given budget. The costs are expressed as \(c_i \times \mathbf{dist}(i,j)\), where \(c_i\) is the unit cost for client \(i\) and \(\mathbf{dist}(i,j)\) is the distance between client \(i\) and facility \(j\).

We formulated the model in four different ways:

  • An MINLP model: this is the simplest model, and is a good start formalizing the problem,
  • an MIQCP model: this model needs a non-convex solver, so limited progress compared to the previous model, 
  • an MISOCP model: this model can be solved with convex solvers, but at the expense of being more complex,
  • and an approximation: facilities are only to be placed on grid points instead of on any point \((x,y)\). This yielded a large but simple MIP model.

The MIP model was by far the easiest to solve. This is not unusual: quadratic (and second order) models with discrete variables are often much more difficult to solve than linear ones.

References


  1. Uncapacitated facility location problem with at most k facilities and budget b, https://stackoverflow.com/questions/59327296/uncapacitated-facility-location-problem-with-at-most-k-facilities-and-budget-b
  2. Solving a facility location problem as an MIQCP, https://yetanothermathprogrammingconsultant.blogspot.com/2018/01/solving-facility-location-problem-as.html
  3. Jared Erickson, Robert Fourer, Detection and Transformation of Second-Order Cone Programming Problems in a General-Purpose Algebraic Modeling Language, 2019,  http://www.optimization-online.org/DB_HTML/2019/05/7194.html

Appendix: MIPstart can disappoint


The MISOCP model is difficult. The final optimal solution of 20 covered clients is found very late in the search. Below is the plot of the MIP bounds: the best found integer solution and the best possible objective:

MISOCP model: MIP bounds
 We see the integer solutions (light curve) jumping by one (the objective is integer valued). The optimal solution of 20 is not found before spending around 1800 seconds.

In the next picture, I show what happens if we find a solution with an objective of 19 very quickly using the approximation model discussed in the post, and use a MIPstart to feed this into the MISOCP solver. Finding the objective of 19 took us 1000 seconds when solving the MISOCP model from scratch. See  the picture above. So this should be a promising approach. Unfortunately, what we see is:


MISOCP model solved with MIPSTART
The total solution time has increased. We see the best integer solution of 19 is accepted. However the time to get to 20 is longer than expected. MIPstarts often help in increasing performance, but, as shown here, not always.

Monday, December 16, 2019

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