Added truncated SVD (aimacode#587)

Author: Chipe1

tests/test_utils.py

@@ -139,6 +139,12 @@ def test_normalize():
     assert normalize([1, 2, 1]) == [0.25, 0.5, 0.25]
 
 
+def test_norm():
+    assert isclose(norm([1, 2, 1], 1), 4)
+    assert isclose(norm([3, 4], 2), 5)
+    assert isclose(norm([-1, 1, 2], 4), 18**0.25)
+
+
 def test_clip():
     assert [clip(x, 0, 1) for x in [-1, 0.5, 10]] == [0, 0.5, 1]
 
@@ -155,6 +161,25 @@ def test_gaussian():
     assert gaussian(3,1,3) == 0.3989422804014327
 
 
+def test_truncated_svd():
+    test_mat = [[17, 0],
+                [0, 11]]
+    _, _, eival = truncated_svd(test_mat)
+    assert isclose(eival, 17)
+
+    test_mat = [[17, 0],
+                [0, -34]]
+    _, _, eival = truncated_svd(test_mat)
+    assert isclose(eival, -34)
+
+    test_mat = [[1, 0, 0, 0, 2],
+                [0, 0, 3, 0, 0],
+                [0, 0, 0, 0, 0],
+                [0, 2, 0, 0, 0]]
+    _, _, eival = truncated_svd(test_mat)
+    assert isclose(eival, 3)
+
+
 def test_sigmoid_derivative():
     value = 1
     assert sigmoid_derivative(value) == 0

utils.py

@@ -246,6 +246,11 @@ def normalize(dist):
     return [(n / total) for n in dist]
 
 
+def norm(X, n=2):
+    """Return the n-norm of vector X"""
+    return sum([x**n for x in X])**(1/n)
+
+
 def clip(x, lowest, highest):
     """Return x clipped to the range [lowest..highest]."""
     return max(lowest, min(x, highest))
@@ -270,6 +275,41 @@ def gaussian(mean, st_dev, x):
     return 1/(math.sqrt(2*math.pi)*st_dev)*math.e**(-0.5*(float(x-mean)/st_dev)**2)
 
 
+def truncated_svd(X, max_iter=1000):
+    """Computes the first component of SVD"""
+
+    def normalize_vec(X, n = 2):
+        """Returns normalized vector"""
+        norm_X = norm(X, n)
+        Y = [x/norm_X for x in X]
+        return Y
+
+    m, n = len(X), len(X[0])
+    A = [[0 for _ in range(n + m)] for _ in range(n + m)]
+    for i in range(m):
+        for j in range(n):
+            A[i][j] = A[m + j][n + i] = X[i][j]
+
+    X = [random.random() for _ in range(n + m)]
+    X = normalize_vec(X)
+    for _ in range(max_iter):
+        old_X = X
+        X = matrix_multiplication(A, [[x] for x in X])
+        X = [x[0] for x in X]
+        X = normalize_vec(X)
+        # check for convergence
+        if norm([x1 - x2 for x1, x2 in zip(old_X, X)]) <= 1e-10:
+            break
+
+    projected_X = matrix_multiplication(A, [[x] for x in X])
+    projected_X = [x[0] for x in projected_X]
+    eival = norm(projected_X, 1)/norm(X, 1)
+    eivec_n = normalize_vec(X[:n])
+    eivec_m = normalize_vec(X[n:])
+
+    return (eivec_m, eivec_n, eival)
+
+
 try:  # math.isclose was added in Python 3.5; but we might be in 3.4
     from math import isclose
 except ImportError: