Commit 15d8b60e authored by Kerry He's avatar Kerry He
Browse files

Added direct support for KLdiv constraint

parent d48b4c4a

picos/constraints/con_kldiv.py

+4 −0
Original line number Diff line number Diff line
@@ -123,6 +123,10 @@ class KullbackLeiblerConstraint(Constraint): 
    def _str(self):

        return glyphs.le(self.divergence.string, self.upperBound.string)



    def _get_size(self):

        n = len(self.numerator)

        return (1 + 2*n, 1)



    def _get_slack(self):

        return self.upperBound.safe_value - self.divergence.safe_value

picos/solvers/solver_qics.py

+39 −2
Original line number Diff line number Diff line
@@ -25,7 +25,7 @@ import numpy 


from ..apidoc import api_end, api_start

from ..constraints import (AffineConstraint, DummyConstraint, RSOCConstraint, SOCConstraint,

                           LMIConstraint, ComplexLMIConstraint, ExpConeConstraint)

                           LMIConstraint, ComplexLMIConstraint, ExpConeConstraint, KullbackLeiblerConstraint)

from ..expressions import CONTINUOUS_VARTYPES, AffineExpression

from ..modeling.footprint import Specification

from ..modeling.solution import (PS_FEASIBLE, PS_INFEASIBLE, PS_UNBOUNDED,
@@ -52,7 +52,8 @@ class QICSSolver(Solver): 
            RSOCConstraint,

            LMIConstraint,

            ComplexLMIConstraint,

            ExpConeConstraint])

            ExpConeConstraint,

            KullbackLeiblerConstraint])



    @classmethod

    def supports(cls, footprint, explain=False):
@@ -255,6 +256,30 @@ class QICSSolver(Solver): 
        

        self.int["cones"] += [qics.cones.ClassEntr(1)]

        

    def _import_kldiv_constraint(self, constraint):

        assert isinstance(constraint, KullbackLeiblerConstraint)



        (Gt, ht) = self._Gh(constraint.upperBound)

        (Gx, hx) = self._Gh(constraint.numerator)

        (Gy, hy) = self._Gh(constraint.denominator)

        

        self._qicsConIndices[constraint] = len(self.int["cones"])

        

        # Check if we can reduce to entropy

        if (hy == hy[0]).all:

            Gy_dense = Gy.toarray()

            if (Gy_dense == Gy_dense[0]).all():

                self.int["G"] = self.stack(self.int["G"], -Gt, -Gy[0, :], -Gx)

                self.int["h"] = self.stack(self.int["h"],  ht,  hy[[0]],   hx)

                

                self.int["cones"] += [qics.cones.ClassEntr(len(hx))]

                return

                

        self.int["G"] = self.stack(self.int["G"], -Gt, -Gx, -Gy)

        self.int["h"] = self.stack(self.int["h"],  ht,  hx,  hy)

        

        self.int["cones"] += [qics.cones.ClassRelEntr(len(hx))]



    def _import_lmi_constraint(self, constraint):

        assert isinstance(constraint, LMIConstraint)

        iscomplex = isinstance(constraint, ComplexLMIConstraint)
@@ -301,6 +326,8 @@ class QICSSolver(Solver): 
            self._import_lmi_constraint(constraint)

        elif isinstance(constraint, ExpConeConstraint):

            self._import_expcone_constraint(constraint)

        elif isinstance(constraint, KullbackLeiblerConstraint):

            self._import_kldiv_constraint(constraint)

        else:

            assert isinstance(constraint, DummyConstraint), \

                "Unexpected constraint type: {}".format(
@@ -420,6 +447,16 @@ class QICSSolver(Solver): 
                elif isinstance(constraint, ExpConeConstraint):

                    zxy  = result["z_opt"][indices]

                    dual = cvxopt.matrix([zxy[1][0, 0], zxy[2][0, 0], -zxy[0][0, 0]])

                elif isinstance(constraint, KullbackLeiblerConstraint):

                    dual = result["z_opt"][indices]

                    if len(dual[1]) == 1:

                        # CRE was cast as a CE cone, so transform duals back to CRE

                        t = dual[0]

                        x = dual[2]

                        y = dual[1] * numpy.ones_like(dual[2]) / dual[2].size

                        dual = cvxopt.matrix(numpy.vstack((t, x, y)).ravel())

                    else:

                        dual = cvxopt.matrix(numpy.vstack(dual).ravel())



                duals[constraint] = dual