Bounds investigation (#33)

* div canon stuff

* some test problems in sandbox

* preparing for merge

* some progress

* weird value check

* stress test ML zero mean passes

* very subtle implementation details. MLE stress test (with and without zero mean) work

* minor

* fixed dimension error inside div-canon and added Sharpe ratio as a test

* fixing some linting issues

* more fixing imports

* fixing final imports, giving up

* added entropy related atoms

* merge fix

* merge fix

* More LICQ-friendly abs formulation

---------

Co-authored-by: Daniel <[email protected]>
Co-authored-by: William Zijie Zhang <[email protected]>

Author: Transurgeon

cvxpy/atoms/elementwise/entr.py

@@ -93,3 +93,6 @@ def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0]
+    
+    def point_in_domain(self):
+        return np.ones(self.shape)

cvxpy/atoms/elementwise/kl_div.py

@@ -94,3 +94,6 @@ def _domain(self) -> List[Constraint]:
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0, self.args[1] >= 0]
+    
+    def point_in_domain(self, argument=0):
+        return np.ones(self.args[argument].shape)

cvxpy/atoms/elementwise/log.py

@@ -99,4 +99,4 @@ def _domain(self) -> List[Constraint]:
     def point_in_domain(self) -> np.ndarray:
         """Returns a point in the domain of the node.
         """
-        return np.ones(self.size)
+        return np.ones(self.shape)

cvxpy/atoms/elementwise/rel_entr.py

@@ -95,3 +95,6 @@ def _domain(self):
         """Returns constraints describing the domain of the node.
         """
         return [self.args[0] >= 0, self.args[1] >= 0]
+
+    def point_in_domain(self, argument=0):
+        return np.ones(self.args[argument].shape)

cvxpy/reductions/expr2smooth/canonicalizers/__init__.py

@@ -15,6 +15,9 @@
 """
 from cvxpy.atoms import maximum
 from cvxpy.atoms.elementwise.log import log
+from cvxpy.atoms.elementwise.entr import entr
+from cvxpy.atoms.elementwise.rel_entr import rel_entr
+from cvxpy.atoms.elementwise.kl_div import kl_div
 from cvxpy.atoms.elementwise.minimum import minimum
 from cvxpy.atoms.elementwise.power import power
 from cvxpy.atoms.pnorm import Pnorm
@@ -27,6 +30,9 @@
 from cvxpy.reductions.expr2smooth.canonicalizers.pnorm_canon import pnorm_canon
 from cvxpy.reductions.expr2smooth.canonicalizers.power_canon import power_canon
 from cvxpy.reductions.expr2smooth.canonicalizers.maximum_canon import maximum_canon
+from cvxpy.reductions.expr2smooth.canonicalizers.entr_canon import entr_canon
+from cvxpy.reductions.expr2smooth.canonicalizers.rel_entr_canon import rel_entr_canon
+from cvxpy.reductions.expr2smooth.canonicalizers.kl_div_canon import kl_div_canon
 
 CANON_METHODS = {
     abs: abs_canon,
@@ -35,5 +41,8 @@
     log: log_canon,
     power: power_canon,
     Pnorm : pnorm_canon,
-    DivExpression: div_canon
+    DivExpression: div_canon,
+    entr: entr_canon,
+    rel_entr: rel_entr_canon,
+    kl_div: kl_div_canon,
 }

cvxpy/reductions/expr2smooth/canonicalizers/abs_canon.py

@@ -16,12 +16,38 @@
 
 import numpy as np
 
+
 from cvxpy.expressions.variable import Variable
+from cvxpy.atoms.elementwise.power import power
+
+# TODO (DCED): ask William if this the multiplication we want to use
+from cvxpy.atoms.affine.binary_operators import multiply
+from cvxpy.reductions.expr2smooth.canonicalizers.power_canon import power_canon
 
+#def abs_canon(expr, args):
+#    shape = expr.shape
+#    t1 = Variable(shape, bounds = [0, None])
+#    if expr.value is not None:
+#       #t1.value = np.sqrt(expr.value**2)
+#        t1.value = np.abs(args[0].value)
+#
+#    #return t1, [t1**2 == args[0] ** 2]
+#    square_expr = power(args[0], 2)
+#    t2, constr_sq = power_canon(square_expr, square_expr.args)
+#    return t1, [t1**2 == t2] + constr_sq
 
 def abs_canon(expr, args):
     shape = expr.shape
-    t = Variable(shape)
-    if expr.value is not None:
-        t.value = np.sqrt(expr.value**2)
-    return t, [t**2 == expr**2, t >= 0]
+    t1 = Variable(shape, bounds = [0, None])
+    y = Variable(shape, bounds = [-1.01, 1.01])
+    if args[0].value is not None:
+        #t1.value = np.sqrt(expr.value**2)
+        t1.value = np.abs(args[0].value)
+        y.value = np.sign(args[0].value)
+    
+    t1.value = np.ones(shape)
+    y.value = np.zeros(shape)
+
+    # TODO (DCED): check how multiply is canonicalized. We don't want to introduce a new variable for
+    # y inside multiply. But args[0] should potentially be canonicalized further?
+    return t1, [y ** 2 == np.ones(shape), t1 == multiply(y, args[0])]

cvxpy/reductions/expr2smooth/canonicalizers/entr_canon.py

@@ -0,0 +1,28 @@
+"""
+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
+import numpy as np
+
+def entr_canon(expr, args):
+    t = Variable(args[0].shape, bounds=[0, None])
+    if args[0].value is not None and np.all(args[0].value >= 1):
+        t.value = args[0].value
+    else:
+        t.value = expr.point_in_domain()
+
+    return expr.copy([t]), [t==args[0]]
+

cvxpy/reductions/expr2smooth/canonicalizers/kl_div_canon.py

@@ -0,0 +1,45 @@
+"""
+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
+import numpy as np
+
+def kl_div_canon(expr, args):
+    constraints = []
+
+    if not args[0].is_constant():
+        t1 = Variable(args[0].shape, bounds=[0, None])
+        constraints.append(t1 == args[0])
+
+        if args[0].value is not None and np.all(args[0].value >= 1):
+            t1.value = args[0].value
+        else:
+            t1.value = expr.point_in_domain(argument=0)
+    else:
+        t1 = args[0]
+
+    if not args[1].is_constant():
+        t2 = Variable(args[1].shape, bounds=[0, None])
+        constraints.append(t2 == args[1])
+
+        if args[1].value is not None and np.all(args[1].value >= 1):
+            t2.value = args[1].value
+        else:
+            t2.value = expr.point_in_domain(argument=1)
+    else:
+        t2 = args[1]
+
+    return expr.copy([t1, t2]), constraints

cvxpy/reductions/expr2smooth/canonicalizers/log_canon.py

@@ -18,7 +18,24 @@
 
 from cvxpy.expressions.constants import Constant
 from cvxpy.expressions.variable import Variable
+from cvxpy.expressions.constants import Constant
+import numpy as np
+from cvxpy.atoms.elementwise.exp import exp
+
+def collect_constant_and_variable(expr, constants, variable):
+    if isinstance(expr, Constant):
+        constants.append(expr)
+    elif isinstance(expr, Variable):
+        variable.append(expr)
+    elif hasattr(expr, "args"):
+        for subexpr in expr.args:
+            collect_constant_and_variable(subexpr, constants, variable)
 
+    assert(len(variable) <= 1)
+
+# DCED: Without this lower bound the stress test for ML Gaussian non-zero mean fails.
+# Perhaps this should be a parameter exposed to the user?
+LOWER_BOUND = 1e-5
 
 def collect_constant_and_variable(expr, constants, variable):
     if isinstance(expr, Constant):
@@ -36,7 +53,7 @@ def collect_constant_and_variable(expr, constants, variable):
 LOWER_BOUND = 1e-5
 
 def log_canon(expr, args):
-    t = Variable(args[0].size, bounds=[LOWER_BOUND, None], name='t')
+    t = Variable(args[0].size, bounds=[LOWER_BOUND, None])
 
     # DCED: if args[0] is a * x for a constant scalar or vector 'a' 
     # and a vector variable 'x', we want to add bounds to x if x
@@ -71,3 +88,15 @@ def log_canon(expr, args):
         t.value = args[0].value
 
     return expr.copy([t]), [t==args[0]]
+
+# TODO (DCED): On some problems this canonicalization seems to work better.
+#              We should investigate this further when we have more benchmarks
+#              involving log.
+#def log_canon(expr, args):
+#    t = Variable(args[0].size)
+#    if args[0].value is not None and np.all(args[0].value > 0):
+#        t.value = np.log(args[0].value)
+#    else:
+#        t.value = expr.point_in_domain()
+
+#    return t, [exp(t) == args[0]]

cvxpy/reductions/expr2smooth/canonicalizers/quad_over_lin.py

@@ -0,0 +1,21 @@
+"""
+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 quad_over_lin_canon(expr, args):
+    assert(False)
+    pass