making analytic center of ellipsoid work

Author: Transurgeon

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -16,7 +16,6 @@
 
 import numpy as np
 import scipy.sparse as sp
-import torch
 
 import cvxpy.settings as s
 from cvxpy.constraints import (
@@ -32,7 +31,6 @@
     nonpos2nonneg,
 )
 from cvxpy.utilities.citations import CITATION_DICT
-from cvxtorch import TorchExpression
 
 
 class IPOPT(NLPsolver):
@@ -146,7 +144,7 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         nlp = cyipopt.Problem(
         n=len(x0),
         m=len(bounds.cl),
-        problem_obj=self.Oracles(bounds.new_problem),
+        problem_obj=self.Oracles(bounds.new_problem, x0),
         lb=bounds.lb,
         ub=bounds.ub,
         cl=bounds.cl,
@@ -169,9 +167,10 @@ def cite(self, data):
         return CITATION_DICT["IPOPT"]
 
     class Oracles():
-        def __init__(self, problem):
+        def __init__(self, problem, inital_point):
             self.problem = problem
             self.main_var = []
+            self.initial_point = inital_point
             for var in self.problem.variables():
                 self.main_var.append(var)
 
@@ -201,7 +200,10 @@ def gradient(self, x):
             for var in self.main_var:
                 size = var.size
                 if var in grad_dict:
-                    grad[grad_offset:grad_offset+size] = grad_dict[var]
+                    array = grad_dict[var]
+                    if sp.issparse(array):
+                        array = array.toarray().flatten(order='F')
+                    grad[grad_offset:grad_offset+size] = array
                 grad_offset += size
             return grad
 
@@ -231,23 +233,24 @@ def jacobian(self, x):
 
             values = []
             for constraint in self.problem.constraints:
+                # get the jacobian of the constraint
                 grad_dict = constraint.expr.grad
                 for var in self.main_var:
                     if var in grad_dict:
-                        array = grad_dict[var]
-                        if sp.issparse(array):
-                            array = array.toarray().flatten(order='F')
-                            values.append(np.atleast_1d(array))
+                        jacobian = grad_dict[var].T
+                        if sp.issparse(jacobian):
+                            jacobian = jacobian.toarray().flatten(order='F')
+                            values.append(np.atleast_1d(jacobian))
                         else:
-                            values.append(np.atleast_1d(array))
+                            values.append(np.atleast_1d(jacobian))
             return np.concatenate(values)
 
         def jacobianstructure(self):
             """Returns the sparsity structure of the Jacobian."""
             # Set dummy values to get gradient structure
             offset = 0
             for var in self.main_var:
-                var.value = np.zeros(var.shape)
+                var.value = self.initial_point[offset:offset + var.size]
                 offset += var.size
             rows, cols = [], []
             row_offset = 0
@@ -256,12 +259,16 @@ def jacobianstructure(self):
                 col_offset = 0
                 for var in self.main_var:
                     if var in grad_dict:
-                        array = grad_dict[var]
-                        if sp.issparse(array):
-                            array = array.toarray()
-                        # get row and column indices
-                        rows.extend(np.ones(array.size)*row_offset)
-                        cols.extend(np.arange(col_offset, col_offset + var.size))
+                        jacobian = grad_dict[var].T
+                        if sp.issparse(jacobian):
+                            row_indices, col_indices = np.indices(jacobian.shape)
+                            rows.extend(row_indices.flatten(order='F') + row_offset)
+                            cols.extend(col_indices.flatten(order='F') + col_offset)
+                        """
+                        else:
+                            rows.extend(np.ones(jacobian.size)*row_offset)
+                            cols.extend(np.arange(col_offset, col_offset + var.size))
+                        """
                     col_offset += var.size
                 row_offset += constraint.size
 

cvxpy/sandbox/sum-log.py

@@ -5,14 +5,25 @@
 
 # Generate random data
 np.random.seed(0)
-m, n = 1, 2
-A = np.random.randn(m, n)
-b = np.random.randn(m)
+m, n = 500, 40
+b = np.ones(m)
+rand = np.random.randn(m - 2*n, n)
+A = np.vstack((rand, np.eye(n), np.eye(n) * -1))
+"""
+m, n = 5, 2
+A = np.array([[1, 0], [-1, 0], [0, 1], [0, -1], [-0.5, 1]])
+b = np.array([1, 1, 1, 1, 0.5])
+"""
 # Define the variable
 x = cp.Variable(n)
 # set initial value for x
-objective = cp.Minimize(-cp.sum(cp.log(A @ x - b)))
+objective = cp.Minimize(-cp.sum(cp.log(-A @ x + b)))
 problem = cp.Problem(objective, [])
 # Solve the problem
 problem.solve(solver=cp.IPOPT, nlp=True)
 print("Optimal value:", problem.value)
+print("Optimal x:", x.value)
+
+problem.solve(solver=cp.CLARABEL)
+print("Optimal value:", problem.value)
+print("Optimal x:", x.value)