Source code for geomstats.geometry.hermitian_matrices

"""The vector space of Hermitian matrices.

Lead author: Yann Cabanes.
import logging

import geomstats.backend as gs
from geomstats import algebra_utils as utils
from geomstats.geometry.base import ComplexVectorSpace
from geomstats.geometry.complex_matrices import ComplexMatrices, ComplexMatricesMetric
from geomstats.geometry.matrices import Matrices

[docs] def expmh(mat): """Compute the matrix exponential for a Hermitian matrix. Parameters ---------- mat : array_like, shape=[..., n, n] Symmetric matrix. Returns ------- exponential : array_like, shape=[..., n, n] Exponential of mat. """ n = mat.shape[-1] dim_3_mat = gs.reshape(mat, [-1, n, n]) expm = apply_func_to_eigvalsh(dim_3_mat, gs.exp) return gs.reshape(expm, mat.shape)
[docs] def powermh(mat, power): """Compute the matrix power for a Hermitian matrix. Parameters ---------- mat : array_like, shape=[..., n, n] Symmetric matrix with non-negative eigenvalues. power : float, list Power at which mat will be raised. If a list of powers is passed, a list of results will be returned. Returns ------- powerm : array_like or list of arrays, shape=[..., n, n] Matrix power of mat. """ if isinstance(power, list): power_ = [lambda ev, p=p: gs.power(ev, p) for p in power] else: def power_(ev): return gs.power(ev, power) return apply_func_to_eigvalsh(mat, power_, check_positive=False)
[docs] def apply_func_to_eigvalsh(mat, function, check_positive=False): """Apply function to eigenvalues and reconstruct the matrix. Parameters ---------- mat : array_like, shape=[..., n, n] Hermitian matrix. function : callable, list of callables Function to apply to eigenvalues. If a list of functions is passed, a list of results will be returned. check_positive : bool Whether to check positivity of the eigenvalues. Optional. Default: False. Returns ------- mat : array_like, shape=[..., n, n] Hermitian matrix. """ eigvals, eigvecs = gs.linalg.eigh(mat) if check_positive and gs.any(gs.cast(eigvals, gs.get_default_dtype()) < 0.0): try: name = function.__name__ except AttributeError: name = function[0].__name__ logging.warning("Negative eigenvalue encountered in %s", name) return_list = True if not isinstance(function, list): function = [function] return_list = False reconstruction = [] if gs.is_complex(mat): transp_eigvecs = ComplexMatrices.transconjugate(eigvecs) else: transp_eigvecs = Matrices.transpose(eigvecs) for fun in function: eigvals_f = fun(eigvals) eigvals_f = utils.from_vector_to_diagonal_matrix(eigvals_f) reconstruction.append(Matrices.mul(eigvecs, eigvals_f, transp_eigvecs)) return reconstruction if return_list else reconstruction[0]
[docs] class HermitianMatrices(ComplexVectorSpace): """Class for the vector space of Hermitian 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=n * (n + 1) - n, shape=(n, n), equip=equip) self.n = n
[docs] @staticmethod def default_metric(): """Metric to equip the space with if equip is True.""" return ComplexMatricesMetric
def _create_basis(self): """Compute the basis of the vector space of symmetric matrices. Returns ------- basis : array-like, shape=[dim, n, n] """ diagonal = [] real_part = [] complex_part = [] for row in gs.arange(self.n): for col in gs.arange(row, self.n): if row == col: indices = [(row, row)] values = [1.0 + 0j] diagonal.append( gs.array_from_sparse(indices, values, (self.n,) * 2) ) else: indices = [(row, col), (col, row)] values = [1.0 + 0j, 1.0 + 0j] real_part.append( gs.array_from_sparse(indices, values, (self.n,) * 2) ) values = [1j, -1j] complex_part.append( gs.array_from_sparse(indices, values, (self.n,) * 2) ) return gs.vstack( [gs.stack(diagonal), gs.stack(real_part), gs.stack(complex_part)] )
[docs] def belongs(self, point, atol=gs.atol): """Evaluate if a matrix is Hermitian. 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_hermitian = ComplexMatrices.is_hermitian(point, atol) return gs.logical_and(belongs, is_hermitian) return belongs
[docs] @staticmethod def projection(point): """Make a matrix Hermitian, by averaging with its transconjugate. Parameters ---------- point : array-like, shape=[..., n, n] Matrix. Returns ------- herm : array-like, shape=[..., n, n] Symmetric matrix. """ return ComplexMatrices.to_hermitian(point)
[docs] def random_point(self, n_samples=1, bound=1.0): """Sample a Hermitian matrix. Points are generated by sampling complex matrices from a uniform distribution in a box and averaging with the transconjugate. 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. """ cdtype = gs.get_default_cdtype() size = self.shape if n_samples != 1: size = (n_samples,) + self.shape point = gs.cast( bound * (gs.random.rand(*size) - 0.5) * 2**0.5, dtype=cdtype, ) + 1j * gs.cast( bound * (gs.random.rand(*size) - 0.5) * 2**0.5, dtype=cdtype, ) return ComplexMatrices.to_hermitian(point)
[docs] @staticmethod def to_vector(point): """Convert a Hermitian matrix into a vector. Parameters ---------- point : array-like, shape=[..., n, n] Matrix. Returns ------- vec : array-like, shape=[..., n(n+1)/2] Vector. """ diag = Matrices.diagonal(point) up_triang = gs.triu_to_vec(point, k=1) real_part = gs.real(up_triang) complex_part = gs.imag(up_triang) vec = gs.hstack([diag, real_part, complex_part]) return vec
[docs] def from_vector(self, vec): """Convert a vector into a Hermitian matrix. Parameters ---------- vec : array-like, shape=[..., dim] Vector. Returns ------- mat : array-like, shape=[..., n, n] Hermitian matrix. """ return gs.einsum("...k,...kij->...ij", vec, self.basis)