# Source code for pymanopt.manifolds.oblique

import numpy as np

from pymanopt.manifolds.manifold import RiemannianSubmanifold

[docs]class Oblique(RiemannianSubmanifold):
r"""Manifold of matrices with unit-norm columns.

The oblique manifold deals with matrices of size m x n such that each
column has unit Euclidean norm, i.e., is a point on the unit sphere in
:math:\R^m.
The metric is such that the oblique manifold is a Riemannian submanifold of
the space of m x n matrices with the usual trace inner product.

Args:
m: The number of rows of each matrix.
n: The number of columns of each matrix.
"""

def __init__(self, m: int, n: int):
self._m = m
self._n = n
name = f"Oblique manifold OB({m},{n})"
dimension = (m - 1) * n
super().__init__(name, dimension)

@property
def typical_dist(self):
return np.pi * np.sqrt(self._n)

[docs]    def inner_product(self, point, tangent_vector_a, tangent_vector_b):
return np.tensordot(
tangent_vector_a, tangent_vector_b, axes=tangent_vector_a.ndim
)

[docs]    def norm(self, point, tangent_vector):
return np.linalg.norm(tangent_vector)

[docs]    def dist(self, point_a, point_b):
XY = (point_a * point_b).sum(0)
XY[XY > 1] = 1
return np.linalg.norm(np.arccos(XY))

[docs]    def projection(self, point, vector):
return vector - point * ((point * vector).sum(0)[np.newaxis, :])

to_tangent_space = projection

[docs]    def euclidean_to_riemannian_hessian(
self, point, euclidean_gradient, euclidean_hessian, tangent_vector
):
PXehess = self.projection(point, euclidean_hessian)
return PXehess - tangent_vector * (
(point * euclidean_gradient).sum(0)[np.newaxis, :]
)

[docs]    def exp(self, point, tangent_vector):
norm = np.sqrt((tangent_vector**2).sum(0))[np.newaxis, :]
target_point = point * np.cos(norm) + tangent_vector * np.sinc(
norm / np.pi
)
return target_point

[docs]    def retraction(self, point, tangent_vector):
return self._normalize_columns(point + tangent_vector)

[docs]    def log(self, point_a, point_b):
vector = self.projection(point_a, point_b - point_a)
distances = np.arccos((point_a * point_b).sum(0))
norms = np.sqrt((vector**2).sum(0)).real
# Try to avoid zero-division when both distances and norms are almost
# zero.
epsilon = np.finfo(np.float64).eps
factors = (distances + epsilon) / (norms + epsilon)
return vector * factors

[docs]    def random_point(self):
return self._normalize_columns(
np.random.normal(size=(self._m, self._n))
)

[docs]    def random_tangent_vector(self, point):
vector = np.random.normal(size=point.shape)
tangent_vector = self.projection(point, vector)
return tangent_vector / self.norm(point, tangent_vector)

[docs]    def transport(self, point_a, point_b, tangent_vector_a):
return self.projection(point_b, tangent_vector_a)

[docs]    def pair_mean(self, point_a, point_b):
return self._normalize_columns(point_a + point_b)

[docs]    def zero_vector(self, point):
return np.zeros((self._m, self._n))

def _normalize_columns(self, array):
return array / np.linalg.norm(array, axis=0)[np.newaxis, :]