problem with SDP

Issue summary

The code below solves the SDP correctly in some cases. On some problems we are seeing large discrepancies between CVXPY and CVX (called from matlab). We have good reasons to believe the matlab code. Do you have any suggestions?

Reproduce the problem

First we import and load some test cases

import picos as pic
import cvxopt

# Download pickle here: 
# https://www.dropbox.com/sh/j7045mtxk0evl5n/AACx6B26hCOstF5QHHsnBoOWa?dl=0

import pickle

with open('sdp_data.pickle', 'rb') as f:
    data = pickle.load(f)

# load problems and matlab solutions as a list of np.arrays
input_matrix = [data.get('C1'), data.get('C2'), data.get('C3')]
matlab_cvx_answer = [data.get('X1'), data.get('X2'), data.get('X3')]
p_dim = data.get('pdim')
c_dim = data.get('cdim')

Now use a function to define the SDP

def picos_sdp(Cin: Union[np.ndarray, sps.csr_matrix], phdim: int, codim:int):
    # params
    C = cvxopt.matrix(Cin)
    Id = cvxopt.matrix(np.eye(phdim))
    Id = pic.new_param('Id', Id)
    C = pic.new_param('C', C)
    
    #create the problem in picos
    F = pic.Problem()
    X = F.add_variable('X',(Cin.shape[0], Cin.shape[0]),'hermitian')
    
    # Form objective
    F.set_objective('max', (pic.trace(X*C)).real)    
    
    F.add_constraint( X >> 0 )
    F.add_constraint(pic.partial_trace(X, k=1, dim=[phdim, codim])==Id)
    F.solve(verbose = 0, maxit=200000)
    print(F)
    return X, F

Loop over inputs and compare solution to matlab's solution

# Loop over three problems and compare picos to cvx

for C, Xmatlab, pdim, cdim in zip(input_matrix, matlab_cvx_answer, p_dim, c_dim):
    # solve the SDP
    X, prob = picos_sdp(C.toarray(), pdim, cdim)
    # see how close
    print('allclose', np.allclose(X.value, Xmatlab, atol=1e-07))
    print('max(abs(R-Rans))',np.max(np.abs(X.value - Xmatlab)))
    print('')

Output

---------------------
optimization problem (Complex SDP):
300 variables, 4 affine constraints, 300 vars in 1 SD cones

X 	: (24, 24), hermitian

	maximize 0.5·(trace(X·C) + trace(X·C).conj)
such that
  X ≽ 0
  trace_1(X) = Id
---------------------
allclose False
max(abs(R-Rans)) 0.5000000000296454

---------------------
optimization problem (Complex SDP):
300 variables, 4 affine constraints, 300 vars in 1 SD cones

X 	: (24, 24), hermitian

	maximize 0.5·(trace(X·C) + trace(X·C).conj)
such that
  X ≽ 0
  trace_1(X) = Id
---------------------
allclose False
max(abs(R-Rans)) 0.9988953132710172

---------------------
optimization problem (Complex SDP):
1081 variables, 4 affine constraints, 1081 vars in 1 SD cones

X 	: (46, 46), hermitian

	maximize 0.5·(trace(X·C) + trace(X·C).conj)
such that
  X ≽ 0
  trace_1(X) = Id
---------------------
allclose False
max(abs(R-Rans)) 0.9814188021784475

Version information

from platform import python_version

print('Python:',python_version())
print('Numpy:', np.version.version)
print(cvxopt.sys.version_info)
print('version: ', cvxopt.sys.version)

Python: 3.7.3
Numpy: 1.16.5
sys.version_info(major=3, minor=7, micro=3, releaselevel='final', serial=0)
version:  3.7.3 (default, Mar 27 2019, 16:54:48) 
[Clang 4.0.1 (tags/RELEASE_401/final)]