diff options
Diffstat (limited to '')
| -rw-r--r-- | sumofsquares/canon.py | 81 |
1 files 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 + |