Cvxpy: Large memory usage in very sparse programs

[angeris]

Describe the bug
When generating large sparse programs, CVXPY spends a huge amount of memory canonicalizing the expression and passing it to a solver. In particular, a quadratic program with at most ~60k nonzeros in its equivalent, handcrafted affine conic map takes up < 10Mb, while CVXPY's memory usage shoots up to around 50Gb.

To Reproduce
Clone this repo and run test.py after replacing n = 1000 with n = 10000.

Version

• OS: Mac OS Catalina
• CVXPY Version: 1.0.25

Initially reported by @ludi1001

[akshayka]

If anyone wants to dig into this, here's a Python profile, visualized with Tuna, at master, n=1000. (The performance at 1.0.25 is similar.) This is a profile of execution time, but it might hint at what's using lots of memory.

And here's a cProfile dump, sorted by cumulative time:

         16153245 function calls (14648932 primitive calls) in 17.424 seconds

   Ordered by: cumulative time
   List reduced from 605 to 10 due to restriction <10>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000   17.424   17.424 problem.py:258(solve)
        1    0.000    0.000   17.424   17.424 problem.py:542(_solve)
        1    0.000    0.000   16.929   16.929 problem.py:325(get_problem_data)
        1    0.001    0.001   15.325   15.325 chain.py:47(apply)
        1    0.007    0.007    6.338    6.338 ecos_conif.py:84(apply)
        1    0.046    0.046    6.250    6.250 conic_solver.py:182(format_constraints)
        1    0.005    0.005    5.591    5.591 cone_matrix_stuffing.py:196(apply)
        2    0.026    0.013    5.123    2.561 coeff_extractor.py:60(affine)
32020/22016    0.202    0.000    3.602    0.000 compressed.py:30(__init__)
628155/74037    0.876    0.000    3.591    0.000 performance_utils.py:43(_compute_once)

& sorted by total time

         16153245 function calls (14648932 primitive calls) in 17.424 seconds

   Ordered by: internal time
   List reduced from 605 to 50 due to restriction <50>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        2    1.182    0.591    1.182    0.591 {built-in method _cvxcore.build_matrix}
  1146279    0.984    0.000    0.989    0.000 {built-in method builtins.getattr}
   128066    0.966    0.000    0.966    0.000 {method 'reduce' of 'numpy.ufunc' objects}
628155/74037    0.876    0.000    3.591    0.000 performance_utils.py:43(_compute_once)
    32020    0.589    0.000    1.783    0.000 compressed.py:138(check_format)
    90070    0.528    0.000    1.533    0.000 sputils.py:120(get_index_dtype)
   180140    0.423    0.000    0.423    0.000 getlimits.py:497(__init__)
    47009    0.419    0.000    1.169    0.000 canonInterface.py:345(make_linC_from_linPy)
   435293    0.389    0.000    0.393    0.000 {built-in method numpy.array}
    19001    0.378    0.000    0.476    0.000 leaf.py:85(__init__)
286193/286187    0.354    0.000    1.569    0.000 {built-in method numpy.core._multiarray_umath.implement_array_function}
    12011    0.307    0.000    0.878    0.000 coo.py:267(_check)
    52000    0.297    0.000    0.297    0.000 {built-in method builtins.pow}
   838282    0.255    0.000    0.421    0.000 {built-in method builtins.isinstance}
    74022    0.216    0.000    0.999    0.000 fromnumeric.py:73(_wrapreduction)
   342078    0.206    0.000    0.206    0.000 {built-in method builtins.setattr}

[SteveDiamond]

@angeris did you vectorize the code?

[rileyjmurray]

[Edited because I misunderstood @angeris's test script.]

@SteveDiamond their model is vectorized in the sense that they have a single vector-valued variable. However the number of atoms in the model is very large (about 3 * n atoms with scalar arguments, where the recommended test is with n=10000).

I don't think this model should be used as a benchmark for compile time. However, I do think this model raises an interesting question of memory pressure. Why is it that so much memory is allocated at once? What's going on with garbage collection, either at the C++ or Python level?

One possibility is that each atom ends up duplicating the stored constants, even when the same constant is used across many atoms. In @angeris's test script each atom ends up involving the same sparse matrix, with about 3*n nonzero entries. If you duplicate that matrix 3*n times with n=10000, you've got (30000)^2 * 4 bytes, or equivalently about 4GB just to store the nonzero entries in these matrices involved. Add in the overhead of sparse matrix storage, and whatever overhead is associated with cvxpy's sparse matrix conversions, and the order-of-magnitude 50GB RAM seems more reasonable. @SteveDiamond @akshayka can you comment on if duplication of the same sparse matrix across multiple atoms might be an issue?

[angeris]

This is not vectorized, but either way it's still a relatively small number of constraints (as discussed with @akshayka )—the main reason for no vectorization is that the code was supposed to be clear and follow the exact derivation in the paper.

I can add a simple vectorized version for testing purposes as well, but it would be surprising if looping over maybe 40k-80k constraints (even in Python) would yield such a large slowdown?/large memory consumption?

I haven't benchmarked anything yet since I haven't gotten the chance to dig a little deeper into the code.

Edits based on updated comments: I think that @rileyjmurray has nailed it on the head. I think things are being copied either into cvxcore or copied across atoms, instead of being passed as references and then copied only when a change needs to be performed to the matrix (although even this I'm not convinced is the perfect solution, since you may perform the same change over and over again to the same reference, e.g.).

[akshayka]

@rileyjmurray Duplication of the same sparse matrix across multiple atoms is very likely an issue. As @angeris said, I think our Python-to-C++ interface probably copies data.

In fact you can see the copy in cvxcore/src/LinOp.hpp:set_sparse_data:127

    for (int idx = 0; idx < data_len; idx++) {
      tripletList.push_back(
          Triplet(int(row_idxs[idx]), int(col_idxs[idx]), data[idx]));
    }

The binding from Python to C++ is not pretty. It was written many years ago, and the authors are no longer associated with cvxgroup.

At some point, we should probably drop the swig wrapper and use PyBind11, which makes it easy to convert scipy matrices / numpy arrays to Eigen matrices without copies (https://pybind11.readthedocs.io/en/stable/advanced/cast/eigen.html).

[akshayka]

We likely won't switch to PyBind11 anytime soon.

In the meantime, we could rewrite cvxcore/src/LinOp.hpp:set_sparse_data to map the data directly into a SparseMatrix, in order to avoid a copy (this is what set_dense_data does). See https://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html, "Mapping external buffers".

In principle, this should be easy: all you'd need to do is pass a buffer holding the data, and buffers for what SciPy calls indices and indptr (Eigen calls these InnerIndices and outerIndexPtr, respectively). We'd also need to modify the signature of set_sparse_data and pass in the indices/indptr, & modify the swig binding (cvxcore/python/cvxcore.i).

I wouldn't be surprised if there were more wasteful copies hiding in the code. But fixing this function is an important first step. Contributions are welcome!

[akshayka]

I did some of the initial work in this branch: https://github.com/cvxgrp/cvxpy/tree/fewer_copies_in_cvxcore

If anyone does want to take a go at reducing memory consumption (@angeris ?), you could work off that branch as a starting point.

(The change looks bigger than it is, because it updates Eigen.)

The relevant changes are in LinOp.hpp and cvxcore.i. The change maps sparse data held by atoms directly into cvxcore, reusing the SciPy sparse matrix buffer. There are many more copies in LinOpOperations.cpp, so this change alone might not get you enough memory savings (the functions get_constant_data, sparse_reshape_to_vec, and get_mul_mat, for example, do unnecessary copies). I wouldn't be surprised if there were copies earlier on in Python as well.

It would be a good idea to make similar changes for dense data, at some point.

[akshayka]

@angeris , I looked into this a bit more. I ran the second lower bound problem in your test script at master, with n=10000 (the one that runs after the call to complex_to_real). A few observations.

1. Canonicalization was very slow (around 15 minutes).
2. Memory consumption didn't spike to 50 GB (I only have 16 GB of RAM -- how much does your Mac have?). Python seemed to use about 3 GB of RAM at most, but I wasn't looking at top for the whole 16 minutes / I didn't profile memory consumption.
3. The objective is very large. The expression dual_tree has 630,006 nodes in it. I counted nodes via

    def count_nodes(root):
        if not root.args:
            return 1
        return 1 + sum(count_nodes(arg) for arg in root.args)
    num_nodes = count_nodes(dual_obj)
    print('num nodes: ', num_nodes)

After the dcp2cone canonicalization runs, your problem is going to be very very large (as measured by the number of expressions). Meaning the description of your problem is very large, even though the size of the numerical data is not that large.

You're probably going to need to vectorize the construction of your problem so that its description isn't so large.