Source code for pymanopt.manifolds.sphere

import warnings

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

from pymanopt.manifolds.manifold import EuclideanEmbeddedSubmanifold


class _Sphere(EuclideanEmbeddedSubmanifold):
    """Base class for tensors with unit Frobenius norm."""

    def __init__(self, *shape, name, dimension):
        if len(shape) == 0:
            raise TypeError("Need shape parameters.")
        self._shape = shape
        super().__init__(name, dimension)

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

    def inner(self, X, U, V):
        return float(np.tensordot(U, V, axes=U.ndim))

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

    def dist(self, U, V):
        # Make sure inner product is between -1 and 1
        inner = max(min(self.inner(None, U, V), 1), -1)
        return np.arccos(inner)

    def proj(self, X, H):
        return H - self.inner(None, X, H) * X

    def weingarten(self, X, U, V):
        return -self.inner(X, X, V) * U

    def exp(self, X, U):
        norm_U = self.norm(None, U)
        # Check that norm_U isn't too tiny. If very small then
        # sin(norm_U) / norm_U ~= 1 and retr is extremely close to exp.
        if norm_U > 1e-3:
            return X * np.cos(norm_U) + U * np.sin(norm_U) / norm_U
        else:
            return self.retr(X, U)

    def retr(self, X, U):
        Y = X + U
        return self._normalize(Y)

    def log(self, X, Y):
        P = self.proj(X, Y - X)
        dist = self.dist(X, Y)
        # If the two points are "far apart", correct the norm.
        if dist > 1e-6:
            P *= dist / self.norm(None, P)
        return P

    def rand(self):
        Y = rnd.randn(*self._shape)
        return self._normalize(Y)

    def randvec(self, X):
        H = rnd.randn(*self._shape)
        P = self.proj(X, H)
        return self._normalize(P)

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

    def pairmean(self, X, Y):
        return self._normalize(X + Y)

    def zerovec(self, X):
        return np.zeros(self._shape)

    def _normalize(self, X):
        """
        Return a Frobenius-normalized version of the point X in the ambient
        space.
        """
        return X / self.norm(None, X)


[docs]class Sphere(_Sphere): """Manifold of shape n1 x n2 x ... x nk tensors with unit 2-norm. The metric is such that the sphere is a Riemannian submanifold of Euclidean space. Notes ----- The implementation of the Weingarten map is taken from [1]_. References ---------- .. [1] Absil, P-A., Robert Mahony, and Jochen Trumpf. "An extrinsic look at the Riemannian Hessian." International Conference on Geometric Science of Information. Springer, Berlin, Heidelberg, 2013. """ def __init__(self, *shape): if len(shape) == 0: raise TypeError("Need shape parameters.") if len(shape) == 1: name = "Sphere manifold of {}-vectors".format(*shape) elif len(shape) == 2: name = "Sphere manifold of {}x{} matrices".format(*shape) else: name = "Sphere manifold of shape " + str(shape) + " tensors" dimension = np.prod(shape) - 1 super().__init__(*shape, name=name, dimension=dimension)
class _SphereSubspaceIntersectionManifold(_Sphere): def __init__(self, projector, name, dimension): m, n = projector.shape assert m == n, "projection matrix is not square" if dimension == 0: warnings.warn( "Intersected subspace is 1-dimensional! The manifold '{:s}' " "therefore has dimension 0 as it only consists of isolated " "points".format(self._get_class_name())) self._subspace_projector = projector super().__init__(n, name=name, dimension=dimension) def _validate_span_matrix(self, U): if len(U.shape) != 2: raise ValueError("Input array must be 2-dimensional") num_rows, num_columns = U.shape if num_rows < num_columns: raise ValueError( "The span matrix cannot have fewer rows than columns") def proj(self, X, H): Y = super().proj(X, H) return self._subspace_projector.dot(Y) def rand(self): X = super().rand() return self._normalize(self._subspace_projector.dot(X)) def randvec(self, X): Y = super().randvec(X) return self._normalize(self._subspace_projector.dot(Y))
[docs]class SphereSubspaceIntersection(_SphereSubspaceIntersectionManifold): """Manifold of n-dimensional unit 2-norm vectors intersecting the r-dimensional subspace of R^n spanned by the columns of the matrix U. This implementation is based on spheresubspacefactory.m from the Manopt MATLAB package. """ def __init__(self, U): self._validate_span_matrix(U) m = U.shape[0] Q, _ = la.qr(U) projector = Q.dot(Q.T) subspace_dimension = la.matrix_rank(projector) name = ("Sphere manifold of {}-dimensional vectors intersecting a " "{}-dimensional subspace".format(m, subspace_dimension)) dimension = subspace_dimension - 1 super().__init__(projector, name, dimension)
[docs]class SphereSubspaceComplementIntersection( _SphereSubspaceIntersectionManifold): """Manifold of n-dimensional unit 2-norm vectors which are orthogonal to the r-dimensional subspace of R^n spanned by columns of the matrix U. This implementation is based on spheresubspacefactory.m from the Manopt MATLAB package. """ def __init__(self, U): self._validate_span_matrix(U) m = U.shape[0] Q, _ = la.qr(U) projector = np.eye(m) - Q.dot(Q.T) subspace_dimension = la.matrix_rank(projector) name = ("Sphere manifold of {}-dimensional vectors orthogonal " "to a {}-dimensional subspace".format(m, subspace_dimension)) dimension = subspace_dimension - 1 super().__init__(projector, name, dimension)