# Source code for pymanopt.manifolds.stiefel

```import numpy as np
from scipy.linalg import expm

from pymanopt.manifolds.manifold import EuclideanEmbeddedSubmanifold
from pymanopt.tools.multi import multiprod, multisym, multitransp

[docs]class Stiefel(EuclideanEmbeddedSubmanifold):
"""
Factory class for the Stiefel manifold. Instantiation requires the
dimensions n, p to be specified. Optional argument k allows the user to
optimize over the product of k Stiefels.

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!).
"""

def __init__(self, n, p, k=1):
self._n = n
self._p = p
self._k = k

# Check that n is greater than or equal to p
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 = "Stiefel manifold St(%d, %d)" % (n, p)
elif k >= 2:
name = "Product Stiefel manifold St(%d, %d)^%d" % (n, p, k)
dimension = int(k * (n * p - p * (p + 1) / 2))
super().__init__(name, dimension)

@property
def typicaldist(self):
return np.sqrt(self._p * self._k)

[docs]    def inner(self, X, G, H):
# Inner product (Riemannian metric) on the tangent space
# For the stiefel this is the Frobenius inner product.
return np.tensordot(G, H, axes=G.ndim)

[docs]    def proj(self, X, U):
return U - multiprod(X, multisym(multiprod(multitransp(X), U)))

# TODO(nkoep): Implement the weingarten map instead.
[docs]    def ehess2rhess(self, X, egrad, ehess, H):
symXtG = multisym(XtG)
HsymXtG = multiprod(H, symXtG)
return self.proj(X, ehess - HsymXtG)

# Retract to the Stiefel using the qr decomposition of X + G.
[docs]    def retr(self, X, G):
if self._k == 1:
# Calculate 'thin' qr decomposition of X + G
q, r = np.linalg.qr(X + G)
# Unflip any flipped signs
XNew = np.dot(q, np.diag(np.sign(np.sign(np.diag(r)) + 0.5)))
else:
XNew = X + G
for i in range(self._k):
q, r = np.linalg.qr(XNew[i])
XNew[i] = np.dot(
q, np.diag(np.sign(np.sign(np.diag(r)) + 0.5)))
return XNew

[docs]    def norm(self, X, G):
# Norm on the tangent space of the Stiefel is simply the Euclidean
# norm.
return np.linalg.norm(G)

# Generate random Stiefel 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):
# TODO: Simplify these expressions.
if self._k == 1:
W = expm(np.bmat([[X.T.dot(U), -U.T.dot(U)],
[np.eye(self._p), X.T.dot(U)]]))
Z = np.bmat([[expm(-X.T.dot(U))], [np.zeros((self._p, self._p))]])
Y = np.bmat([X, U]).dot(W).dot(Z)
else:
Y = np.zeros(np.shape(X))
for i in range(self._k):
W = expm(np.bmat([[X[i].T.dot(U[i]), -U[i].T.dot(U[i])],
[np.eye(self._p), X[i].T.dot(U[i])]]))
Z = np.bmat([[expm(-X[i].T.dot(U[i]))],
[np.zeros((self._p, self._p))]])
Y[i] = np.bmat([X[i], U[i]]).dot(W).dot(Z)
return Y

[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))
```