Fix bug in LogDet canonmain

1 files changed, 24 insertions, 21 deletions

diff --git a/sumofsquares/canon.py b/sumofsquares/canon.py
index 97a271b..f315de1 100644
--- a/sumofsquares/canon.py
+++ b/sumofsquares/canon.py
@@ -3,7 +3,7 @@ import polymatrix as poly
 
 from abc import abstractmethod
 from dataclassabc import dataclassabc
-from dataclasses import replace
+from dataclasses import replace, field
 from typing import Iterable
 from typing_extensions import override
 
@@ -81,7 +81,7 @@ class PutinarPSatz(Canonicalization):
         constraints, to_process = partition(need_positivstellensatz, self.problem.constraints)
         constraints = tuple(constraints)
         to_process = tuple(to_process)
-
+ 
         # Variables that are involved in the problem
         variable_indices: set[VariableIndex] = set()
 
@@ -192,31 +192,36 @@ class LogDet(Canonicalization):
     
       maximize      t
     
-      subject to    [ A    Z       ]
-                    [ Z.T  diag(Z) ] >= 0
+      subject to    [ A    L       ]
+                    [ L.T  diag(L) ] >= 0
     
-                    Z lower triangular
-                    t <= sum_i log(Z[i,i])
+                    L lower triangular
+                    t <= sum_i log(L[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
+      u_i <= log(L[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
+      exp(u_i) <= L[i,i]
+
+    so
+
+      (u[i], 1, L[i,i]) in Exponential Cone for all i
     
+    using the definition of the exponential cone of SCS.
     Hence we can replace the original problem with
     
-      minimize      - sum_i u[i]
+      maximize      sum_i u[i]
     
-      subject to    [ A    Z       ]
-                    [ Z.T  diag(Z) ] >= 0
+      subject to    [ A    L       ]
+                    [ L.T  diag(L) ] >= 0
     
-                    Z lower triangular
-                    (Z[i,i], 1, u[i]) in ExpCone for all i
+                    L lower triangular
+                    (u[i], 1, L[i,i]) in ExpCone for all i
     """
     problem: SOSProblem
     cost: Expression | None = None
@@ -227,26 +232,24 @@ class LogDet(Canonicalization):
         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?
+        # i.e. what if there is already a variable named u_logdet / L_logdet?
         # TODO: deterministically generate unique names
 
-        t = opt_variable_from_name('t_logdet')
-
         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))
+        L = poly.lower_triangular(opt_variable_from_name('L_logdet', shape=m))
+
+        Q = poly.concatenate(((A, L.T), (L, L.diag())))
+        E = poly.h_stack((u, poly.ones((n, 1)), L.diag()))
 
-        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_cost = - t
         new_constraints = tuple(self.problem.constraints) + (
             PositiveSemiDefinite(Q),
             ExponentialCone(E),
             NonNegative(u),
-            NonNegative(u.T.sum() - t)
         )
 
         return replace(self.problem, cost=new_cost, constraints=new_constraints).apply(state)