adding more examples: ML, logistic regression and rayleigh quotient

cvxpy/atoms/affine/binary_operators.py

@@ -190,6 +190,9 @@ def _grad(self, values):
         X = values[0]
         Y = values[1]
 
+        if np.isscalar(X) or np.isscalar(Y):
+            return [sp.csc_matrix(Y), sp.csc_matrix(X)]
+
         # Get dimensions
         m, n = self.args[0].shape if len(self.args[0].shape) == 2 else (self.args[0].size, 1)
         n2, p = self.args[1].shape if len(self.args[1].shape) == 2 else (self.args[1].size, 1)
@@ -198,7 +201,7 @@ def _grad(self, values):
         assert n == n2, f"Inner dimensions must match for multiplication: {n} != {n2}"
 
         # Compute ∂vec(Z)/∂vec(X) = (Y.T ⊗ I_m).T
-        # This is a (m*n) × (m*p) matrix 
+        # This is a (m*n) × (m*p) matrix
         DX = sp.kron(Y.T, sp.eye(m, format='csc'), format='csc').T
 
         # Compute ∂vec(Z)/∂vec(Y) = (I_p ⊗ X).T

cvxpy/reductions/expr2smooth/canonicalizers/div_canon.py

@@ -14,8 +14,10 @@
 limitations under the License.
 """
 
+from cvxpy.atoms.affine.binary_operators import multiply
 from cvxpy.expressions.variable import Variable
 
+
 # We canonicalize div(x, y) as z * y = x.
 def div_canon(expr, args):
     # TODO: potential bounds here?
@@ -33,4 +35,4 @@ def div_canon(expr, args):
     else:
         z.value = expr.point_in_domain()
     # TODO: should we also include y >= 0 here?
-    return z, [z * y == args[0], y == args[1]]#], #y >= 0, z >= 0]
+    return z, [multiply(z, y) == args[0], y == args[1]]#, y >= 0, z >= 0]

cvxpy/reductions/expr2smooth/canonicalizers/log_canon.py

@@ -18,7 +18,7 @@
 
 
 def log_canon(expr, args):
-    t = Variable(args[0].size)
+    t = Variable(args[0].shape)
     if args[0].value is not None:
         t.value = args[0].value
     else:

cvxpy/reductions/expr2smooth/canonicalizers/pnorm_canon.py

@@ -19,17 +19,15 @@
 
 def pnorm_canon(expr, args):
     x = args[0]
-    p = expr.p_rational
-    w = expr.w
+    p = expr.p
 
     if p == 1:
         return x, []
 
     shape = expr.shape
     t = Variable(shape)
     if p % 2 == 0:
-        summation = [x[i]**p for i in range(len(x))]
+        summation = sum(x[i]**p for i in range(x.size))
         return t, [t**p == summation, t >= 0]
     else:
         z = Variable(shape)
-

cvxpy/reductions/expr2smooth/canonicalizers/trig_canon.py

@@ -0,0 +1,41 @@
+"""
+Copyright 2025 CVXPY developers
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+from cvxpy.expressions.variable import Variable
+
+
+def cos_canon(expr, args):
+    if type(args[0]) is not Variable:
+        t = Variable(args[0].size)
+        t.value = expr.point_in_domain()
+        return expr.copy([t]), [t==args[0]]
+    return expr, []
+
+
+def sin_canon(expr, args):
+    if type(args[0]) is not Variable:
+        t = Variable(args[0].size)
+        t.value = expr.point_in_domain()
+        return expr.copy([t]), [t==args[0]]
+    return expr, []
+
+
+def tan_canon(expr, args):
+    if type(args[0]) is not Variable:
+        t = Variable(args[0].size)
+        t.value = expr.point_in_domain()
+        return expr.copy([t]), [t==args[0]]
+    return expr, []

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -14,9 +14,9 @@
 limitations under the License.
 """
 
 import numpy as np
 import scipy.sparse as sp
 from time import time
 
 import cvxpy.settings as s
 from cvxpy.constraints import (
@@ -244,7 +244,7 @@
                         jacobian = grad_dict[var].T
                         if sp.issparse(jacobian):
                             jacobian = sp.dok_matrix(jacobian)
                             data = np.array([jacobian.get((r, c), 0) for r, c in zip(rows, cols)])        
                             values.append(np.atleast_1d(data))
                         else:
                             values.append(np.atleast_1d(jacobian))
@@ -261,7 +261,7 @@
                     #var.value = self.initial_point[offset]
                     var.value = np.nan
                 else:
-                    var.value = np.nan * np.ones(var.size)
+                    var.value = (np.nan * np.ones(var.size)).reshape(var.shape, order='F')
                     #var.value = np.atleast_1d(self.initial_point[offset:offset + var.size])
                 #offset += var.size
             rows, cols = [], []

cvxpy/sandbox/logistic.py

@@ -0,0 +1,50 @@
+import cvxpy as cp
+import numpy as np
+
+# Generate synthetic data
+np.random.seed(42)
+n_samples = 1000
+n_features = 50
+
+# Create two classes
+X_class0 = np.random.randn(n_samples // 2, n_features) - 1
+X_class1 = np.random.randn(n_samples // 2, n_features) + 1
+X = np.vstack([X_class0, X_class1])
+
+# Labels: -1 for class 0, +1 for class 1
+y = np.concatenate([-np.ones(n_samples // 2), np.ones(n_samples // 2)])
+
+# Add intercept term (bias)
+X_with_intercept = np.column_stack([np.ones(n_samples), X])
+
+# Define CVXPY variables
+n_features_with_intercept = n_features + 1
+w = cp.Variable(n_features_with_intercept)  # weights including bias
+
+# Regularization parameter
+lambda_reg = 0.1
+
+# Logistic regression objective: minimize negative log-likelihood + L2 regularization
+# Using log-sum-exp formulation for numerical stability
+log_likelihood = cp.sum(cp.logistic(-cp.multiply(y, X_with_intercept @ w)))
+regularization = lambda_reg * cp.norm(w[1:], 2)**2  # Don't regularize intercept
+
+# Define the optimization problem
+objective = cp.Minimize(log_likelihood + regularization)
+problem = cp.Problem(objective)
+
+# Solve the problem
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+# Print results
+print(f"Optimization status: {problem.status}")
+print(f"Optimal objective value: {problem.value:.4f}")
+print("Optimal weights (including bias):")
+print(f"  Bias (intercept): {w.value[0]:.4f}")
+for i in range(n_features):
+    print(f"  Weight {i+1}: {w.value[i+1]:.4f}")
+
+# Make predictions on training data
+predictions = np.sign(X_with_intercept @ w.value)
+accuracy = np.mean(predictions == y)
+print(f"\nTraining accuracy: {accuracy:.2%}")

cvxpy/sandbox/mle-canon.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+import cvxpy as cp
+
+n = 2
+np.random.seed(1234)
+data = np.random.randn(n)
+
+mu = cp.Variable(1, name="mu")
+mu.value = np.array([0.0])
+sigma = cp.Variable(1, name="sigma")
+sigma.value = np.array([1.0])
+t1 = cp.Variable(1, name="t1")
+t2 = cp.Variable(n, name="t2")
+t3 = cp.Variable(1, name="t3")
+t4 = cp.Variable(1, name="t4")
+t5 = cp.Variable(1, name="t5")
+
+log_likelihood = ((n / 2) * cp.log(t4) - cp.sum(cp.square(t2)) / (2 * (t1)**2))
+
+constraints = [
+    t1 == sigma,
+    t2 == data - mu,
+    t3 == (2 * np.pi * (t5)**2),
+    t4 == 1 / t3,
+    t5 == sigma,
+]
+t1.value = sigma.value
+t2.value = data - mu.value
+t3.value = (2 * np.pi * (sigma.value)**2)
+t4.value = 1 / t3.value
+t5.value = sigma.value
+
+objective = cp.Maximize(log_likelihood)
+problem = cp.Problem(objective, constraints)
+problem.solve(solver=cp.IPOPT, nlp=True)
+assert problem.status == cp.OPTIMAL
+assert np.allclose(mu.value, np.mean(data))
+assert np.allclose(sigma.value, np.std(data))