Cvxpy: Trace function of certain 1x1 matrix expressions can cause a segfault at solve time

[GrayThomas]

My favorite library maintainers,
This following code snippet quietly segfaults in cvxpy version 1.0.1 (exit code -11):

import cvxpy as cvx
X_C = cvx.Variable((1, 1), PSD=True)*1 # segfaults
# X_C = cvx.Variable((1, 1), PSD=True)*cvx.Constant(1) # also segfaults
# X_C = cvx.Variable((1, 1), PSD=True) # fine
cons = [cvx.trace(X_C)==1]
cvx.Problem(cvx.Maximize(X_C[0,0]), cons).solve(solver="SCS", verbose=True)

All 1817 nose tests are passing. I'm on Ubuntu 16.04 (fresh install) with the anaconda distribution:

    cvxcanon-0.1.1 (cvxgrp)
    cvxpy-1.0.1 (cvxgrp)
    dill-0.2.7.1
    ecos-2.0.5 (cvxgrp)
    libgcc-7.2.0
    multiprocess-0.70.4 (cvxgrp)
    scs-1.2.6 (cvxgrp)

Valid DCP should not segfault, so I thought you'd like to know. Also, in previous versions this trick of multiplying semidefinite variable definitions by hard-coded scalings was able to effectively alter the initial step size, and was the difference between solving and not solving several problems with CVXOPT's interior point solver. (I'd try it with CVXOPT, but I'm still not passing nose tests with CVXOPT)

[SteveDiamond]

Uh oh! We will look into this.

[GrayThomas]

I neglected to check non semi-definite variables. They also segfault:

# X_C = cvx.Variable((1, 1)) # fine
# X_C = cvx.Variable((1, 1))*1 # also segfaults
# X_C = cvx.Variable((1, 1)) * cvx.Constant(1) # also segfaults

[rileyjmurray]

This issue is resolved, at least for Python 2.7 and cvxpy 1.0.5.

Example code below (no segfaults).

(cvxpy1dev) Rileys-MacBook-Pro:research RJMurray$ python
Python 2.7.14 |Anaconda, Inc.| (default, Dec  7 2017, 11:07:58) 
[GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import cvxpy
>>> x = cvxpy.Variable((1,1))
>>> y = 1*x
>>> y = x*1
>>> print(cvxpy.__version__)
1.0.5

@GrayThomas please try again on your machine (after updating to 1.0.5), and let us know if it's fixed for you as well.

[GrayThomas]

After updating to 1.0.5, that works for me (still using Python 3.6.5). But, to my surprise, my original test case still fails. Digging a little deeper, it seems to be related to the trace function---and also seems to be specific to the case where the argument of the trace function is a 1 by 1 matrix, but the error occurs when the problem is solved (from your code snippet above, it looks like I miscommunicated the location of the segfault). Having the trace of a 1x1 matrix that is the result of multiplication within either the objective or an equality constraint (inequality constraints untested) is sufficient to reproduce the segfault (of prob.solve()) behavior. I'll update the name of the issue to reflect the greater clarity I have now. I'm still a little worried that the conditions that can cause a segfault are more general than the example below (for instance, I only ever tried expressions with trace), but I only have a limited amount of time available to diagnose.

import cvxpy
import numpy as np
x = cvxpy.Variable((1,1)) # this problem does not occur with 2x2 matrices
print(cvxpy.__version__)
dangerous_expression = cvxpy.trace(x*1) # segfault on prob.solve()
# dangerous_expression = x*1 # fine
# dangerous_expression = cvxpy.trace(x) # fine
# dangerous_expression = cvxpy.trace(x*1+0) # Also fine!?
# dangerous_expression = cvxpy.trace(x*1+cvxpy.Constant(np.array([[0,],]))) # fine
# dangerous_expression = cvxpy.trace(x*1+x*1) # fine
# dangerous_expression = cvxpy.trace((x*1+x*1)*1) # segfault on prob.solve()
# dangerous_expression = cvxpy.trace(1*(x*1+x*1)) # segfault on prob.solve()
def trace(X): # quick workaround
    if X.shape==(1,1): # trivial 1x1 matrix case
        return X
    return cvxpy.trace(X)
# dangerous_expression = trace(x*1) # fine
# dangerous_expression = trace((x*1+x*1)*1) # fine
safe_expression = x
prob = cvxpy.Problem(cvxpy.Minimize(safe_expression),[dangerous_expression==1]) 
# prob = cvxpy.Problem(cvxpy.Minimize(dangerous_expression),[safe_expression==1]) # segfaults under the same conditions
prob.solve() # the segfault happens on this line, if it happens.

[rileyjmurray]

@GrayThomas thank you for providing more detail! I'm working on this now.

The error occurs on line 75 of canonInterface.py, once build_matrix(lin_vec, id_to_col_C) is called within get_problem_matrix.
Here's the full call stack right before the error.

[rileyjmurray]

You can generate the same segfault on unconstrained problems, without even using the trace function. You also don't need the variable to be PSD; any 1-by-1 variable will suffice.

The example(s) below show that indexing certain 1-by-1 expressions results in a segfault.

import cvxpy as cvx
x = cvx.Variable((1, 1))
expr0 = x * 1  # segfaults iff indexed
expr1 = x * cvx.Constant(1) # segfaults iff indexed
expr2 = x * 1 + 0 # works regardless of indexing

expr = expr0  # adjust this value during testing (try expr0, expr1, expr2, or more if you'd like)

prob = cvx.Problem(cvx.Minimize(expr), [])
prob.solve()
print('first is ok')
other_expr = expr[0,0]
prob = cvx.Problem(cvx.Minimize(other_expr), [])
prob.solve()
print('second is ok')

In order to get a better understanding of these issues, I think it wise to look at expr.canonical_form and other_expr.canonical_form.

@SteveDiamond @akshayka can you guys dig into this for a bit?

[akshayka]

Yikes, i never saw this issue. I dug into the C++ stack trace for Riley's code example:

    frame #0: 0x00000001194a5058 _cvxcore.cpython-37m-darwin.so`add_triplets(tripletList=size=0, slices=size=2, dims=size=0, axis=1, col_offset=0, row_offset=0) at LinOpOperations.cpp:571 [opt]
    frame #1: 0x00000001194a17e6 _cvxcore.cpython-37m-darwin.so`get_index_mat(lin=0x00007ffeefbfd600) at LinOpOperations.cpp:610 [opt]
    frame #2: 0x000000011949e6b8 _cvxcore.cpython-37m-darwin.so`get_func_coeffs(lin=<unavailable>) at LinOpOperations.cpp:83 [opt]
    frame #3: 0x00000001194932cc _cvxcore.cpython-37m-darwin.so`get_coefficient(lin=0x00007ffeefbfd600) at cvxcore.cpp:80 [opt]
    frame #4: 0x0000000119493d7a _cvxcore.cpython-37m-darwin.so`process_constraint(lin=0x00007ffeefbfd600, V=size=0, I=size=0, J=size=0, constant_vec=size=1, vert_offset=0x00007ffeefbfd6fc, id_to_col=size=1, horiz_offset=0x00007ffeefbfd6dc) at cvxcore.cpp:145 [opt]
    frame #5: 0x000000011949497b _cvxcore.cpython-37m-darwin.so`build_matrix(constraints=size=1, id_to_col=size=4294967295) at cvxcore.cpp:224 [opt]

In line 571 of LinOpOperations.cpp, dims is a 0-length vector and axis is 1, hence the segfault.

I'm not familiar with this code. @SteveDiamond , can you take a look when you get a chance? The issue has to do with indexing into a 1 by 1 expression.

[akshayka]

For reference, on a mac, you can inspect the C++ stack trace by dropping into lldb. After copying @rileyjmurray 's script into a file named segf.py, enter the following at the command line:

lldb python segf.py

Then, at the lldb console, type run and push enter:

(lldb) run

to obtain the following output:

Process 73810 launched: '/Users/Akshay/virtualenvs/cvxpy-dev/bin/python' (x86_64)
first is ok
_cvxcore.cpython-37m-darwin.so was compiled with optimization - stepping may behave oddly; variables may not be available.
Process 73810 stopped
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0x4)
    frame #0: 0x0000000119eea058 _cvxcore.cpython-37m-darwin.so`add_triplets(tripletList=size=0, slices=size=2, dims=size=0, axis=1, col_offset=0, row_offset=0) at LinOpOperations.cpp:571 [opt]
   568    int step = slices[axis][2];
   569    int pointer = start;
   570    while (true) {
-> 571      if (pointer < 0 || pointer >= dims[axis]) {
   572        break;
   573      }
   574      int new_offset = col_offset + pointer*vecprod_before(dims, axis);
Target 0: (python) stopped.
(lldb)

Typing bt will print the stack trace. You can print variables using the print command, and travel up and down the stack using up and down.

[SteveDiamond]

Sure, I'll look into this.

[SteveDiamond]

Ok this is fixed now. I added a special case for 1x1.