more fixes for geo_mean canonicalization #1685 (follow-up on PR #1687)

cvxpy/atoms/geo_mean.py

@@ -21,6 +21,7 @@
 
 from cvxpy.atoms.atom import Atom
 from cvxpy.constraints.constraint import Constraint
+from cvxpy.expressions import cvxtypes
 from cvxpy.utilities.power_tools import (approx_error, decompose, fracify,
                                          lower_bound, over_bound, prettydict,)
 
@@ -199,6 +200,12 @@ def __init__(self, x, p: Optional[List[int]] = None, max_denom: int = 1024) -> N
             The number of second order cones used to form this geometric mean
 
         """
+        if p is not None and hasattr(p, '__getitem__'):
+            p = np.array(p)
+            idxs = p > 0
+            Expression = cvxtypes.expression()
+            x = Expression.cast_to_const(x)[idxs]
+            p = p[idxs]
         super(geo_mean, self).__init__(x)
 
         x = self.args[0]

cvxpy/tests/test_power_tools.py

@@ -0,0 +1,113 @@
+"""
+Copyright 2022, the 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.
+"""
+
+import numpy as np
+
+import cvxpy as cp
+from cvxpy.tests.base_test import BaseTest
+
+
+class TestGeoMean(BaseTest):
+
+    def test_multi_step_dyad_completion(self) -> None:
+        """
+        Consider four market equilibrium problems.
+
+        The budgets "b" in these problems are chosen so that canonicalization
+        of geo_mean(u, b) hits a recursive code-path in power_tools.dyad_completion(...).
+
+        The reference solution is computed by taking the log of the geo_mean objective,
+        which has the effect of making the problem ExpCone representable.
+        """
+        if 'MOSEK' in cp.installed_solvers():
+            log_solve_args = {'solver': 'MOSEK'}
+        else:
+            log_solve_args = {'solver': 'ECOS'}
+        n_buyer = 5
+        n_items = 7
+        np.random.seed(0)
+        V = 0.5 * (1 + np.random.rand(n_buyer, n_items))
+        X = cp.Variable(shape=(n_buyer, n_items), nonneg=True)
+        cons = [cp.sum(X, axis=0) <= 1]
+        u = cp.sum(cp.multiply(V, X), axis=1)
+        bs = np.array([
+            [110, 14, 6, 77, 108],
+            [15., 4., 8., 0., 9.],
+            [14., 21., 217., 57., 6.],
+            [3., 36., 77., 8., 8.]
+        ])
+        for i, b in enumerate(bs):
+            log_objective = cp.Maximize(b @ cp.log(u))
+            log_prob = cp.Problem(log_objective, cons)
+            log_prob.solve(**log_solve_args)
+            expect_X = X.value
+
+            geo_objective = cp.Maximize(cp.geo_mean(u, b))
+            geo_prob = cp.Problem(geo_objective, cons)
+            geo_prob.solve()
+            actual_X = X.value
+            try:
+                self.assertItemsAlmostEqual(actual_X, expect_X, places=3)
+            except AssertionError as e:
+                print(f'Failure at index {i} (when b={str(b)}).')
+                log_prob.solve(**log_solve_args, verbose=True)
+                print(X.value)
+                geo_prob.solve(verbose=True)
+                print(X.value)
+                print('The valuation matrix was')
+                print(V)
+                raise e
+
+    def test_3d_power_cone_approx(self):
+        """
+        Use
+            geo_mean((x,y), (alpha, 1-alpha)) >= |z|
+        as a reformulation of
+            PowCone3D(x, y, z, alpha).
+
+        Check validity of the reformulation by solving
+        orthogonal projection problems.
+        """
+        if 'MOSEK' in cp.installed_solvers():
+            proj_solve_args = {'solver': 'MOSEK'}
+        else:
+            proj_solve_args = {'solver': 'SCS', 'eps': 1e-10}
+        min_numerator = 2
+        denominator = 25
+        x = cp.Variable(3)
+        np.random.seed(0)
+        y = 10 * np.random.rand(3)  # the third value doesn't matter
+        for i, numerator in enumerate(range(min_numerator, denominator, 3)):
+            alpha_float = numerator / denominator
+            y[2] = (y[0] ** alpha_float) * (y[1] ** (1 - alpha_float)) + 0.05
+            objective = cp.Minimize(cp.norm(y - x, 2))
+
+            actual_constraints = [cp.constraints.PowCone3D(x[0], x[1], x[2],
+                                                           [alpha_float])]
+            actual_prob = cp.Problem(objective, actual_constraints)
+            actual_prob.solve(**proj_solve_args)
+            actual_x = x.value.copy()
+
+            weights = np.array([alpha_float, 1 - alpha_float])
+            approx_constraints = [cp.geo_mean(x[:2], weights) >= cp.abs(x[2])]
+            approx_prob = cp.Problem(objective, approx_constraints)
+            approx_prob.solve()
+            approx_x = x.value.copy()
+            try:
+                self.assertItemsAlmostEqual(actual_x, approx_x, places=4)
+            except AssertionError as e:
+                print(f'Failure at index {i} (when alpha={alpha_float}).')
+                raise e

cvxpy/tests/test_problem.py

@@ -1640,20 +1640,6 @@ def short_geo_mean(x, p):
         xval = np.array(x.value).flatten()
         self.assertTrue(np.allclose(xval, x_true, 1e-3))
 
-        # GH Issue #1685 (first try Fractions, then try floats)
-        p = np.array([Fraction(1, 8), Fraction(1, 6), Fraction(1, 12),
-                      Fraction(3, 16), Fraction(7, 16)])
-        expr = cp.geo_mean(x, p)
-        prob = cp.Problem(cp.Maximize(cp.min(x)), [0 <= x, expr >= 1, x <= 1])
-        prob.solve()
-        self.assertItemsAlmostEqual(x.value, x_true)  # using Fractions
-
-        p = np.array([1/8, 1/6, 1/12, 3/16, 7/16])
-        expr = cp.geo_mean(x, p)
-        prob = cp.Problem(cp.Maximize(cp.min(x)), [0 <= x, expr >= 1, x <= 1])
-        prob.solve()
-        self.assertItemsAlmostEqual(x.value, x_true)  # using floats
-
     def test_pnorm(self) -> None:
         import numpy as np
 

cvxpy/utilities/power_tools.py

@@ -66,8 +66,9 @@ def gm_constrs(t, x_list, p):
     d = defaultdict(lambda: Variable(t.shape))
     d[w] = t
 
-    if len(x_list) < len(w):
-        x_list += [t]
+    long_w = len(w) - len(x_list)
+    if long_w > 0:
+        x_list += [t]*long_w
 
     assert len(x_list) == len(w)
 
@@ -225,6 +226,14 @@ def is_weight(w) -> bool:
     return valid_elems and sum(w) == 1
 
 
+__EXCEED_DENOMINATOR_LIMIT__ = \
+    """
+    Can't reliably represent the input weight vector.
+    Try increasing `max_denom` or checking the denominators
+    of your input fractions.
+    """
+
+
 def fracify(a, max_denom: int = 1024, force_dyad: bool = False):
     """ Return a valid fractional weight tuple (and its dyadic completion)
         to represent the weights given by ``a``.
@@ -263,6 +272,7 @@ def fracify(a, max_denom: int = 1024, force_dyad: bool = False):
         That is, if w has fractions with denominators that are not a power of 2,
         and ``len(w) == n`` then w_dyad has length n+1, dyadic fractions for elements,
         and ``w_dyad[:-1]/w_dyad[n] == w``.
+        # ^ That isn't always possible, is it?
 
         Alternatively, the ratios between the
         first n elements of ``w_dyad`` are equal to the corresponding ratios between
@@ -362,17 +372,18 @@ def fracify(a, max_denom: int = 1024, force_dyad: bool = False):
         w_frac = tuple(Fraction(v, total) for v in a)
         d = max(v.denominator for v in w_frac)
         if d > max_denom:
-            msg = ("Can't reliably represent the input weight vector."
-                   "\nTry increasing `max_denom` or checking the denominators "
-                   "of your input fractions.")
-            raise ValueError(msg)
+            raise ValueError(__EXCEED_DENOMINATOR_LIMIT__)
     else:
         # fall through code
         w_frac = tuple(Fraction(float(v)/total).limit_denominator(max_denom) for v in a)
         if sum(w_frac) != 1:
             w_frac = make_frac(a, max_denom)
 
-    return w_frac, dyad_completion(w_frac)
+    w_dyad = dyad_completion(w_frac)
+    if max(v.denominator for v in w_dyad) > max_denom:
+        raise ValueError(__EXCEED_DENOMINATOR_LIMIT__)
+
+    return w_frac, w_dyad
 
 
 def make_frac(a, denom):
@@ -427,7 +438,8 @@ def dyad_completion(w):
         # need to add the dummy variable to represent as dyadic
         d = max(non_dyad_dens)
         p = next_pow2(d)
-        return tuple(Fraction(v*d, p) for v in w) + (Fraction(p-d, p),)
+        w_aug = tuple(Fraction(v*d, p) for v in w) + (Fraction(p-d, p),)
+        return dyad_completion(w_aug)
     else:
         return w