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:`SquaredNorm`.""" 

20 

21import operator 

22from collections import namedtuple 

23 

24import cvxopt 

25 

26from .. import glyphs 

27from ..apidoc import api_end, api_start 

28from ..caching import cached_property, cached_unary_operator 

29from ..constraints import SquaredNormConstraint 

30from .data import convert_and_refine_arguments, convert_operands, cvxopt_hcat 

31from .exp_affine import AffineExpression, ComplexAffineExpression 

32from .exp_quadratic import QuadraticExpression 

33from .expression import Expression, refine_operands, validate_prediction 

34 

35_API_START = api_start(globals()) 

36# ------------------------------- 

37 

38 

39class SquaredNorm(QuadraticExpression): 

40 """A squared Euclidean or Frobenius norm. 

41 

42 This is a lightweight wrapper around 

43 :class:`~picos.expressions.QuadraticExpression` that can handle common 

44 constraint formulations more efficiently. 

45 """ 

46 

47 # -------------------------------------------------------------------------- 

48 # Initialization and factory methods. 

49 # -------------------------------------------------------------------------- 

50 

51 @convert_and_refine_arguments("x") 

52 def __init__(self, x): 

53 """Create a squared Euclidean or Frobenius norm. 

54 

55 :param x: 

56 The (complex) affine expression under the squared norm. 

57 :type affinePart: 

58 ~picos.expressions.ComplexAffineExpression 

59 """ 

60 # Validate x. 

61 if not isinstance(x, ComplexAffineExpression): 

62 raise TypeError("Can only form the squared norm of an affine " 

63 "expression, not of {}.".format(type(x).__name__)) 

64 

65 if len(x) == 1: 

66 typeStr = "Squared Scalar" 

67 symbStr = glyphs.squared(x.string) 

68 else: 

69 typeStr = "Squared Norm" 

70 symbStr = glyphs.squared(glyphs.norm(x.string)) 

71 

72 Expression.__init__(self, typeStr, symbStr) 

73 

74 # TODO: Add a nonzero-vectorization to BiaffineExpression that returns 

75 # a scalar zero for an all-zero expression and the vectorization 

76 # of the expression with zero rows removed otherwise. 

77 if x.is0: 

78 self._x = AffineExpression.zero() 

79 else: 

80 # Vectorize and stack real and imaginary parts. 

81 vec = x.vec if x.isreal else x.vec.real // x.vec.imag 

82 

83 # Remove zero rows from the vectorization. 

84 A = abs(cvxopt_hcat(vec._coefs.values())) 

85 a = cvxopt.sparse(sum(A[:, j] for j in range(A.size[1]))) 

86 nonzero = a.I 

87 nnz = len(nonzero) 

88 B = cvxopt.spmatrix([1.0]*nnz, range(nnz), nonzero, (nnz, len(x))) 

89 

90 self._x = B*vec 

91 

92 # -------------------------------------------------------------------------- 

93 # Allow inheriting from QuadraticExpression. 

94 # -------------------------------------------------------------------------- 

95 

96 @cached_property 

97 def _quadratic_form(self): 

98 """The squared norm as a pure quadratic expression. 

99 

100 If the expression under the norm is constant, then this is :obj:`None`. 

101 """ 

102 # HACK: Make a shallow copy of self._x so that the product does not 

103 # recognize that the operation below represents a squared norm. 

104 # This only works as long as the product checks for operand 

105 # equality with the "is" keyword. 

106 result = (self._x.renamed("HACK") | self._x) 

107 

108 if isinstance(result, AffineExpression): 

109 return None 

110 else: 

111 assert isinstance(result, QuadraticExpression) 

112 result._symbStr = self.string 

113 return result 

114 

115 @cached_property 

116 def _quads(self): 

117 if self._quadratic_form: 

118 return self._quadratic_form._quads 

119 else: 

120 return {} 

121 

122 @cached_property 

123 def _aff(self): 

124 if self._quadratic_form: 

125 return self._quadratic_form._aff 

126 else: 

127 refined = self.refined 

128 assert isinstance(refined, ComplexAffineExpression) 

129 return refined 

130 

131 @cached_property 

132 def _sf(self): 

133 if len(self._x) == 1: 

134 return (self._x, self._x) 

135 else: 

136 return None 

137 

138 # -------------------------------------------------------------------------- 

139 # Squared norm specific properties. 

140 # -------------------------------------------------------------------------- 

141 

142 @property 

143 def argdim(self): 

144 """Number of nonzero elements of the expression under the norm.""" 

145 return len(self._x) 

146 

147 # -------------------------------------------------------------------------- 

148 # Abstract method implementations and method overridings, except _predict. 

149 # -------------------------------------------------------------------------- 

150 

151 @cached_unary_operator 

152 def _get_refined(self): 

153 if self._x.constant: 

154 value = self._x.value_as_matrix 

155 return AffineExpression.from_constant( 

156 value.T*value, (1, 1), self._symbStr) 

157 else: 

158 return self 

159 

160 Subtype = namedtuple("Subtype", ("argdim", "quadratic_subtype")) 

161 

162 @cached_unary_operator 

163 def _get_subtype(self): 

164 return self.Subtype( 

165 len(self._x), 

166 self._quadratic_form.subtype if self._quadratic_form else None) 

167 

168 def _get_value(self): 

169 value = self._x._get_value() 

170 return value.T*value 

171 

172 @cached_unary_operator 

173 def _get_variables(self): 

174 return self._x.variables 

175 

176 def _is_convex(self): 

177 return True 

178 

179 def _is_concave(self): 

180 return self._x.constant 

181 

182 def _replace_variables(self, var_map): 

183 return self.__class__(self._x._replace_variables(var_map)) 

184 

185 # -------------------------------------------------------------------------- 

186 # Python special method implementations, except constraint-creating ones. 

187 # -------------------------------------------------------------------------- 

188 

189 def __len__(self): 

190 # Faster version that overrides Expression.__len__. 

191 return 1 

192 

193 @convert_operands(sameShape=True) 

194 @refine_operands() 

195 def __add__(self, other): 

196 """Denote addition from the right hand side.""" 

197 if isinstance(other, SquaredNorm): 

198 result = SquaredNorm(self._x // other._x) 

199 result._symbStr = glyphs.clever_add(self.string, other.string) 

200 return result 

201 else: 

202 return QuadraticExpression.__add__(self, other) 

203 

204 @convert_operands(scalarRHS=True) 

205 @refine_operands() 

206 def __mul__(self, other): 

207 """Denote scaling from the right hand side.""" 

208 if isinstance(other, AffineExpression): 

209 if not other.constant: 

210 raise TypeError("You may only multiply a squared norm with a " 

211 "constant term.") 

212 

213 value = other.value 

214 

215 if value < 0: 

216 return QuadraticExpression.__mul__(self, other) 

217 

218 result = SquaredNorm(self._x*value**0.5) 

219 result._symbStr = glyphs.clever_mul(self.string, other.string) 

220 return result 

221 else: 

222 return QuadraticExpression.__mul__(self, other) 

223 

224 @convert_operands(scalarRHS=True) 

225 @refine_operands() 

226 def __rmul__(self, other): 

227 """Denote scaling from the left hand side.""" 

228 if isinstance(other, AffineExpression): 

229 result = self.__mul__(other) 

230 # NOTE: __mul__ always creates a fresh expression. 

231 result._symbStr = glyphs.clever_mul(other.string, self.string) 

232 return result 

233 else: 

234 return QuadraticExpression.__rmul__(self, other) 

235 

236 @convert_operands(scalarRHS=True) 

237 @refine_operands() 

238 def __truediv__(self, other): 

239 """Denote division by a constant scalar.""" 

240 if isinstance(other, AffineExpression): 

241 if not other.constant: 

242 raise TypeError("You may only divide a squared norm by a " 

243 "constant term.") 

244 

245 result = self.__mul__(1.0 / other.value) 

246 # NOTE: __mul__ always creates a fresh expression. 

247 result._symbStr = glyphs.div(self.string, other.string) 

248 return result 

249 else: 

250 return QuadraticExpression.__truediv__(self, other) 

251 

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

253 # Method overridings for QuadraticExpression. 

254 # -------------------------------------------------------------------------- 

255 

256 @cached_property 

257 def fullroot(self): 

258 """Affine expression whose squared norm equals the expression. 

259 

260 Overrides :meth:`~.exp_quadratic.QuadraticExpression.fullroot` of 

261 :class:`~.exp_quadratic.QuadraticExpression`. 

262 """ 

263 return self._x.renamed(glyphs.Fn("fullroot({})")(self.string)) 

264 

265 @property 

266 def is_squared_norm(self): 

267 """Always :obj:`True` for squared norm instances. 

268 

269 Overrides :meth:`~.exp_quadratic.QuadraticExpression.is_squared_norm` 

270 of :class:`~.exp_quadratic.QuadraticExpression`. 

271 """ 

272 return True 

273 

274 @property 

275 def is0(self): 

276 """Whether the expression is zero. 

277 

278 Overrides :meth:`~.exp_quadratic.QuadraticExpression.is0` of 

279 :class:`~.exp_quadratic.QuadraticExpression`. 

280 """ 

281 return self._x.is0 

282 

283 # -------------------------------------------------------------------------- 

284 # Constraint-creating operators, and _predict. 

285 # -------------------------------------------------------------------------- 

286 

287 @classmethod 

288 def _predict(cls, subtype, relation, other): 

289 assert isinstance(subtype, cls.Subtype) 

290 

291 if relation == operator.__le__: 

292 if issubclass(other.clstype, AffineExpression) \ 

293 and other.subtype.dim == 1: 

294 return SquaredNormConstraint.make_type( 

295 subtype.argdim, other.subtype.constant) 

296 

297 return QuadraticExpression._predict( 

298 subtype.quadratic_subtype, relation, other) 

299 

300 @convert_operands(sameShape=True) 

301 @validate_prediction 

302 @refine_operands() 

303 def __le__(self, other): 

304 if isinstance(other, AffineExpression): 

305 return SquaredNormConstraint(self, other) 

306 

307 # NOTE: The following should handle the case where the upper bound has a 

308 # scalar factorization efficiently by virtue of 

309 # SquaredNorm.fullroot. See ConicQuadraticConstraint. 

310 return QuadraticExpression.__le__(self, other) 

311 

312 

313# -------------------------------------- 

314__all__ = api_end(_API_START, globals())