# Source code for pymanopt.manifolds.oblique

```import numpy as np
import numpy.linalg as la
import numpy.random as rnd

from pymanopt.manifolds.manifold import EuclideanEmbeddedSubmanifold

[docs]class Oblique(EuclideanEmbeddedSubmanifold):
"""
Manifold of matrices w/ unit-norm columns.

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

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

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

[docs]    def inner(self, X, U, V):
return float(np.tensordot(U, V))

[docs]    def norm(self, X, U):
return la.norm(U)

[docs]    def dist(self, X, Y):
XY = (X * Y).sum(0)
XY[XY > 1] = 1
U = np.arccos(XY)
return la.norm(U)

[docs]    def proj(self, X, H):
return H - X * ((X * H).sum(0)[np.newaxis, :])

[docs]    def ehess2rhess(self, X, egrad, ehess, U):
PXehess = self.proj(X, ehess)
# TODO(nkoep): Move the second summand to the 'weingarten' method
return PXehess - U * ((X * egrad).sum(0)[np.newaxis, :])

[docs]    def exp(self, X, U):
norm_U = np.sqrt((U ** 2).sum(0))[np.newaxis, :]

Y = X * np.cos(norm_U) + U * (np.sin(norm_U) / norm_U)

# For those columns where the step is too small, use a retraction.
exclude = np.nonzero(norm_U <= 4.5e-8)[-1]
Y[:, exclude] = self._normalize_columns(X[:, exclude] + U[:, exclude])

return Y

[docs]    def retr(self, X, U):
return self._normalize_columns(X + U)

[docs]    def log(self, X, Y):
V = self.proj(X, Y - X)
dists = np.arccos((X * Y).sum(0))
norms = np.sqrt((V ** 2).sum(0)).real
factors = dists / norms
# For very close points, dists is almost equal to norms, but because
# they are both almost zero, the division above can return NaN's. To
# avoid that, we force those ratios to 1.
factors[dists <= 1e-6] = 1

return V * factors

[docs]    def rand(self):
return self._normalize_columns(rnd.randn(self._m, self._n))

[docs]    def randvec(self, X):
H = rnd.randn(*X.shape)
P = self.proj(X, H)
return P / self.norm(X, P)

[docs]    def transp(self, X, Y, U):
return self.proj(Y, U)

[docs]    def pairmean(self, X, Y):
return self._normalize_columns(X + Y)

[docs]    def zerovec(self, X):
return np.zeros((self._m, self._n))

def _normalize_columns(self, X):
"""Return an l2-column-normalized copy of the matrix X."""
return X / la.norm(X, axis=0)[np.newaxis, :]
```