Source code for geomstats.geometry.symmetric_matrices

"""The vector space of symmetric matrices.

Lead author: Yann Thanwerdas.
"""

import geomstats.backend as gs
from geomstats.geometry.base import VectorSpace
from geomstats.geometry.matrices import Matrices, MatricesMetric


[docs] class SymmetricMatrices(VectorSpace): """Class for the vector space of symmetric matrices of size n. Parameters ---------- n : int Integer representing the shapes of the matrices: n x n. """ def __init__(self, n, equip=True): super().__init__(dim=int(n * (n + 1) / 2), shape=(n, n), equip=equip) self.n = n
[docs] @staticmethod def default_metric(): """Metric to equip the space with if equip is True.""" return MatricesMetric
def _create_basis(self): """Compute the basis of the vector space of symmetric matrices.""" indices, values = [], [] k = -1 for row in range(self.n): for col in range(row, self.n): k += 1 if row == col: indices.append((k, row, row)) values.append(1.0) else: indices.extend([(k, row, col), (k, col, row)]) values.extend([1.0, 1.0]) return gs.array_from_sparse(indices, values, (k + 1, self.n, self.n))
[docs] def belongs(self, point, atol=gs.atol): """Evaluate if a matrix is symmetric. Parameters ---------- point : array-like, shape=[.., n, n] Point to test. atol : float Tolerance to evaluate equality with the transpose. Returns ------- belongs : array-like, shape=[...,] Boolean evaluating if point belongs to the space. """ belongs = super().belongs(point) if gs.any(belongs): is_symmetric = Matrices.is_symmetric(point, atol) return gs.logical_and(belongs, is_symmetric) return belongs
[docs] def projection(self, point): """Make a matrix symmetric, by averaging with its transpose. Parameters ---------- point : array-like, shape=[..., n, n] Matrix. Returns ------- sym : array-like, shape=[..., n, n] Symmetric matrix. """ return Matrices.to_symmetric(point)
[docs] def random_point(self, n_samples=1, bound=1.0): """Sample a symmetric matrix. Samples from a uniform distribution in a box and then converts to symmetric. Parameters ---------- n_samples : int Number of samples. Optional, default: 1. bound : float Side of hypercube support of the uniform distribution. Optional, default: 1.0 Returns ------- point : array-like, shape=[..., n, n] Sample. """ sample = super().random_point(n_samples, bound) return Matrices.to_symmetric(sample)
[docs] @staticmethod def to_vector(point): """Convert a symmetric matrix into a vector. Parameters ---------- mat : array-like, shape=[..., n, n] Matrix. Returns ------- vec : array-like, shape=[..., n(n+1)/2] Vector. """ return gs.triu_to_vec(point)
[docs] @staticmethod def from_vector(vec): """Convert a vector into a symmetric matrix. Parameters ---------- vec : array-like, shape=[..., n(n+1)/2] Vector. Returns ------- mat : array-like, shape=[..., n, n] Symmetric matrix. """ vec_dim = vec.shape[-1] mat_dim = (gs.sqrt(8.0 * vec_dim + 1) - 1) / 2 if mat_dim != int(mat_dim): raise ValueError( "Invalid input dimension, it must be of the form" "(n_samples, n * (n + 1) / 2)" ) mat_dim = int(mat_dim) shape = (mat_dim, mat_dim) mask = 2 * gs.ones(shape) - gs.eye(mat_dim) indices = list(zip(*gs.triu_indices(mat_dim))) if gs.ndim(vec) == 1: upper_triangular = gs.array_from_sparse(indices, vec, shape) else: upper_triangular = gs.stack( [gs.array_from_sparse(indices, data, shape) for data in vec] ) mat = Matrices.to_symmetric(upper_triangular) * mask return mat