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# ------------------------------------------------------------------------------
20r"""Implement special matrix vectorization formats.
22These formats are used to efficiently store structured mutable types such as
23symmetric matrix variables in the form of real vectors.
24"""
26import random
27from abc import ABC, abstractmethod
28from copy import copy
30import cvxopt
32from .. import glyphs, settings
33from ..apidoc import api_end, api_start
34from .data import cvxopt_equals, cvxopt_hcat, cvxopt_vcat, load_shape
36_API_START = api_start(globals())
37# -------------------------------
40#: Number of instances to cache per vectorization format.
41CACHE_SIZE = 100
43#: Number of cached instances to drop at random when the cache is full.
44CACHE_BULK_REMOVE = CACHE_SIZE // 4
47class BaseVectorization(ABC, object):
48 """Abstract base class for special matrix vectorization formats.
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.
57 .. warning::
59 Due to how caching is implemented, derived classes may not inherit from
60 each other but only from :class:`BaseVectorization` directly!
61 """
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())
67 if not hasattr(cls, "_cache"):
68 cls._cache = {}
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)
76 cls._cache[shape] = object.__new__(cls)
78 return cls._cache[shape]
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()
86 def __len__(self):
87 return self._shape[0] * self._shape[1]
89 @property
90 def shape(self):
91 """The shape of matrices being vectorized."""
92 return self._shape
94 @property
95 def dim(self):
96 """The length of the vectorization.
98 This corresponds to the dimension of a matrix mutable being vectorized.
99 """
100 return self._special2full.size[1]
102 @classmethod
103 @abstractmethod
104 def _square_input(cls):
105 """Whether input matrices must be square."""
106 pass
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.")
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)))
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.")
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.")
130 if vector.typecode == "z":
131 raise TypeError("Cannot devectorize a complex vector: All "
132 "vectorizations are expected to be real.")
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)))
139 @abstractmethod
140 def _make_special_to_full(self):
141 """Return a mapping from the special to the full vectorization."""
142 pass
144 @abstractmethod
145 def _make_full_to_special(self):
146 """Return a mapping from the full to the special vectorization.
148 Returns :obj:`None` if the input matrix is complex. Then,
149 :meth:`vectorize` needs to be overridden.
150 """
151 pass
153 @property
154 def identity(self):
155 """A linear mapping from the special to the full vectorization.
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
165 def vectorize(self, matrix):
166 """Given a matrix, return its special vectorization.
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[:]
175 def devectorize(self, vector):
176 """Given a special vectorization, return the corresponding matrix.
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
187class FullVectorization(BaseVectorization):
188 """A basic column-major matrix vectorization."""
190 @classmethod
191 def _square_input(cls):
192 return False
194 def _validate_matrix(self, matrix):
195 BaseVectorization._validate_matrix(self, matrix)
196 self._ensure_real(matrix)
198 def _make_full_to_special(self):
199 n = len(self)
200 return cvxopt.spmatrix([1]*n, range(n), range(n), tc="d")
202 def vectorize(self, matrix):
203 """Override :meth:`BaseVectorization.vectorize` for speed reasons."""
204 self._validate_matrix(matrix)
206 if matrix.typecode == "d":
207 return matrix[:]
208 else:
209 assert matrix.typecode == "i"
210 return self._full2special*matrix[:]
212 def _make_special_to_full(self):
213 # Not actually used, see devectorize.
214 return self._make_full_to_special()
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
224class ComplexVectorization(BaseVectorization):
225 """An isometric vectorization that stacks real and imaginary parts."""
227 @classmethod
228 def _square_input(cls):
229 return False
231 def _validate_matrix(self, matrix):
232 BaseVectorization._validate_matrix(self, matrix)
234 def _make_full_to_special(self):
235 # No such matrix exists as the full vectorization is complex.
236 return None
238 def vectorize(self, matrix):
239 """Override :meth:`BaseVectorization.vectorize`.
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()])
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")
256class SymmetricVectorization(BaseVectorization):
257 """An isometric symmetric matrix vectorization.
259 See [svec]_ for the precise vectorization used.
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 """
267 @classmethod
268 def _square_input(cls):
269 return True
271 def _validate_matrix(self, matrix):
272 BaseVectorization._validate_matrix(self, matrix)
273 self._ensure_real(matrix)
275 if not cvxopt_equals(matrix, matrix.T,
276 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE):
277 raise ValueError("The given matrix is not numerically symmetric.")
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
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])
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)
304 return cvxopt.spmatrix(V, I, J, (m, n**2), tc="d")
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")
318class SkewSymmetricVectorization(BaseVectorization):
319 """An isometric skew-symmetric matrix vectorization."""
321 @classmethod
322 def _square_input(cls):
323 return True
325 def _validate_matrix(self, matrix):
326 BaseVectorization._validate_matrix(self, matrix)
327 self._ensure_real(matrix)
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.")
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")
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")
365class HermitianVectorization(BaseVectorization):
366 r"""An isometric hermitian matrix vectorization.
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 """
373 def __init__(self, shape):
374 """Initialize a vectorization format for hermitian matrices.
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)
383 @classmethod
384 def _square_input(cls):
385 return True
387 def _validate_matrix(self, matrix):
388 BaseVectorization._validate_matrix(self, matrix)
390 if not cvxopt_equals(matrix, matrix.H,
391 relTol=settings.RELATIVE_HERMITIANNESS_TOLERANCE):
392 raise ValueError("The given matrix is not numerically hermitian.")
394 def _make_full_to_special(self):
395 # No such matrix exists as the full vectorization is complex.
396 return None
398 def vectorize(self, matrix):
399 """Override :meth:`BaseVectorization.vectorize`.
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()])
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
421class LowerTriangularVectorization(BaseVectorization):
422 """An isometric lower triangular matrix vectorization."""
424 @classmethod
425 def _square_input(cls):
426 return True
428 def _validate_matrix(self, matrix):
429 BaseVectorization._validate_matrix(self, matrix)
430 self._ensure_real(matrix)
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.")
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")
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")
464class UpperTriangularVectorization(BaseVectorization):
465 """An isometric upper triangular matrix vectorization."""
467 @classmethod
468 def _square_input(cls):
469 return True
471 def _validate_matrix(self, matrix):
472 BaseVectorization._validate_matrix(self, matrix)
473 self._ensure_real(matrix)
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.")
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")
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")
507# --------------------------------------
508__all__ = api_end(_API_START, globals())