Fix imports in orthogonalize and add unit tests.

Author: tboggs

spectral/algorithms/algorithms.py

@@ -1271,22 +1271,18 @@ def orthogonalize(vecs, start=0):
 
         A new `CxB` containing an orthonormal basis for the given vectors.
     '''
-
     (M, N) = vecs.shape
     basis = np.array(np.transpose(vecs))
-    eye = identity(N).astype(float)
-    if start == 0:
-        basis[:, 0] /= math.sqrt(np.dot(basis[:, 0], basis[:, 0]))
-        start = 1
-
+    eye = np.identity(N).astype(float)
     for i in range(start, M):
-        v = basis[:, i] / math.sqrt(np.dot(basis[:, i], basis[:, i]))
+        if i == 0:
+            basis[:, 0] /= np.linalg.norm(basis[:, 0])
+            continue
+        v = basis[:, i] / np.linalg.norm(basis[:, i])
         U = basis[:, :i]
-        P = eye - np.dot(U, np.dot(np.inv(np.dot(np.transpose(U), U)),
-                                   np.transpose(U)))
-        basis[:, i] = np.dot(P, v)
-        basis[:, i] /= math.sqrt(np.dot(basis[:, i], basis[:, i]))
-
+        P = eye - U.dot(np.linalg.inv(U.T.dot(U)).dot(U.T))
+        basis[:, i] = P.dot(v)
+        basis[:, i] /= np.linalg.norm(basis[:, i])
     return np.transpose(basis)
 
 

spectral/tests/dimensionality.py

@@ -126,6 +126,36 @@ def test_pca_runs_from_stats(self):
         stats = spy.calc_stats(data)
         xdata = spy.principal_components(stats).transform(data)
 
+    def test_orthogonalize(self):
+        '''Can correctly create an orthogonal basis from vectors.'''
+        x = np.linspace(0, np.pi, 1001)
+        # Create sin and cos vectors of unit length
+        sin_h = np.sin(x)
+        sin_h /= np.linalg.norm(sin_h)
+        cos_h = np.cos(x)
+        cos_h /= np.linalg.norm(cos_h)
+
+        X = np.array([50 * sin_h, 75 * cos_h])
+        Y = spy.orthogonalize(X)
+        assert(np.allclose(Y.dot(Y.T), np.array([[1, 0], [0, 1]])))
+        assert(np.allclose(X.dot(Y.T), np.array([[50, 0], [0, 75]])))
+
+    def test_orthogonalize_subset(self):
+        '''Can correctly create an orthogonal basis from vector subset.'''
+        x = np.linspace(0, np.pi, 1001)
+        # Create sin and cos vectors of unit length
+        sin_h = np.sin(x)
+        sin_h /= np.linalg.norm(sin_h)
+        cos_h = np.cos(x)
+        cos_h /= np.linalg.norm(cos_h)
+
+        # First vector in X will already be a unit vector
+        X = np.array([sin_h, 75 * cos_h])
+        Y = spy.orthogonalize(X, start=1)
+        assert(np.allclose(Y.dot(Y.T), np.array([[1, 0], [0, 1]])))
+        assert(np.allclose(X.dot(Y.T), np.array([[1, 0], [0, 75]])))
+
+
 def run():
     print('\n' + '-' * 72)
     print('Running dimensionality tests.')