2
$\begingroup$

I am familiarizing myself with CVXPY, and encountered a strange problem. I have the following simple toy optimization problem:

import numpy as np
import cvxpy as cp

A=np.array([[1,0,0],[0,1,0], [0,0,1]])
y=np.array([1,1,1])

# Upper bound for the constraint term
upper=1
# Solve the optimization problem using CVXPY
x = cp.Variable(3)
objective = cp.Minimize(cp.sum_squares(x))
constraint = [cp.sum_squares(A*x - y) <= upper]
prob = cp.Problem(objective, constraint)
prob.solve()
optimal_x = x.value

print('Value of constraint at optimal x:' + str(np.linalg.norm(A*optimal_x - y)**2))

Now, I expect my output number to be samller than upper=1, but what I get is the following:

Value of constraint at optimal x:3.0000000068183947

I am very confused about how this could be true. Am I using the function cp.sum_squares incorrectly? Am I just setting up the optimization in a wrong way? Any help is appreciated!!

$\endgroup$
  • $\begingroup$ I suspect that this will be flagged as off-topic, see scicomp.stackexchange.com/help/on-topic, since it is "essentially about debugging a code sample, or about low-level language syntax". But anyway, one thing to check is whether A*x in cvxpy and A*optimal_x in NumPy mean the same thing. $\endgroup$ – LonelyProf Dec 26 '18 at 9:32
2
$\begingroup$

You are doing a broadcast (A*x), rather than matrix multiplication (A@x), so your code should look like this:

#!/usr/bin/env python3

import numpy as np
import cvxpy as cp

A=np.array([[1,0,0],[0,1,0], [0,0,1]])
y=np.array([1,1,1])

# Upper bound for the constraint term
upper=1
# Solve the optimization problem using CVXPY
x          = cp.Variable(3)
objective  = cp.Minimize(cp.sum_squares(x))
constraint = [cp.sum_squares(A@x - y) <= upper]
prob       = cp.Problem(objective, constraint)
prob.solve()
optimal_x = x.value

print('Value of constraint at optimal x: {0}'.format(np.linalg.norm(A@optimal_x - y)**2))

Which gives

Value of constraint at optimal x: 1.000000002699704
$\endgroup$

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.