Source code for geomstats.geometry.full_rank_matrices

r"""Full rank Euclidean matrices :math:R_*^{m\times n}.

Lead author: Anna Calissano.
"""

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

[docs]
class FullRankMatrices(VectorSpaceOpenSet):
r"""Class for :math:R_*^{n\times k} matrices of dimension n x k and full rank.

Parameters
----------
n : int
Integer representing the shape of the matrices: n x k
k : int
Integer representing the shape of the matrices: n x k
"""

def __init__(self, n, k, equip=True):
super().__init__(dim=n * k, embedding_space=Matrices(n, k), equip=equip)
self.rank = min(n, k)
self.n = n
self.k = k

[docs]
@staticmethod
def default_metric():
"""Metric to equip the space with if equip is True."""
return MatricesMetric

[docs]
def belongs(self, point, atol=gs.atol):
r"""Check if the matrix belongs to :math:R_*^{n \times k}.

Parameters
----------
point : array-like, shape=[..., n, k]
Matrix to be checked.
atol : float
Unused.

Returns
-------
belongs : Boolean
Denoting if point is in :math:R_*^{m\times n}.
"""
has_right_size = self.embedding_space.belongs(point)
has_right_rank = gs.where(
gs.linalg.matrix_rank(point) == self.rank, True, False
)
return gs.logical_and(gs.array(has_right_size), has_right_rank)

[docs]
def projection(self, point):
r"""Project a matrix to the set of full rank matrices.

As the space of full rank matrices is dense in the space of matrices,
this is not a projection per se, but a regularization if the matrix input X
is not already full rank: :math:X + \epsilon [I_{rank}, 0] is returned
where :math:\epsilon=gs.atol

Parameters
----------
point : array-like, shape=[..., n, k]
Point in embedding manifold.

Returns
-------
projected : array-like, shape=[..., n, k]
Projected point.
"""
belongs = self.belongs(point)
regularization = gs.einsum(
"...,ij->...ij",
gs.where(~belongs, gs.atol, 0.0),
gs.eye(self.embedding_space.shape[0], self.embedding_space.shape[1]),
)
projected = point + regularization
return projected

[docs]
def random_point(self, n_samples=1, bound=1.0, n_iter=100):
r"""Sample in :math:R_*^{n\times k} from a normal distribution.

Parameters
----------
n_samples : int
Number of samples.
Optional, default: 1.
bound: float
This parameter is ignored.
n_iter : int
Maximum number of trials to sample a matrix with full rank.
Optional, default: 100.

Returns
-------
samples : array-like, shape=[..., n, k]
Point sampled on :math:R_*^{n\times k}.
"""
sample = []
n_accepted, iteration = 0, 0
while n_accepted < n_samples and iteration < n_iter:
raw_samples = gs.random.normal(
size=(n_samples - n_accepted, self.n, self.k)
)
ranks = gs.linalg.matrix_rank(raw_samples)
selected = ranks == self.rank
sample.append(raw_samples[selected])
n_accepted += gs.sum(selected)
iteration += 1
if n_samples == 1:
return sample[0][0]
return gs.concatenate(sample)