Coverage for picos/expressions/vectorizations.py: 78.15%

270 statements  

« prev     ^ index     » next       coverage.py v6.5.0, created at 2023-03-26 07:46 +0000

1# ------------------------------------------------------------------------------ 

2# Copyright (C) 2019 Maximilian Stahlberg 

3# Based on the original picos.expressions module by Guillaume Sagnol. 

4# 

5# This file is part of PICOS. 

6# 

7# PICOS is free software: you can redistribute it and/or modify it under the 

8# terms of the GNU General Public License as published by the Free Software 

9# Foundation, either version 3 of the License, or (at your option) any later 

10# version. 

11# 

12# PICOS is distributed in the hope that it will be useful, but WITHOUT ANY 

13# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

14# A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

15# 

16# You should have received a copy of the GNU General Public License along with 

17# this program. If not, see <http://www.gnu.org/licenses/>. 

18# ------------------------------------------------------------------------------ 

19 

20r"""Implement special matrix vectorization formats. 

21 

22These formats are used to efficiently store structured mutable types such as 

23symmetric matrix variables in the form of real vectors. 

24""" 

25 

26import random 

27from abc import ABC, abstractmethod 

28from copy import copy 

29 

30import cvxopt 

31 

32from .. import glyphs, settings 

33from ..apidoc import api_end, api_start 

34from .data import cvxopt_equals, cvxopt_hcat, cvxopt_vcat, load_shape 

35 

36_API_START = api_start(globals()) 

37# ------------------------------- 

38 

39 

40#: Number of instances to cache per vectorization format. 

41CACHE_SIZE = 100 

42 

43#: Number of cached instances to drop at random when the cache is full. 

44CACHE_BULK_REMOVE = CACHE_SIZE // 4 

45 

46 

47class BaseVectorization(ABC, object): 

48 """Abstract base class for special matrix vectorization formats. 

49 

50 Subclass instances are cached: If multiple instances of the same 

51 vectorization format and concerning matrices of the same shape are requested 

52 successively, then the instance created to serve the first request is 

53 retrieved from a cache on successive requests. The module attributes 

54 :data:`CACHE_SIZE` and :data:`CACHE_BULK_REMOVE` control the size of the 

55 cache for each vectorization format. 

56 

57 .. warning:: 

58 

59 Due to how caching is implemented, derived classes may not inherit from 

60 each other but only from :class:`BaseVectorization` directly! 

61 """ 

62 

63 def __new__(cls, shape): 

64 """Lookup or create a vectorization format for a fixed matrix shape.""" 

65 shape = load_shape(shape, squareMatrix=cls._square_input()) 

66 

67 if not hasattr(cls, "_cache"): 

68 cls._cache = {} 

69 

70 if shape not in cls._cache: 

71 if len(cls._cache) >= CACHE_SIZE: 

72 for remove in random.sample( 

73 list(cls._cache.keys()), CACHE_BULK_REMOVE): 

74 cls._cache.pop(remove) 

75 

76 cls._cache[shape] = object.__new__(cls) 

77 

78 return cls._cache[shape] 

79 

80 def __init__(self, shape): 

81 """Initialize a vectorization format for a fixed matrix shape.""" 

82 self._shape = load_shape(shape, squareMatrix=self._square_input()) 

83 self._special2full = self._make_special_to_full() 

84 self._full2special = self._make_full_to_special() 

85 

86 def __len__(self): 

87 return self._shape[0] * self._shape[1] 

88 

89 @property 

90 def shape(self): 

91 """The shape of matrices being vectorized.""" 

92 return self._shape 

93 

94 @property 

95 def dim(self): 

96 """The length of the vectorization. 

97 

98 This corresponds to the dimension of a matrix mutable being vectorized. 

99 """ 

100 return self._special2full.size[1] 

101 

102 @classmethod 

103 @abstractmethod 

104 def _square_input(cls): 

105 """Whether input matrices must be square.""" 

106 pass 

107 

108 @abstractmethod 

109 def _validate_matrix(self, matrix): 

110 """Raise an exception if the given matrix cannot be vectorized.""" 

111 if not isinstance(matrix, (cvxopt.matrix, cvxopt.spmatrix)): 

112 raise TypeError("May only vectorize CVXOPT matrix types.") 

113 

114 if matrix.size != self._shape: 

115 raise TypeError("Cannot vectorize a matrix of shape {} according " 

116 "to a vectorization recipe for {} matrices." 

117 .format(glyphs.shape(matrix.size), glyphs.shape(self._shape))) 

118 

119 def _ensure_real(self, matrix): 

120 """Raise a :exc:`TypeError` if the given matrix is not real.""" 

121 if matrix.typecode == "z": 

122 raise TypeError( 

123 "The vectorization format does not support complex input.") 

124 

125 def _validate_vector(self, vector): 

126 """Raise an exception if the given vector cannot be devectorized.""" 

127 if not isinstance(vector, (cvxopt.matrix, cvxopt.spmatrix)): 

128 raise TypeError("May only devectorize CVXOPT matrix types.") 

129 

130 if vector.typecode == "z": 

131 raise TypeError("Cannot devectorize a complex vector: All " 

132 "vectorizations are expected to be real.") 

133 

134 if vector.size != (self.dim, 1): 

135 raise TypeError( 

136 "Invalid shape of vectorized data: Expected {} but got {}." 

137 .format(glyphs.shape((self.dim, 1)), glyphs.shape(vector.size))) 

138 

139 @abstractmethod 

140 def _make_special_to_full(self): 

141 """Return a mapping from the special to the full vectorization.""" 

142 pass 

143 

144 @abstractmethod 

145 def _make_full_to_special(self): 

146 """Return a mapping from the full to the special vectorization. 

147 

148 Returns :obj:`None` if the input matrix is complex. Then, 

149 :meth:`vectorize` needs to be overridden. 

150 """ 

151 pass 

152 

153 @property 

154 def identity(self): 

155 """A linear mapping from the special to the full vectorization. 

156 

157 The term *identity* comes from the fact that these matrices are used 

158 as the coefficients that map the internal (vectorized) representation of 

159 a :class:`~.mutable.Mutable` object to the 

160 :class:`~.exp_biaffine.BiaffineExpression` instance that represents the 

161 mutable in algebraic operations. 

162 """ 

163 return self._special2full 

164 

165 def vectorize(self, matrix): 

166 """Given a matrix, return its special vectorization. 

167 

168 :raises TypeError: If the input isn't a CVXOPT matrix or does not have 

169 the expected numeric type or shape. 

170 :raises ValueError: If the matrix does not have the expected structure. 

171 """ 

172 self._validate_matrix(matrix) 

173 return self._full2special*matrix[:] 

174 

175 def devectorize(self, vector): 

176 """Given a special vectorization, return the corresponding matrix. 

177 

178 :raises TypeError: If the input isn't a CVXOPT column vector or does not 

179 have the expected numeric type or length. 

180 """ 

181 self._validate_vector(vector) 

182 M = self._special2full*vector 

183 M.size = self._shape 

184 return M 

185 

186 

187class FullVectorization(BaseVectorization): 

188 """A basic column-major matrix vectorization.""" 

189 

190 @classmethod 

191 def _square_input(cls): 

192 return False 

193 

194 def _validate_matrix(self, matrix): 

195 BaseVectorization._validate_matrix(self, matrix) 

196 self._ensure_real(matrix) 

197 

198 def _make_full_to_special(self): 

199 n = len(self) 

200 return cvxopt.spmatrix([1]*n, range(n), range(n), tc="d") 

201 

202 def vectorize(self, matrix): 

203 """Override :meth:`BaseVectorization.vectorize` for speed reasons.""" 

204 self._validate_matrix(matrix) 

205 

206 if matrix.typecode == "d": 

207 return matrix[:] 

208 else: 

209 assert matrix.typecode == "i" 

210 return self._full2special*matrix[:] 

211 

212 def _make_special_to_full(self): 

213 # Not actually used, see devectorize. 

214 return self._make_full_to_special() 

215 

216 def devectorize(self, vector): 

217 """Override :meth:`BaseVectorization.devectorize` for speed reasons.""" 

218 self._validate_vector(vector) 

219 M = copy(vector) 

220 M.size = self._shape 

221 return M 

222 

223 

224class ComplexVectorization(BaseVectorization): 

225 """An isometric vectorization that stacks real and imaginary parts.""" 

226 

227 @classmethod 

228 def _square_input(cls): 

229 return False 

230 

231 def _validate_matrix(self, matrix): 

232 BaseVectorization._validate_matrix(self, matrix) 

233 

234 def _make_full_to_special(self): 

235 # No such matrix exists as the full vectorization is complex. 

236 return None 

237 

238 def vectorize(self, matrix): 

239 """Override :meth:`BaseVectorization.vectorize`. 

240 

241 This is necessary because extracting the real and the imaginary part 

242 cannot be done with a linear transformation matrix as is done for other 

243 vectorization formats. 

244 """ 

245 self._validate_matrix(matrix) 

246 fullVectorization = matrix[:] 

247 return cvxopt_vcat([fullVectorization.real(), fullVectorization.imag()]) 

248 

249 def _make_special_to_full(self): 

250 cRank = len(self) # Rank on the complex field. 

251 rRank = 2*cRank # Rank on the real field. 

252 return cvxopt.spmatrix([1]*cRank + [1j]*cRank, 

253 list(range(cRank))*2, range(rRank), tc="z") 

254 

255 

256class SymmetricVectorization(BaseVectorization): 

257 """An isometric symmetric matrix vectorization. 

258 

259 See [svec]_ for the precise vectorization used. 

260 

261 .. [svec] Dattorro, J. (2018). Isomorphism of symmetric matrix subspace. 

262 In *Convex Optimization & Euclidean Distance Geometry (2nd ed.)* 

263 (pp. 47f.). California, Meboo Publishing USA. Retrieved from 

264 `<https://meboo.convexoptimization.com/Meboo.html>`_. 

265 """ 

266 

267 @classmethod 

268 def _square_input(cls): 

269 return True 

270 

271 def _validate_matrix(self, matrix): 

272 BaseVectorization._validate_matrix(self, matrix) 

273 self._ensure_real(matrix) 

274 

275 if not cvxopt_equals(matrix, matrix.T, 

276 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE): 

277 raise ValueError("The given matrix is not numerically symmetric.") 

278 

279 def _make_full_to_special(self): 

280 n = self._shape[0] 

281 m = n*(n + 1) // 2 

282 I, J, V = [], [], [] 

283 for i in range(m): 

284 c = int((1 + 8*i)**0.5 - 1) // 2 

285 r = i - c*(c + 1) // 2 

286 

287 # Version 1: Average elements in lower and upper triangular part. 

288 # This could work around noisy input data. 

289 if c == r: 

290 I.append(i) 

291 J.append(c*n + r) 

292 V.append(1) 

293 else: 

294 I.extend([i, i]) 

295 J.extend([c*n + r, r*n + c]) 

296 V.extend([2**0.5 / 2, 2**0.5 / 2]) 

297 

298 # Version 2: Ignore elements below the main diagonal. 

299 # This is the faster approach. 

300 # I.append(i) 

301 # J.append(c*n + r) 

302 # V.append(1 if c == r else 2**0.5) 

303 

304 return cvxopt.spmatrix(V, I, J, (m, n**2), tc="d") 

305 

306 def _make_special_to_full(self): 

307 n = self._shape[0] 

308 I, J, V = range(n**2), [], [] 

309 for i in I: 

310 c, r = divmod(i, n) 

311 if c < r: # Entries below the diagonal are infered. 

312 c, r = r, c 

313 J.append(c*(c + 1) // 2 + r) 

314 V.append(1 if c == r else 1 / 2**0.5) 

315 return cvxopt.spmatrix(V, I, J, (n**2, n*(n + 1) // 2), tc="d") 

316 

317 

318class SkewSymmetricVectorization(BaseVectorization): 

319 """An isometric skew-symmetric matrix vectorization.""" 

320 

321 @classmethod 

322 def _square_input(cls): 

323 return True 

324 

325 def _validate_matrix(self, matrix): 

326 BaseVectorization._validate_matrix(self, matrix) 

327 self._ensure_real(matrix) 

328 

329 # NOTE: Use hermitianess tolerance. 

330 if not cvxopt_equals(matrix, -matrix.T, 

331 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE): 

332 raise ValueError( 

333 "The given matrix is not numerically skew-symmetric.") 

334 

335 def _make_full_to_special(self): 

336 n = self._shape[0] 

337 m = n*(n - 1) // 2 

338 I, J, V = [], [], [] 

339 for i in range(m): 

340 c = int((1 + 8*i)**0.5 + 1) // 2 

341 r = i - c*(c - 1) // 2 

342 I.append(i) 

343 J.append(c*n + r) 

344 V.append(2**0.5) 

345 return cvxopt.spmatrix(V, I, J, (m, n**2), tc="d") 

346 

347 def _make_special_to_full(self): 

348 n = self._shape[0] 

349 I, J, V = [], [], [] 

350 for i in range(n**2): 

351 c, r = divmod(i, n) 

352 if c == r: # Entries on the diagonal are zero. 

353 continue 

354 elif c < r: # Entries below the diagonal are infered. 

355 I.append(i) 

356 J.append(r*(r - 1) // 2 + c) 

357 V.append(-1 / 2**0.5) 

358 else: 

359 I.append(i) 

360 J.append(c*(c - 1) // 2 + r) 

361 V.append(1 / 2**0.5) 

362 return cvxopt.spmatrix(V, I, J, (n**2, n*(n - 1) // 2), tc="d") 

363 

364 

365class HermitianVectorization(BaseVectorization): 

366 r"""An isometric hermitian matrix vectorization. 

367 

368 The vectorization is isometric with respect to the Hermitian inner product 

369 :math:`\langle A, B \rangle = \operatorname{tr}(B^H A)` on the matrices and 

370 the real dot product on their vectorizations. 

371 """ 

372 

373 def __init__(self, shape): 

374 """Initialize a vectorization format for hermitian matrices. 

375 

376 Uses :class:`SymmetricVectorization` (for the real part) and 

377 :class:`SkewSymmetricVectorization` (for the imaginary part) internally. 

378 """ 

379 self._sym = SymmetricVectorization(shape) 

380 self._skw = SkewSymmetricVectorization(shape) 

381 super(HermitianVectorization, self).__init__(shape) 

382 

383 @classmethod 

384 def _square_input(cls): 

385 return True 

386 

387 def _validate_matrix(self, matrix): 

388 BaseVectorization._validate_matrix(self, matrix) 

389 

390 if not cvxopt_equals(matrix, matrix.H, 

391 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE): 

392 raise ValueError("The given matrix is not numerically hermitian.") 

393 

394 def _make_full_to_special(self): 

395 # No such matrix exists as the full vectorization is complex. 

396 return None 

397 

398 def vectorize(self, matrix): 

399 """Override :meth:`BaseVectorization.vectorize`. 

400 

401 This is necessary because extracting the real and the imaginary part 

402 cannot be done with a linear transformation matrix as is done for other 

403 vectorization formats. 

404 """ 

405 self._validate_matrix(matrix) 

406 fullVectorization = matrix[:] 

407 return cvxopt_vcat([ 

408 self._sym._full2special*fullVectorization.real(), 

409 self._skw._full2special*fullVectorization.imag()]) 

410 

411 def _make_special_to_full(self): 

412 A = cvxopt_hcat([ 

413 self._sym._special2full, 

414 cvxopt.spmatrix([], [], [], self._skw._special2full.size)]) 

415 B = cvxopt_hcat([ 

416 cvxopt.spmatrix([], [], [], self._sym._special2full.size), 

417 self._skw._special2full]) 

418 return A + 1j*B 

419 

420 

421class LowerTriangularVectorization(BaseVectorization): 

422 """An isometric lower triangular matrix vectorization.""" 

423 

424 @classmethod 

425 def _square_input(cls): 

426 return True 

427 

428 def _validate_matrix(self, matrix): 

429 BaseVectorization._validate_matrix(self, matrix) 

430 self._ensure_real(matrix) 

431 

432 n = self._shape[0] 

433 for i in range(n): 

434 for j in range(i + 1, n): 

435 if matrix[i, j]: 

436 raise ValueError( 

437 "The given matrix is not lower triangular.") 

438 

439 def _make_full_to_special(self): 

440 n = self._shape[0] 

441 m = n*(n + 1) // 2 

442 I, J, V = [], [], [] 

443 for i in range(m): 

444 c = int((1 + 8*i)**0.5 - 1) // 2 

445 r = i - c*(c + 1) // 2 

446 I.append(i) 

447 J.append(r*n + c) 

448 V.append(1) 

449 return cvxopt.spmatrix(V, I, J, (m, n**2), tc="d") 

450 

451 def _make_special_to_full(self): 

452 n = self._shape[0] 

453 m = n*(n + 1) // 2 

454 I, J, V = [], [], [] 

455 for j in range(m): 

456 c = int((1 + 8*j)**0.5 - 1) // 2 

457 r = j - c*(c + 1) // 2 

458 I.append(r*n + c) 

459 J.append(j) 

460 V.append(1) 

461 return cvxopt.spmatrix(V, I, J, (n**2, m), tc="d") 

462 

463 

464class UpperTriangularVectorization(BaseVectorization): 

465 """An isometric upper triangular matrix vectorization.""" 

466 

467 @classmethod 

468 def _square_input(cls): 

469 return True 

470 

471 def _validate_matrix(self, matrix): 

472 BaseVectorization._validate_matrix(self, matrix) 

473 self._ensure_real(matrix) 

474 

475 n = self._shape[0] 

476 for i in range(n): 

477 for j in range(i): 

478 if matrix[i, j]: 

479 raise ValueError( 

480 "The given matrix is not upper triangular.") 

481 

482 def _make_full_to_special(self): 

483 n = self._shape[0] 

484 m = n*(n + 1) // 2 

485 I, J, V = [], [], [] 

486 for i in range(m): 

487 c = int((1 + 8*i)**0.5 - 1) // 2 

488 r = i - c*(c + 1) // 2 

489 I.append(i) 

490 J.append(c*n + r) 

491 V.append(1) 

492 return cvxopt.spmatrix(V, I, J, (m, n**2), tc="d") 

493 

494 def _make_special_to_full(self): 

495 n = self._shape[0] 

496 m = n*(n + 1) // 2 

497 I, J, V = [], [], [] 

498 for j in range(m): 

499 c = int((1 + 8*j)**0.5 - 1) // 2 

500 r = j - c*(c + 1) // 2 

501 I.append(c*n + r) 

502 J.append(j) 

503 V.append(1) 

504 return cvxopt.spmatrix(V, I, J, (n**2, m), tc="d") 

505 

506 

507# -------------------------------------- 

508__all__ = api_end(_API_START, globals())