Source code for pymanopt.manifolds.grassmann

import numpy as np
from numpy.linalg import svd

from pymanopt.manifolds.manifold import Manifold
from import multiprod, multitransp

[docs]class Grassmann(Manifold): """ Factory class for the Grassmann manifold. This is the manifold of p- dimensional subspaces of n dimensional real vector space. Initiation requires the dimensions n, p to be specified. Optional argument k allows the user to optimize over the product of k Grassmanns. Elements are represented as n x p matrices (if k == 1), and as k x n x p matrices if k > 1 (Note that this is different to manopt!). """ # I have chaned the retraction to one using the polar decomp as am now # implementing vector transport. See comment below (JT) # April 17, 2013 (NB) : # Retraction changed to the polar decomposition, so that the vector # transport is now correct, in the sense that it is compatible with # the retraction, i.e., transporting a tangent vector G from U to V # where V = Retr(U, H) will give Z, and transporting GQ from UQ to VQ # will give ZQ: there is no dependence on the representation, which # is as it should be. Notice that the polar factorization requires an # SVD whereas the qfactor retraction requires a QR decomposition, # which is cheaper. Hence, if the retraction happens to be a # bottleneck in your application and you are not using vector # transports, you may want to replace the retraction with a qfactor. def __init__(self, n, p, k=1): self._n = n self._p = p self._k = k if n < p or p < 1: raise ValueError("Need n >= p >= 1. Values supplied were n = %d " "and p = %d." % (n, p)) if k < 1: raise ValueError("Need k >= 1. Value supplied was k = %d." % k) if k == 1: name = "Grassmann manifold Gr({:d}, {:d})".format(n, p) elif k >= 2: name = "Product Grassmann manifold Gr({:d}, {:d})^{:d}".format( n, p, k) dimension = int(k * (n * p - p ** 2)) super().__init__(name, dimension) @property def typicaldist(self): return np.sqrt(self._p * self._k) # Geodesic distance for Grassmann
[docs] def dist(self, X, Y): u, s, v = svd(multiprod(multitransp(X), Y)) s[s > 1] = 1 s = np.arccos(s) return np.linalg.norm(s)
[docs] def inner(self, X, G, H): # Inner product (Riemannian metric) on the tangent space # For the Grassmann this is the Frobenius inner product. return np.tensordot(G, H, axes=G.ndim)
[docs] def proj(self, X, U): return U - multiprod(X, multiprod(multitransp(X), U))
egrad2rgrad = proj
[docs] def ehess2rhess(self, X, egrad, ehess, H): # Convert Euclidean into Riemannian Hessian. PXehess = self.proj(X, ehess) XtG = multiprod(multitransp(X), egrad) HXtG = multiprod(H, XtG) return PXehess - HXtG
[docs] def retr(self, X, G): # Calculate 'thin' qr decomposition of X + G # XNew, r = np.linalg.qr(X + G) # We do not need to worry about flipping signs of columns here, # since only the column space is important, not the actual # columns. Compare this with the Stiefel manifold. # Compute the polar factorization of Y = X+G u, s, vt = svd(X + G, full_matrices=False) return multiprod(u, vt)
[docs] def norm(self, X, G): # Norm on the tangent space is simply the Euclidean norm. return np.linalg.norm(G)
# Generate random Grassmann point using qr of random normally distributed # matrix.
[docs] def rand(self): if self._k == 1: X = np.random.randn(self._n, self._p) q, r = np.linalg.qr(X) return q X = np.zeros((self._k, self._n, self._p)) for i in range(self._k): X[i], r = np.linalg.qr(np.random.randn(self._n, self._p)) return X
[docs] def randvec(self, X): U = np.random.randn(*np.shape(X)) U = self.proj(X, U) U = U / np.linalg.norm(U) return U
[docs] def transp(self, x1, x2, d): return self.proj(x2, d)
[docs] def exp(self, X, U): u, s, vt = svd(U, full_matrices=False) cos_s = np.expand_dims(np.cos(s), -2) sin_s = np.expand_dims(np.sin(s), -2) Y = (multiprod(multiprod(X, multitransp(vt) * cos_s), vt) + multiprod(u * sin_s, vt)) # From numerical experiments, it seems necessary to # re-orthonormalize. This is overall quite expensive. if self._k == 1: Y, unused = np.linalg.qr(Y) return Y else: for i in range(self._k): Y[i], unused = np.linalg.qr(Y[i]) return Y
[docs] def log(self, X, Y): ytx = multiprod(multitransp(Y), X) At = multitransp(Y) - multiprod(ytx, multitransp(X)) Bt = np.linalg.solve(ytx, At) u, s, vt = svd(multitransp(Bt), full_matrices=False) arctan_s = np.expand_dims(np.arctan(s), -2) U = multiprod(u * arctan_s, vt) return U
[docs] def zerovec(self, X): if self._k == 1: return np.zeros((self._n, self._p)) return np.zeros((self._k, self._n, self._p))