diff options
| -rw-r--r-- | sumofsquares/canon.py | 45 |
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) |