Efficient substitutions of moment equalities if one side of the equal…

…ity is the moment of a monomial, and the other side is a constant

doc/source/revision_history.rst

@@ -1,6 +1,10 @@
 Revision History
 ****************
 
+**Since 1.10.1**
+  - New: Very efficient substitutions of moment equalities if one side of the equality is the moment of a monomial, and the other side is a constant.
+
+
 **Version 1.10.1 (2016-01-29)** 
   - Fixed: The moment equalities are removed correctly if asked.
 

ncpol2sdpa/nc_utils.py

@@ -634,6 +634,33 @@ def find_variable_set(variable_sets, polynomial):
     return -1
 
 
+def check_simple_substitution(equality):
+    if isinstance(equality, str):
+        return (0, 0)
+    elif equality.is_Relational:
+        eq = convert_relational(equality)
+    else:
+        eq = equality
+    if eq.is_Mul or eq.is_Pow:
+        monomial, _ = separate_scalar_factor(eq)
+        return (monomial, 0)
+    if eq.is_Add:
+        components = eq.as_ordered_terms()
+        if len(components) == 2:
+            n_constant = 0
+            constant = 0
+            coeff = 0
+            for component in components:
+                if is_number_type(component):
+                    n_constant += 1
+                    constant = component
+                else:
+                    monomial, coeff = separate_scalar_factor(component)
+            if n_constant == 1:
+                return (monomial, float(-constant/coeff))
+    return (0, 0)
+
+
 def is_number_type(exp):
     return isinstance(exp, (int, float, complex, Number))
 

ncpol2sdpa/sdp_relaxation.py

@@ -18,7 +18,8 @@
 except ImportError:
     from .sparse_utils import lil_matrix
 from .nc_utils import apply_substitutions, build_permutation_matrix, \
-                      convert_relational, find_variable_set, flatten, \
+                      check_simple_substitution, convert_relational, \
+                      find_variable_set, flatten, \
                       get_all_monomials, is_number_type, \
                       is_pure_substitution_rule, iscomplex, ncdegree, \
                       pick_monomials_up_to_degree, save_monomial_index, \
@@ -155,6 +156,7 @@ def __init__(self, variables, parameters=None, verbose=0, normalized=True):
         self.monomial_sets = []
         self.pure_substitution_rules = True
         self.constraints = []
+        self.moment_substitutions = {}
         self.complex_matrix = False
         n_noncommutative_hermitian = 0
         n_noncommutative_nonhermitian = 0
@@ -221,17 +223,22 @@ def _process_monomial(self, monomial, n_vars):
         """
         monomial, coeff = separate_scalar_factor(monomial)
         k = 0
-        # Have we seen this monomial before?
+        # Are we substituting this for a constant?
         try:
-            # If yes, then we improve sparsity by reusing the
-            # previous variable to denote this entry in the matrix
-            k = self.monomial_index[monomial]
+            coeff2 = self.moment_substitutions[monomial]
+            coeff *= coeff2
         except KeyError:
-            # Otherwise we define a new entry in the associated
-            # array recording the monomials, and add an entry in
-            # the moment matrix
-            k = n_vars + 1
-            self.monomial_index[monomial] = k
+            # Have we seen this monomial before?
+            try:
+                # If yes, then we improve sparsity by reusing the
+                # previous variable to denote this entry in the matrix
+                k = self.monomial_index[monomial]
+            except KeyError:
+                # Otherwise we define a new entry in the associated
+                # array recording the monomials, and add an entry in
+                # the moment matrix
+                k = n_vars + 1
+                self.monomial_index[monomial] = k
         return k, coeff
 
     def _push_monomial(self, monomial, n_vars, row_offset, rowA, columnA, N,
@@ -378,16 +385,20 @@ def _get_index_of_monomial(self, element, enablesubstitution=True,
                     monomial = -monomial
                     coeff = -1.0 * coeff
             try:
-                k = self.monomial_index[monomial]
-                result.append((k, coeff))
+                coeff0 = self.moment_substitutions[monomial]
+                result.append((0, coeff0*coeff))
             except KeyError:
-                if not daggered:
-                    dag_result = self._get_index_of_monomial(monomial.adjoint(),
-                                                             daggered=True)
-                    result += [(k, coeff0*coeff) for k, coeff0 in dag_result]
-                else:
-                    raise RuntimeError("The requested monomial " +
-                                       str(monomial) + " could not be found.")
+                try:
+                    k = self.monomial_index[monomial]
+                    result.append((k, coeff))
+                except KeyError:
+                    if not daggered:
+                        dag_result = self._get_index_of_monomial(monomial.adjoint(),
+                                                                 daggered=True)
+                        result += [(k, coeff0*coeff) for k, coeff0 in dag_result]
+                    else:
+                        raise RuntimeError("The requested monomial " +
+                                           str(monomial) + " could not be found.")
         return result
 
     def __push_facvar_sparse(self, polynomial, block_index, row_offset, i, j):
@@ -535,6 +546,8 @@ def __remove_equalities(self, equalities, momentequalities):
         """Attempt to remove equalities by solving the linear equations.
         """
         A = self.__process_equalities(equalities, momentequalities)
+        if A.shape[0] == 0:
+            return
         c = np.array(self.obj_facvar)
         Q, R, P = qr(np.transpose(A[:, 1:]), pivoting=True)
         E = build_permutation_matrix(P)
@@ -738,10 +751,16 @@ def _calculate_block_structure(self, inequalities, equalities, bounds,
             for _ in momentinequalities:
                 monomial_sets.append([S.One])
                 self.block_struct.append(1)
-        if not removeequalities and momentequalities is not None:
-            for _ in momentequalities:
-                monomial_sets += [[S.One], [S.One]]
-                self.block_struct += [1, 1]
+        if momentequalities is not None:
+            for moment_eq in momentequalities:
+                substitution = check_simple_substitution(moment_eq)
+                if substitution == (0, 0):
+                    if not removeequalities:
+                        monomial_sets += [[S.One], [S.One]]
+                        self.block_struct += [1, 1]
+                else:
+                    self.moment_substitutions[substitution[0]] = \
+                        substitution[1]
         self.localizing_monomial_sets = monomial_sets
 
     def __generate_monomial_sets(self, extramonomials):
@@ -859,19 +878,27 @@ def process_constraints(self, inequalities=None, equalities=None,
         if momentinequalities is not None:
             for mineq in momentinequalities:
                 self.constraints.append(mineq)
-        if not removeequalities and momentequalities is not None:
+        reduced_moment_equalities = []
+        if momentequalities is not None:
             for meq in momentequalities:
-                self.constraints.append(meq)
-                if isinstance(meq, str):
-                    tmp = meq.replace("+", "p")
-                    tmp = tmp.replace("-", "+")
-                    tmp = tmp.replace("p", "-")
-                    self.constraints.append(tmp)
-                else:
-                    self.constraints.append(-meq)
+                substitution = check_simple_substitution(meq)
+                if substitution == (0, 0):
+                    if not removeequalities:
+                        self.constraints.append(meq)
+                        if isinstance(meq, str):
+                            tmp = meq.replace("+", "p")
+                            tmp = tmp.replace("-", "+")
+                            tmp = tmp.replace("p", "-")
+                            self.constraints.append(tmp)
+                        else:
+                            self.constraints.append(-meq)
+                    else:
+                        reduced_moment_equalities.append(meq)
         self.__process_inequalities(block_index)
+        if reduced_moment_equalities == []:
+            reduced_moment_equalities = None
         if removeequalities:
-            self.__remove_equalities(equalities, momentequalities)
+            self.__remove_equalities(equalities, reduced_moment_equalities)
 
     def set_objective(self, objective, extraobjexpr=None):
         """Set or change the objective function of the polynomial optimization