From 3fe39825bf36cbad77bd563c32da257838d7aff3 Mon Sep 17 00:00:00 2001
From: Nao Pross <np@0hm.ch>
Date: Sat, 25 May 2024 12:25:13 +0200
Subject: Sketch (non-working) LogDet canonicalization

---
 sumofsquares/canon.py | 81 ++++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 73 insertions(+), 8 deletions(-)

diff --git a/sumofsquares/canon.py b/sumofsquares/canon.py
index e5960e7..72cab7a 100644
--- a/sumofsquares/canon.py
+++ b/sumofsquares/canon.py
@@ -14,12 +14,6 @@ from polymatrix.statemonad import init_state_monad
 from polymatrix.polymatrix.mixins import PolyMatrixMixin
 from polymatrix.variable import Variable
 
-from polymatrix.expression.init import (
-        init_concatenate_expr,
-        init_diag_expr,
-        init_lower_triangular_expr,
-        init_transpose_expr)
-
 from .abc import Constraint, Result
 from .constraints import NonNegative, PositiveSemiDefinite, ExponentialCone
 from .error import AlgebraicError
@@ -29,7 +23,9 @@ from .utils import partition
 
 
 class Canonicalization:
-    """ Transform a problem into another problem. """
+    """
+    Transform a problem into another problem.
+    """
 
     def __str__(self):
         s = f"Apply canonicalization procedure {self.__class__.__qualname__} to\n"
@@ -54,6 +50,12 @@ class Canonicalization:
 
 @dataclassabc(frozen=True)
 class PutinarPSatz(Canonicalization):
+    """
+    Apply Putinar's Positivstellensatz to replace the domain-restricted
+    non-negativity constraints of a problem, with equivalent global
+    non-negativity constraints by adding langrange-multiplier-type polynomials.
+    """
+
     problem: SOSProblem
 
     @staticmethod
@@ -116,6 +118,7 @@ class PutinarPSatz(Canonicalization):
         for i, constr in enumerate(to_process):
             new_constr = constr.expression
             for domain_poly in constr.domain.polynomials:
+                # FIXME: need to check that there is not a \gamma_{i}!
                 m = PutinarPSatz.make_multiplier(prob.variables,
                                                  prob.polynomial_variables,
                                                  constr, domain_poly, rf"\gamma_{i}")
@@ -125,7 +128,69 @@ class PutinarPSatz(Canonicalization):
                 new_constraints.append(NonNegative(m))
 
             new_constraints.append(NonNegative(new_constr))
-
         return replace(self.problem, constraints=new_constraints).apply(state)
 
 
+@dataclassabc(frozen=True)
+class LogDet(Canonicalization):
+    """
+    The problem
+    
+      maximize logdet(A)
+    
+    is equivalent to solving
+    
+      maximize      t
+    
+      subject to    [ A    Z       ]
+                    [ Z.T  diag(Z) ] >= 0
+    
+                    Z lower triangular
+                    t <= sum_i log(Z[i,i])
+    
+    and the last constraint of the above is equivalent to
+    
+      t <= sum_i u[i]
+      u_i <= log(Z[i, i])  for all i
+    
+    And finally to get rid of the log the latter constraint one is
+    equivalent to
+    
+      (Z[i,i], 1, u[i]) in Exponential Cone for all i
+    
+    Hence we can replace the original problem with
+    
+      minimize      - sum_i u[i]
+    
+      subject to    [ A    Z       ]
+                    [ Z.T  diag(Z) ] >= 0
+    
+                    Z lower triangular
+                    (Z[i,i], 1, u[i]) in ExpCone for all i
+    """
+    problem: SOSProblem
+
+    @override
+    def apply(self, state: ExpressionState) -> tuple[ExpressionState, InternalSOSProblem]:
+        A = self.problem.cost
+
+        # FIXME: should check for name clashes in state object
+        # i.e. what if there is already a variable named u_logdet / Z_logdet?
+        # TODO: deterministically generate unique names
+
+        n = A.shape[0]
+        u = opt_variable_from_name('u_logdet', shape=(n, 1))
+
+        m = poly.v_stack((n * (n + 1) // 2, 1))
+        Z = poly.lower_triangular(opt_variable_from_name('Z_logdet', shape=m))
+
+        # we call the new big matrix Q
+        Q = poly.concatenate(((A, Z), (Z.T, Z.diag())))
+
+        E = poly.h_stack((Z.diag(), poly.ones((n, 1)), u))
+
+        new_cost = - u.T.sum()
+        new_constraints = [PositiveSemiDefinite(Q), ExponentialCone(E)]
+
+        raise NotImplementedError
+
-- 
cgit v1.2.1