-
Notifications
You must be signed in to change notification settings - Fork 44
Open
Description
Hello all,
I was following this discussion here: https://github.com/orgs/osqp/discussions/681
import numpy as np
import scipy.sparse as spa
from osqp import OSQP
class HuberExample:
def __init__(self, n, seed=1):
np.random.seed(seed)
self.n = int(n) # Number of features
self.m = int(self.n * 100) # Number of data-points
self.Ad = spa.random(self.m, self.n, density=0.15, data_rvs=np.random.randn)
self.x_true = np.random.randn(n) / np.sqrt(n)
ind95 = (np.random.rand(self.m) < 0.95).astype(float)
self.bd = (
self.Ad.dot(self.x_true)
+ np.multiply(0.5 * np.random.randn(self.m), ind95)
+ np.multiply(10.0 * np.random.rand(self.m), 1.0 - ind95)
)
def get_problem(self):
# Construct the problem
# minimize 1/2 z.T * z + np.ones(m).T * (r + s)
# subject to Ax - b - z = r - s
# r >= 0
# s >= 0
# The problem reformulation follows from Eq. (24) of the following paper:
# https://doi.org/10.1109/34.877518
# x_solver = (x, z, r, s)
Im = spa.eye(self.m)
P = spa.block_diag(
(
spa.csc_matrix((self.n, self.n)),
Im,
spa.csc_matrix((2 * self.m, 2 * self.m)),
),
format="csc",
)
q = np.hstack([np.zeros(self.n + self.m), np.ones(2 * self.m)])
A = spa.bmat(
[[self.Ad, -Im, -Im, Im], [None, None, Im, None], [None, None, None, Im]],
format="csc",
)
l = np.hstack([self.bd, np.zeros(2 * self.m)])
u = np.hstack([self.bd, np.inf * np.ones(2 * self.m)])
# Constraints without bounds
A_nobounds = spa.hstack([self.Ad, -Im, -Im, Im], format="csc")
l_nobounds = self.bd
u_nobounds = self.bd
# Bounds
lx = np.hstack([-np.inf * np.ones(self.n + self.m), np.zeros(2 * self.m)])
ux = np.inf * np.ones(self.n + 3 * self.m)
bounds_idx = np.arange(self.n + self.m, self.n + 3 * self.m)
problem = {}
problem["P"] = P
problem["q"] = q
problem["A"] = A
problem["l"] = l
problem["u"] = u
problem["m"] = A.shape[0]
problem["n"] = A.shape[1]
problem["A_nobounds"] = A_nobounds
problem["l_nobounds"] = l_nobounds
problem["u_nobounds"] = u_nobounds
problem["bounds_idx"] = bounds_idx
problem["lx"] = lx
problem["ux"] = ux
return problem
# if __name__ == "__main__":
settings = {
"eps_abs": 0.001,
"eps_rel": 0.001,
"polish": False,
"max_iter": 1000000000,
"eps_prim_inf": 1e-15,
"eps_dual_inf": 1e-15,
"verbose": False,
"time_limit": 1000.0,
"adaptive_rho_tolerance": 5,
"verbose": True
}
problem = HuberExample(200).get_problem()
solvers = [
OSQP(),
OSQP(algebra="cuda")
]
for solver in solvers:
solver.setup(
problem['P'],
problem['q'],
problem['A'],
problem['l'],
problem['u'],
**settings
)
print(f"Solving using {solver}")
solver.solve()
Does anyone know why I am getting this run time error:
Here is my device info:

Metadata
Metadata
Assignees
Labels
No labels