Cvxpy: cvxcore limits nnz of constraint matrix to 2^32

Created on 25 Sep 2019 · 31Comments · Source: cvxgrp/cvxpy

Describe the bug
Hi everyone. First of all, thank you all for the great tool you have built. I am facing a problem trying to run a fairly large optimization problem in my opinion. You can see the code below. The variable b size is 500 x 96. What I am trying to do is to match a sum of timeseries profiles (351236 x 15 min timesteps) with a bigger profile by minimizing their difference. With the same formulation and a much smaller problem (672 timesteps and a b variable of the size 10 x 5) the problem is solved in under 2 seconds without a problem. But when I am running it for the full scale problem I get the error you see below.

I am running this on Jupyter Lab and python 3.7.4. The python installation is done with conda.

To Reproduce

X_opt = cp.Constant(np.asarray(X.iloc[:,:500])) # the array size is (35136,500)
K_opt = cp.Constant(np.asarray(K.YearlyDemand)) # the vector size is 96
b = cp.Variable((500,96),boolean = True, value = np.zeros((500,96)))
Y_opt = cp.Constant(np.asarray(y)) # the vector size is 35136

constraints = []

constraints.append( cp.sum(b, axis = 0) == 1 ) # the sum of the elements of every column of b must be equal to 1
constraints.append( cp.sum(b, axis = 1) <= 1 ) # the sum of the elements of every row of b must be smaller or equal to 1

objective = cp.Minimize(cp.sum(cp.abs(Y_opt-cp.sum((cp.diag(K_opt)*((X_opt@b).T)).T, axis = 1))))

prob = cp.Problem(objective, constraints)

prob.solve(solver = cp.GLPK_MI, verbose = True)

Expected behavior
I would expect the problem to solve as with the much smaller problem. But when I run this one, RAM usage explodes up to 100 GB (about 99% of the available RAM on the server). After a while the RAM usage goes down and then a periodical swinging begins (RAM goes up and down from 50% to 100% every few minutes). From the error and after a lot of googling my suspicion is that the problem is too big for the memory and that at some point data is getting broken down to smaller pieces. I do not think it reaches to the point, where the solver does its work. I tried to optimize the code by vectorizing everything (current version) and trying not to have loops etc. in the formulation. But this did not change anything. Do you guys have any clue if this is a bug or a limitation? Or do you maybe have an idea on how to solve this?

Output

ValueError Traceback (most recent call last)
in

D:Anaconda3envspy37DuALlibsite-packagescvxpyproblemsproblem.py in solve(self, args, *kwargs)
287 else:
288 solve_func = Problem._solve
--> 289 return solve_func(self, args, *kwargs)
290
291 @classmethod

D:Anaconda3envspy37DuALlibsite-packagescvxpyproblemsproblem.py in _solve(self, solver, warm_start, verbose, parallel, gp, qcp, **kwargs)
567 self._construct_chains(solver=solver, gp=gp)
568 data, solving_inverse_data = self._solving_chain.apply(
--> 569 self._intermediate_problem)
570 solution = self._solving_chain.solve_via_data(
571 self, data, warm_start, verbose, kwargs)

D:Anaconda3envspy37DuALlibsite-packagescvxpyreductionschain.py in apply(self, problem)
63 inverse_data = []
64 for r in self.reductions:
---> 65 problem, inv = r.apply(problem)
66 inverse_data.append(inv)
67 return problem, inverse_data

D:Anaconda3envspy37DuALlibsite-packagescvxpyreductionsmatrix_stuffing.py in apply(self, problem)
98 # Batch expressions together, then split apart.
99 expr_list = [arg for c in cons for arg in c.args]
--> 100 Afull, bfull = extractor.affine(expr_list)
101 if 0 not in Afull.shape and 0 not in bfull.shape:
102 Afull = cvxtypes.constant()(Afull)

D:Anaconda3envspy37DuALlibsite-packagescvxpyutilitiescoeff_extractor.py in affine(self, expr)
76 size = sum([e.size for e in expr_list])
77 op_list = [e.canonical_form[0] for e in expr_list]
---> 78 V, I, J, b = canonInterface.get_problem_matrix(op_list, self.id_map)
79 A = sp.csr_matrix((V, (I, J)), shape=(size, self.N))
80 return A, b.flatten()

D:Anaconda3envspy37DuALlibsite-packagescvxpycvxcorepythoncanonInterface.py in get_problem_matrix(linOps, id_to_col, constr_offsets)
65
66 # Unpacking
---> 67 V = problemData.getV(len(problemData.V))
68 I = problemData.getI(len(problemData.I))
69 J = problemData.getJ(len(problemData.J))

D:Anaconda3envspy37DuALlibsite-packagescvxpycvxcorepythoncvxcore.py in getV(self, values)
320
321 def getV(self, values):
--> 322 return _cvxcore.ProblemData_getV(self, values)
323
324 def getI(self, values):

ValueError: negative dimensions are not allowed

Version

OS: Windows Server 2012 R2
CVXPY Version: 1.0.25

Additional context
Add any other context about the problem here.

Source

aselviar

Most helpful comment

If you are in a hurry, then your model is sufficiently simple that could do in Mosek optimizer API for Python.

erling-d-andersen on 9 Oct 2019

👍2

All 31 comments

I tried to run the program on a VM on the cloud with 964 GB of RAM. Same error occurs. May not have to do with RAM itself, but may have to do with the limits of something, which I cannot get.

aselviar on 27 Sep 2019

A negative variable dimension could easily be the result of an integer overflow somewhere. In particular, problemData.V could contain essentially all nonzero entries of the constraint matrix in the compiled version of your optimization problem. The number of such entries can be several orders of magnitude larger than the number of variables in your problem, and could overflow a 32bit int.

For the VM with 964GB RAM, was that a Windows server as well? The default compiler on Windows turns the long C datatype into 32bit ints, while the default compilers for Linux and OSX keep the long C datatypes as 64 bits. (See https://stackoverflow.com/questions/384502/what-is-the-bit-size-of-long-on-64-bit-windows for details.)

rileyjmurray on 27 Sep 2019

@rileyjmurray the vm instance was a machine on google cloud running on Ubuntu 18.04 LTS.

aselviar on 30 Sep 2019

@SteveDiamond @akshayka it seems that the C++ implementations of getV, getI, and getJ accept ints, not longs.

https://github.com/cvxgrp/cvxpy/blob/b3820156336e035b7f48be6af322f20c2c2432d3/cvxpy/cvxcore/src/ProblemData.hpp#L60

As a result, highly vectorized models (such as the one used by @aselviar) might produce runtime errors when the canonicalized problem's constraint matrix has more than 2^32 nonzero entries. Although this only happens for large, vectorized problems, such a scale is of practical importance. I suggest changing cvxcore so that these functions accept long inputs. Such a change should resolve this issue on Linux and OSX, as well as Windows machines with decent C++ compilers.

Edit: I ran the following code with a Jupyter notebook, attached to my compute server (which has 256GB RAM). It didn't explicitly print any error, but the kernel did die shortly after reaching 50% RAM utilization.

import cvxpy
import numpy as np

n = 2**3
m = int(2**32 / n) + 1

A = np.ones(shape=(m, n))
x = cvxpy.Variable(shape=(n,))
cons = [A @ x >= 0]
prob = cvxpy.Problem(cvxpy.Maximize(0), cons)

rileyjmurray on 30 Sep 2019

@aselviar If you really want to dig into this, you can modify the source code for cvxpy\cvxcore\python\canonInterface.get_problem_matrix(linOps, id_to_col, constr_offsets), and have it print out len(problemData.V) just prior to line 67. I suspect it's larger than 2^32.

That said-- even if this is the issue, and even if Steven or Akshay push a fix for this, it might not be possible to solve your problem. It seems that not only is your problem large (requiring tons of RAM), but it also has 50,000 binary variables. It is extremely difficult to provably solve general discrete problems at that scale.

rileyjmurray on 30 Sep 2019

If glpk fails then try one of the commercial options.

erling-d-andersen on 30 Sep 2019

@rileyjmurray I printed out len(problemData.V) and the number is 3'373'222'272 which is indeed > than the max int32 (2,147,483,647). If you could fix this problem and the problem is "compiled" (don't know if I use the word right), then I could also try with commercial options, as @erling-d-andersen said.

aselviar on 30 Sep 2019

👍1

Thanks for the update @aselviar . The thing is, I've never personally touched the cvxcore code before. This kind of change should probably be handled by Steven Diamond.

@SteveDiamond it would be great if a small fix for this could be pushed on the 1.0 branch. If you are pretty swamped right now, then maybe a better goal is to fix this with version 1.1.

rileyjmurray on 30 Sep 2019

👍1

Hi everyone! Is there any update on the topic?

aselviar on 8 Oct 2019

If you are in a hurry, then your model is sufficiently simple that could do in Mosek optimizer API for Python.

erling-d-andersen on 9 Oct 2019

👍2

@aselviar, would you like to send a pull request that changes the integers to longs?

akshayka on 10 Oct 2019

Careful with long integers and anaconda numpy. This might be related: https://github.com/oxfordcontrol/osqp-python/issues/14#issuecomment-494023914

bstellato on 11 Oct 2019

@rileyjmurray @bstellato this issue has a bit of a history. It's come up before. There may be a reason ints are 32 bit now. I guess we can try the change and see what happens with the CI.

SteveDiamond on 11 Oct 2019

@aselviar can you share your exact problem data with us? I'm running my cooked-up example, but the error it encountered wasn't actually the same as yours. (Mine creates a segfault somewhere within Eigen.)

rileyjmurray on 30 Mar 2020

Hello there @rileyjmurray and others. I'm encountering the same issue with int32 overflowing. I didn't think my problem was all that large: I'm trying to run a vectorized optimisation over hourly data for a year (8760 elements in each vector), and I don't have an excessive amount of constraints, or so I thought.

I thought I'd add some code that reproduces the error, since that was the last comment on here, and I'd definitely appreciate a resolution. Other than that, this package has been great, thank you!

import numpy as np
import cvxpy as cp

# Timestamps
t_start = 0
T = 4018                                          # change to T = 4019 for error
timestamp = np.arange(8760)[t_start:t_start+T]
SHAPE = timestamp.shape

# Primary input vectors
R = np.sin(timestamp)[t_start:t_start+T]
C = np.random.normal(0.5, 0.3, SHAPE)[t_start:t_start+T]

X_max = 100
X_min = 50
X_initial = X_max

# Decision variables
BG = cp.Variable(SHAPE, name="BG", nonneg=True)
GB = cp.Variable(SHAPE, name="GB", nonneg=True)
PG = cp.Variable(SHAPE, name="PG", nonneg=True)
PB = cp.Variable(SHAPE, name="PB", nonneg=True)
X = cp.Variable(SHAPE, name="X", nonneg=True)
B = cp.Variable(SHAPE, name="B")
c = cp.Variable(SHAPE, name="c", integer=True)

# Constraints
constr = []
constr.append(PG + PB >= 0)
constr.append(PG + PB <= R)
constr.append(BG <= 1 - c)
constr.append(GB <= c)
constr.append(PB <= c)
constr.append(PG <= 1 - c)
constr.append(BG + PG <= 1)
constr.append(GB <= 1)
constr.append(X >= X_max)
constr.append(X <= X_min)
constr.append(c >= 0)
constr.append(c <= 1)
constr.append(GB + PB - BG == B)
constr.extend([X[t-1] + B[t] == X[t] if t>0 else X_initial + B[t] == X[t] for t in range(X.size)])

# Solve
opt_fun = cp.Maximize(cp.sum((cp.multiply(BG, C) + PG) - (cp.multiply(GB, C) + PB)))
prob = cp.Problem(opt_fun, constr)
result = prob.solve(solver='GLPK_MI')

I've also added the following line at line 369 in cvxpy\cvxcore\python\canonInterface.get_problem_matrix(), just before the line A = scipy.sparse.csc_matrix(.... This prints out the calculation of the size of the problem, from what I can see:

print("constr_length*(var_length+1), param_size_plus_one\n",
      "= ({}*({}+1), {})\n =".format(constr_length, var_length, param_size_plus_one),
      (constr_length*(var_length+1), param_size_plus_one))

When I run the code, I get the following output and it completes without error (though the solution is infeasible):

constr_length*(var_length+1), param_size_plus_one
 = (1*(28126+1), 1)
 = (28127, 1)
constr_length*(var_length+1), param_size_plus_one
 = (76342*(28126+1), 1)
 = (2147271434, 1)

But when I change T = 4019 at the beginning, I get the following output, showing the int32 overflow.

constr_length*(var_length+1), param_size_plus_one
 = (1*(28133+1), 1)
 = (28134, 1)
constr_length*(var_length+1), param_size_plus_one
 = (76361*(28133+1), 1)
 = (-2146626922, 1)

This also then results in the error on the following line, yielding:
ValueError: 'shape' elements cannot be negative
as was also described in #792 .

I'm working on testing this in linux, but it's a bit longer to do that at my work. I am currently on Windows 10.

justnigel on 17 Jun 2020

Do you have the latest cvxpy? Try installing from source.

SteveDiamond on 18 Jun 2020

Ok it'll have to wait till next week to test out. I don't typically install from source. And it seems like it worked on linux.
Thanks for your reply.

justnigel on 22 Jun 2020

I have also met this issue and already solved it. The problem is that when calculating the shape for scipy.sparse.csc_matrix, it always gets a negative value, since the shape is numpy.int32 and if the shape is too large, we get overflow problem. The idea to solve it is very simple, just convert it into numpy.int64.

So do the following steps:
1,find canonInterface.py in cvxpycvxcorepython and find get_problem_matrix function. you will see around 368line, use this instead A = scipy.sparse.csc_matrix((V, (I, J)), shape=(np.int64(constr_length)*np.int64(var_length+1), param_size_plus_one))
2, find conic_solver.py in cvxpyreductionssolversconic_solvers and find format_constraints function. you will see around 225 line and use this instead restructured_A = restructured_A.reshape(np.int64(restruct_mat.shape[0]) * np.int64(problem.x.size + 1),problem.A.shape[1], order='F')

YongchangHu on 24 Jun 2020

Already solved it.
See https://github.com/cvxgrp/cvxpy/issues/826#issuecomment-648618636

YongchangHu on 24 Jun 2020

@rileyjmurray unfortunately I cannot upload the data and the problem formulation, since they are both confidential i.a. I am not the data owner

aselviar on 24 Jun 2020

I have also met this issue and already solved it. The problem is that when calculating the shape for scipy.sparse.csc_matrix, it always gets a negative value, since the shape is numpy.int32 and if the shape is too large, we get overflow problem. The idea to solve it is very simple, just convert it into numpy.int64.

So do the following steps:
1,find canonInterface.py in cvxpycvxcorepython and find get_problem_matrix function. you will see around 368line, use this instead A = scipy.sparse.csc_matrix((V, (I, J)), shape=(np.int64(constr_length)*np.int64(var_length+1), param_size_plus_one))
2, find conic_solver.py in cvxpyreductionssolversconic_solvers and find format_constraints function. you will see around 225 line and use this instead restructured_A = restructured_A.reshape(np.int64(restruct_mat.shape[0]) * np.int64(problem.x.size + 1),problem.A.shape[1], order='F')

I have been struggling with this error for some time. First, creating a 3.6 conda environment solved the problem. Second, YongchangHu's advice solved my issue with the 3.7 environment. If it hasn't already been done, it would be great to include the changes mentioned above in the next revision.

teyber on 21 Jul 2020

@teyber the solution mentioned by @YongchangHu has been implemented on CVXPY's master branch. The next release will contain this fix. I haven't closed this issue, because I'm not entirely convinced that the issue is resolved in CVXPY's C++ backend.

rileyjmurray on 21 Jul 2020

Similar error for quadratic programming

Hello, I am trying to solve a quadratic program which is very similar to the portfolio optimization you have here:
https://www.cvxpy.org/examples/basic/quadratic_program.html
the only difference is that my matrix Sigma is 64000x64000 sparse matrix of type '<class 'numpy.float64'>' with 1533280 stored elements in Compressed Sparse Column format.

The program works fine when the matrix has dimension 27000x27000

It seems like the error comes from the qp_matrix_stuffing.py and coeff_extractor.py

I was wondering if there is a fix similar to the one proposed by @YongchangHu for quadratic programming. it seems that the fix you have is only for conic_solvers

Below are the complete error messages:

First one is a Warning
C:\Users\Eli\Anaconda3\lib\site-packages\cvxpy\utilities\coeff_extractor.py:266: RuntimeWarning: overflow encountered in long_scalars acc_height += P_height * P_shape[1]

The second one is the one that is similar to the one encountered in this post:
`---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
in
1 t = time.time()
----> 2 prob.solve()
3 elapsed = time.time() - t
4
5 #print output

~Anaconda3libsite-packagescvxpyproblemsproblem.py in solve(self, args, *kwargs)
394 else:
395 solve_func = Problem._solve
--> 396 return solve_func(self, args, *kwargs)
397
398 @classmethod

~Anaconda3libsite-packagescvxpyproblemsproblem.py in _solve(self, solver, warm_start, verbose, gp, qcp, requires_grad, enforce_dpp, **kwargs)
743
744 data, solving_chain, inverse_data = self.get_problem_data(
--> 745 solver, gp, enforce_dpp)
746 solution = solving_chain.solve_via_data(
747 self, data, warm_start, verbose, kwargs)

~Anaconda3libsite-packagescvxpyproblemsproblem.py in get_problem_data(self, solver, gp, enforce_dpp)
523 inverse_data = self._cache.inverse_data + [solver_inverse_data]
524 else:
--> 525 data, inverse_data = solving_chain.apply(self)
526 safe_to_cache = (
527 isinstance(data, dict)

~Anaconda3libsite-packagescvxpyreductionschain.py in apply(self, problem)
69 inverse_data = []
70 for r in self.reductions:
---> 71 problem, inv = r.apply(problem)
72 inverse_data.append(inv)
73 return problem, inverse_data

~Anaconda3libsite-packagescvxpyreductionsqp2quad_formqp_matrix_stuffing.py in apply(self, problem)
261 extractor = CoeffExtractor(inverse_data)
262 params_to_P, params_to_q, flattened_variable = self.stuffed_objective(
--> 263 problem, extractor)
264 # Lower equality and inequality to Zero and NonPos.
265 cons = []

~Anaconda3libsite-packagescvxpyreductionsqp2quad_formqp_matrix_stuffing.py in stuffed_objective(self, problem, extractor)
245 # extract to 0.5 * x.T * P * x + q.T * x + r
246 expr = problem.objective.expr.copy()
--> 247 params_to_P, params_to_q = extractor.quad_form(expr)
248 # Handle 0.5 factor.
249 params_to_P = 2*params_to_P

~Anaconda3libsite-packagescvxpyutilitiescoeff_extractor.py in quad_form(self, expr)
267
268 # Stitch together Ps and qs and constant.
--> 269 P = sp.coo_matrix((vals, (rows, cols)), shape=(acc_height, num_params))
270 # Stack q with constant offset as last row.
271 q = np.vstack(q_list)

~Anaconda3libsite-packagesscipysparsecoo.py in __init__(self, arg1, shape, dtype, copy)
152 # Use 2 steps to ensure shape has length 2.
153 M, N = shape
--> 154 self._shape = check_shape((M, N))
155
156 idx_dtype = get_index_dtype(maxval=max(self.shape))

~Anaconda3libsite-packagesscipysparsesputils.py in check_shape(args, current_shape)
283 raise ValueError('shape must be a 2-tuple of positive integers')
284 elif new_shape[0] < 0 or new_shape[1] < 0:
--> 285 raise ValueError("'shape' elements cannot be negative")
286
287 else:

ValueError: 'shape' elements cannot be negative`

Any help is appreciated! Thank you

enegrini on 28 Sep 2020

Are you able to try the fix yourself? Just write acc_height += np.int64(P_height) * np.int64(P_shape[1]).

SteveDiamond on 28 Sep 2020

Are you able to try the fix yourself? Just write acc_height += np.int64(P_height) * np.int64(P_shape[1]).

Yup, I just changed it, I'll try to run it again!
Thank you!

enegrini on 28 Sep 2020

Are you able to try the fix yourself? Just write acc_height += np.int64(P_height) * np.int64(P_shape[1]).

Yup, I just changed it, I'll try to run it again!
Thank you!

So I ran it again and got the following:
`---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
in
1 t = time.time()
----> 2 prob.solve()
3 elapsed = time.time() - t
4
5 #print output

~Anaconda3libsite-packagescvxpyproblemsproblem.py in solve(self, args, *kwargs)
394 else:
395 solve_func = Problem._solve
--> 396 return solve_func(self, args, *kwargs)
397
398 @classmethod

~Anaconda3libsite-packagesscipysparsecoo.py in __init__(self, arg1, shape, dtype, copy)
190 self.data = self.data.astype(dtype, copy=False)
191
--> 192 self._check()
193
194 def reshape(self, args, *kwargs):

~Anaconda3libsite-packagesscipysparsecoo.py in _check(self)
281 raise ValueError('column index exceeds matrix dimensions')
282 if self.row.min() < 0:
--> 283 raise ValueError('negative row index found')
284 if self.col.min() < 0:
285 raise ValueError('negative column index found')

ValueError: negative row index found`

Thank you for helping out!

enegrini on 28 Sep 2020

It would be most helpful if you can post a code example that produces the bug. Otherwise you can just keep checking the different values that go into the construction of P, and see if they're 32 bit ints.

SteveDiamond on 28 Sep 2020

Here is the code, I will also look what goes into the construction of P.
This was a Jupyter notebook, I tried to upload the file but it was not allowed, so sorry for the terrible impagination.

!/usr/bin/env python

coding: utf-8

In[1]:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import copy
from scipy import sparse
from scipy.linalg import toeplitz
from scipy.sparse import lil_matrix, csr_matrix, identity, coo_matrix
from scipy.sparse.linalg import LinearOperator
from scipy.sparse.linalg import eigs
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.animation as animation
import cvxpy as cp
import time
from numpy import linalg as LA

Define matrix A to minimize $frac{1}{2} x^{T}A^{T}Ax$

I use the code from the file FPE.m until line 60

In[2]:

Define prameters

Nx = 40 #original code had Nx=50, but it takes more to run
Lx = 14.
Ly = 20.
Lz = 20.

dx = Lx/(Nx-1)
dy = Ly/(Nx-1)
dz = Lz/(Nx-1)

vx = np.linspace(-Lx/2.,Lx/2.,Nx)
vy = np.linspace(-Ly/2.,Ly/2.,Nx)
vz = np.linspace(25.-Lz/2.,25.+Lz/2.,Nx)

sigma = 3.
rho = 26.5
beta = 0.16
gamma = 0.02

In[3]:

iden = identity(Nx, dtype='float', format='csr')

N = np.linspace(0,Nx-1,Nx) #start from 0 instead of 1 because of different counting in MATLAB and Py

In[4]:

xcoord = coo_matrix((vx, (N, N)), shape=(Nx, Nx))
ycoord = coo_matrix((vy, (N, N)), shape=(Nx, Nx))
zcoord = coo_matrix((vz, (N, N)), shape=(Nx, Nx))

In[5]:

A = np.zeros((1,Nx))
A[0,1] = 1
D1 = toeplitz(-A,A)
D1 = coo_matrix(D1)

B = np.zeros((1,Nx))
B[0,0] = -2
B[0,1] = 1
D2 = toeplitz(B,B)
D2 = coo_matrix(D2)

In[6]:

Kronecker products

use coordinate format like MATLAB code

identity = sparse.kron(sparse.kron(iden, iden, format = 'coo'), iden, format = 'coo')

x = sparse.kron(sparse.kron(xcoord, iden, format = 'coo'), iden, format = 'coo')
y = sparse.kron(sparse.kron(iden, ycoord, format = 'coo'), iden, format = 'coo')
z = sparse.kron(sparse.kron(iden, iden, format = 'coo'), zcoord, format = 'coo')

Dx = sparse.kron(sparse.kron(D1, iden, format = 'coo'), iden, format = 'coo')/(2.dx)
Dy = sparse.kron(sparse.kron(iden, D1, format = 'coo'), iden, format = 'coo')/(2.dy)
Dz = sparse.kron(sparse.kron(iden, iden, format = 'coo'), D1, format = 'coo')/(2.*dz)

D2x = sparse.kron(sparse.kron(D2, iden, format = 'coo'), iden, format = 'coo')/dx2
D2y = sparse.kron(sparse.kron(iden, D2, format = 'coo'), iden, format = 'coo')/dy2
D2z = sparse.kron(sparse.kron(iden, iden, format = 'coo'), D2, format = 'coo')/dz**2

In[7]:

laplacian = D2x + D2y + D2z

Define FPE operator for Lorenz attractor, this is our matrix A

In[8]:

H = Dx(sigma(x-y)) + Dy(y-rhox+xz) + Dz(betaz-xy) + gamma*laplacian
H

Use CVXPY to solve the problem:

Minimize $frac{1}{2}x^{T}Sigma x - r^{T}x $

Subject to $x geq 0$ and $textbf{1}^{T}x = 1$

In our case we have $ Sigma = A^{T}A$ and $r = textbf{0}$

see example here https://www.cvxpy.org/examples/basic/quadratic_program.html

In[9]:

Sigma = H.transpose()*H
Sigma.tocoo()

In[10]:

n = Sigma.shape[0]
x = cp.Variable(n)
prob = cp.Problem(cp.Minimize((1/2)*cp.quad_form(x, Sigma)),
[ x >= 0, #x should have non negative entries
cp.sum(x) == 1]) #x should sum to one

In[11]:

import scipy.sparse.sputils
import numpy as np

scipy.sparse.sputils.get_index_dtype((Sigma,))

In[12]:

t = time.time()
prob.solve()
elapsed = time.time() - t

print output

print('Time needed to solve problem:', elapsed)
print("nThe optimal value is", prob.value)
print("A solution x is")
print(x.value)
print("A dual solution corresponding to the inequality constraints is")
print(prob.constraints[0].dual_value)

enegrini on 29 Sep 2020

Really sorry @enegrini I'm not sure why I asked for a code sample. Your problem is on windows, so I can't actually replicate it :(
I'm afraid you'll have to try casting things in coeff_extractor to int64 and see when it goes away.

SteveDiamond on 29 Sep 2020

No worries! I'll try that! thank you!

enegrini on 30 Sep 2020

You might want to use the latest master. I made some changes to address #1135

SteveDiamond on 30 Sep 2020

Was this page helpful?

0 / 5 - 0 ratings