// Geodesic distance for Grassmann
def dist(self, X, Y):
if self._k == 1:
u, s, v = np.linalg.svd(np.dot(X.T, Y))
s[s > 1] = 1
s = np.arccos(s)
return np.linalg.norm(s)
else:
XtY = multiprod(multitransp(X), Y)
square_d = 0
for i in xrange(self._k):
s = np.linalg.svd(XtY[i], compute_uv=False)
// Ensure that -1 <= s <= 1
s = np.fmin(s, [1])
s = np.fmax(s, [-1])
After Change
// Geodesic distance for Grassmann
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)