Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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# ------------------------------------------------------------------------------ 

18 

19"""Implements :class:`WassersteinAmbiguitySet`.""" 

20 

21from collections import namedtuple 

22 

23import numpy 

24 

25from ... import glyphs 

26from ...apidoc import api_end, api_start 

27from ..data import cvx2np 

28from ..exp_affine import AffineExpression 

29from ..samples import Samples 

30from .perturbation import Perturbation, PerturbationUniverse 

31 

32_API_START = api_start(globals()) 

33# ------------------------------- 

34 

35 

36class WassersteinAmbiguitySet(PerturbationUniverse): 

37 r"""A wasserstein ambiguity set centered at a discrete distribution. 

38 

39 :Model of uncertainty: 

40 

41 As a distributional ambiguity set, an instance :math:`\mathcal{P}` of this 

42 class 

43 

44 1. represents a safety region for a partially known (ambiguous) probability 

45 distribution :math:`\Xi \in \mathcal{P}` and 

46 2. provides a random, ambiguously distributed perturbation parameter 

47 :math:`\xi \sim \Xi` that can be used to define worst-case-expectation 

48 expressions of the form 

49 

50 .. math:: 

51 

52 \mathop{(\max\;\textit{or}\;\min)}_{\Xi \in \mathcal{P}} 

53 \mathbb{E}_\Xi[f(x, \xi)] 

54 

55 for a selection of functions :math:`f` and a decision variable :math:`x`. 

56 

57 :Definition: 

58 

59 Formally, this class can describe discrepancy-based ambiguity sets of the 

60 form 

61 

62 .. math:: 

63 

64 \mathcal{P} = \left\{ 

65 \Xi \in \mathcal{M} ~\middle|~ 

66 \operatorname{W}_p(\Xi, \Xi_\text{N}) \leq \epsilon 

67 \right\} 

68 

69 where discrepancy from the discrete nominal distribution 

70 

71 .. math:: 

72 

73 \Xi_\text{N} = \sum_{i = 1}^N w_i \delta_{\xi_{(i)}} \in \mathcal{M} 

74 

75 is measured with respect to the Wasserstein distance of order 

76 :math:`p \geq 1`, 

77 

78 .. math:: 

79 

80 \operatorname{W}_p(\Xi, \Xi') = 

81 {\left( 

82 \inf_{\Phi \in \Pi(\Xi, \Xi')} 

83 \int_{\mathbb{R}^m \times \mathbb{R}^m} 

84 \lVert \phi - \phi' \rVert^p \; 

85 \Phi( 

86 \mathop{}\!\mathrm{d} \phi 

87 \times 

88 \mathop{}\!\mathrm{d} \phi') 

89 \right)}^{\frac{1}{p}}, 

90 

91 where 

92 

93 1. :math:`\mathcal{M}` is the set of all Borel probability measures on 

94 :math:`\mathbb{R}^n` for some :math:`n \in \mathbb{Z}_{\geq 1}`, 

95 2. :math:`\Pi(\Xi, \Xi')` denotes the set of all couplings of :math:`\Xi` 

96 and :math:`\Xi'`, 

97 3. :math:`\xi_{(i)} \in \mathbb{R}^n` for all :math:`i \in [N]` are the 

98 :math:`N \in \mathbb{Z}_{\geq 1}` *samples* comprising the support of 

99 :math:`\Xi_\text{N}`, 

100 4. :math:`w_i \in \mathbb{R}_{\geq 0}` are *weights* denoting the nominal 

101 probabilitiy mass at :math:`\xi_{(i)}` for all :math:`i \in [N]`, 

102 5. :math:`\delta_{\xi_{(i)}}` denotes the Dirac delta function with unit 

103 mass at :math:`\xi_{(i)}` for all :math:`i \in [N]` and where 

104 6. :math:`\epsilon \in \mathbb{R}_{\geq 0}` controls the radius of the 

105 ambiguity set. 

106 

107 :Supported functions: 

108 

109 For :math:`p = 1`: 

110 

111 1. A convex piecewise linear function :math:`f(x, \xi) = max_{i=1}^k a_i(x, 

112 \xi)` where :math:`a` is biaffine in :math:`x` and :math:`\xi` for all 

113 :math:`i \in [k]`. This can be written as ``picos.max([a_1, ..., a_k])`` 

114 in Python. 

115 2. A concave piecewise linear function :math:`f(x, \xi) = min_{i=1}^k a_i(x, 

116 \xi)` where :math:`a` is biaffine in :math:`x` and :math:`\xi` for all 

117 :math:`i \in [k]`. This can be written as ``picos.min([a_1, ..., a_k])`` 

118 in Python. 

119 

120 For :math:`p = 2`: 

121 

122 1. A squared norm :math:`f(x, \xi) = \lVert A(x, \xi) \rVert_F^2` where 

123 :math:`A` is biaffine in :math:`x` and :math:`\xi`. This can be written 

124 as ``abs(A)**2`` in Python. 

125 """ 

126 

127 def __init__(self, parameter_name, p, eps, samples, weights=1): 

128 r"""Create a :class:`WassersteinAmbiguitySet`. 

129 

130 :param str parameter_name: 

131 Name of the random parameter :math:`\xi`. 

132 

133 :param float p: 

134 The Wasserstein type/order parameter :math:`p`. 

135 

136 :param float eps: 

137 The Wasserstein ball radius :math:`\epsilon`. 

138 

139 :param samples: 

140 The support of the discrete distribution :math:`\Xi_\text{D}` given 

141 as the *samples* :math:`\xi_{(i)}`. The original shape of the 

142 samples determines the shape of :math:`\xi`. 

143 :type samples: 

144 aynthing recognized by :class:`~.samples.Samples` 

145 

146 :param weights: 

147 A vector denoting the nonnegative weight (e.g. frequency or 

148 probability) of each sample. Its length must match the number of 

149 samples provided. The argument will be normalized such that its 

150 entries sum to one. Entries of zero will be dropped alongside their 

151 associated sample. The default value of ``1`` denotes the empirical 

152 distribution on the samples. 

153 

154 .. warning:: 

155 

156 Duplicate samples are not detected and can impact performance. If 

157 duplicate samples are likely, make sure to detect them and encode 

158 their frequency in the weight vector. 

159 """ 

160 # Load p. 

161 self._p = float(p) 

162 

163 if self._p < 1: 

164 raise ValueError("The Wasserstein parameter p must be >= 1.") 

165 

166 supported_p = (1, 2) 

167 if self._p not in supported_p: 

168 raise NotImplementedError("Currently, Wasserstein DRO is only " 

169 "supported for p in {}.".format(set(supported_p))) 

170 

171 # Load epsilon. 

172 self._eps = float(eps) 

173 

174 if self._eps < 0: 

175 raise ValueError("The Wasserstein ball radius must be nonnegative.") 

176 

177 # Load the samples. 

178 self._samples = Samples(samples) 

179 

180 # Load the normalized weights. 

181 w = AffineExpression.from_constant(weights, (len(self._samples), 1)) 

182 w_np = numpy.ravel(cvx2np(w.value_as_matrix)) 

183 

184 if any(w_np < 0): 

185 raise ValueError( 

186 "The weight vector must be nonnegative everywhere.") 

187 

188 if any(w_np == 0): 

189 if all(w_np == 0): 

190 raise ValueError("The weight vector must be nonzero.") 

191 

192 nonzero = numpy.where(w_np != 0)[0].tolist() 

193 w = w[nonzero] 

194 self._samples = self._samples.select(nonzero) 

195 

196 self._weights = (w / (w | 1)).renamed("w") 

197 

198 assert len(self._samples) == len(self._weights) 

199 

200 # Create the perturbation parameter. 

201 self._parameter = Perturbation( 

202 self, parameter_name, self._samples.original_shape) 

203 

204 @property 

205 def p(self): 

206 """The Wasserstein order :math:`p`.""" 

207 return self._p 

208 

209 @property 

210 def eps(self): 

211 r"""The Wasserstein ball radius :math:`\epsilon`.""" 

212 return self._eps 

213 

214 @property 

215 def samples(self): 

216 """The registered samples as a :class:`~.samples.Samples` object.""" 

217 return self._samples 

218 

219 @property 

220 def weights(self): 

221 """The sample weights a constant PICOS vector.""" 

222 return self._weights 

223 

224 Subtype = namedtuple("Subtype", ("sample_dim", "sample_num", "p")) 

225 

226 def _subtype(self): 

227 return self.Subtype(self._samples.dim, self._samples.num, self._p) 

228 

229 def __str__(self): 

230 return "WAS(p={}, eps={}, N={})".format( 

231 self._p, self._eps, self._samples.num) 

232 

233 @classmethod 

234 def _get_type_string_base(cls): 

235 return "Wasserstein Ambiguity Set" 

236 

237 def __repr__(self): 

238 return glyphs.repr2("{} {}".format(glyphs.shape(self._parameter.shape), 

239 self._get_type_string_base()), self.__str__()) 

240 

241 @property 

242 def distributional(self): 

243 """Implement for :class:`~.perturbation.PerturbationUniverse`.""" 

244 return True 

245 

246 @property 

247 def parameter(self): 

248 r"""The random perturbation parameter :math:`\xi`.""" 

249 return self._parameter 

250 

251 

252# -------------------------------------- 

253__all__ = api_end(_API_START, globals())