Coverage for picos/constraints/uncertain/ucon_ws_sqnorm.py: 96.10%
77 statements
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
« prev ^ index » next coverage.py v6.5.0, created at 2023-03-26 07:46 +0000
1# ------------------------------------------------------------------------------
2# Copyright (C) 2020 Maximilian Stahlberg
3#
4# This file is part of PICOS.
5#
6# PICOS is free software: you can redistribute it and/or modify it under the
7# terms of the GNU General Public License as published by the Free Software
8# Foundation, either version 3 of the License, or (at your option) any later
9# version.
10#
11# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY
12# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14#
15# You should have received a copy of the GNU General Public License along with
16# this program. If not, see <http://www.gnu.org/licenses/>.
17# ------------------------------------------------------------------------------
19"""Implements :class:`WassersteinAmbiguousSquaredNormConstraint`."""
21from collections import namedtuple
23from ... import glyphs
24from ...apidoc import api_end, api_start
25from ..constraint import Constraint, ConstraintConversion
27_API_START = api_start(globals())
28# -------------------------------
31class WassersteinAmbiguousSquaredNormConstraint(Constraint):
32 """A bound on a Wasserstein-ambiguous expected value of a squared norm."""
34 class DistributionallyRobustConversion(ConstraintConversion):
35 """Distributionally robust counterpart conversion."""
37 @classmethod
38 def predict(cls, subtype, options):
39 """Implement :meth:`~.constraint.ConstraintConversion.predict`."""
40 from ...expressions import RealVariable, SymmetricVariable
41 from .. import AffineConstraint, LMIConstraint
43 k = subtype.sqnorm_argdim
44 m = subtype.universe_subtype.sample_dim
45 N = subtype.universe_subtype.sample_num
47 yield ("var", RealVariable.make_var_type(dim=1, bnd=1), 1) # gamma
48 yield ("var", RealVariable.make_var_type(dim=N, bnd=0), 1) # s
49 yield ("var", SymmetricVariable.make_var_type( # U
50 dim=(m * (m + 1)) // 2, bnd=0), 1)
51 yield ("var", RealVariable.make_var_type(dim=m, bnd=0), 1) # u
52 yield ("var", RealVariable.make_var_type(dim=1, bnd=0), 1) # mu
54 yield ("con", AffineConstraint.make_type(dim=1, eq=False), 1)
55 yield ("con", LMIConstraint.make_type(diag=(k + m + 1)), 1)
56 yield ("con", LMIConstraint.make_type(diag=(m + 1)), N)
58 @classmethod
59 def convert(cls, con, options):
60 """Implement :meth:`~.constraint.ConstraintConversion.convert`."""
61 # The recipe is found in "Robust conic optimization in
62 # Python" (Stahlberg 2020) and extends a result in "Wasserstein
63 # distributionally robust optimization: Theory and applications in
64 # machine learning" (Kuhn, Esfahani, Nguyen and Shafieezadeh-Abadeh
65 # 2019).
66 from ...expressions import RealVariable, SymmetricVariable
67 from ...expressions.algebra import block
68 from ...modeling import Problem
70 problem = Problem()
72 # Load the uncertain suqared norm.
73 a = con.sqnorm.x
74 B, b = a.factor_out(a.perturbation)
76 # Load the ambiguity set.
77 WAS = con.sqnorm.universe
78 S = WAS.samples
79 m = S.dim
80 N = S.num
81 w = WAS.weights
82 eps = WAS.eps
84 # Load the upper bound.
85 omega = con.ub
87 # Introduce auxiliary variables.
88 gamma = RealVariable("__gamma", lower=0)
89 s = RealVariable("__s", N)
90 U = SymmetricVariable("__U", m)
91 u = RealVariable("__u", m)
92 mu = RealVariable("__mu")
94 # Compute redundant terms that appear in constraints.
95 h1 = gamma.dupdiag(m) - U
96 h2 = tuple(gamma*S[i] + u for i in range(N))
97 h3 = tuple(gamma*abs(S[i])**2 + s[i] - mu for i in range(N))
99 # Add constraints.
100 problem.add_constraint(
101 gamma*eps**2 + w.T*s <= omega)
102 problem.add_constraint(
103 block([["I", B, b],
104 [B.T, U, u],
105 [b.T, u.T, mu]]) >> 0)
106 problem.add_list_of_constraints([
107 block([[h1, h2[i]],
108 [h2[i].T, h3[i]]]) >> 0
109 for i in range(N)])
111 return problem
113 def __init__(self, sqnorm, upper_bound):
114 """Construct a :class:`WassersteinAmbiguousSquaredNormConstraint`.
116 :param ~picos.expressions.UncertainSquaredNorm sqnorm:
117 Uncertain squared norm to upper bound the expectation of.
119 :param ~picos.expressions.AffineExpression upper_bound:
120 Upper bound on the expected value.
121 """
122 from ...expressions import AffineExpression, UncertainSquaredNorm
123 from ...expressions.uncertain.pert_wasserstein import (
124 WassersteinAmbiguitySet)
126 assert isinstance(sqnorm, UncertainSquaredNorm)
127 assert isinstance(sqnorm.universe, WassersteinAmbiguitySet)
128 assert sqnorm.universe.p == 2
129 assert isinstance(upper_bound, AffineExpression)
130 assert upper_bound.scalar
132 self.sqnorm = sqnorm
133 self.ub = upper_bound
135 super(WassersteinAmbiguousSquaredNormConstraint, self).__init__(
136 "Wasserstein-ambiguous Expected Squared Norm", printSize=True)
138 Subtype = namedtuple("Subtype", ("sqnorm_argdim", "universe_subtype"))
140 def _subtype(self):
141 return self.Subtype(len(self.sqnorm.x), self.sqnorm.universe.subtype)
143 @classmethod
144 def _cost(cls, subtype):
145 return float("inf")
147 def _expression_names(self):
148 yield "sqnorm"
149 yield "ub"
151 def _str(self):
152 return glyphs.le(self.sqnorm.worst_case_string("max"), self.ub.string)
154 def _get_size(self):
155 return (1, 1)
157 def _get_slack(self):
158 return self.ub.value - self.sqnorm.worst_case_value(direction="max")
161# --------------------------------------
162__all__ = api_end(_API_START, globals())