Added support for sparsity in the console application.

irlb/scipy_bench.py

@@ -1,8 +1,10 @@
 from __future__ import division, print_function
 from irlb import *
 from argparse import *
+import scipy.sparse as sp
+import scipy.sparse.linalg as splinalg
 
-import timeit
+import time
 
 def parse_args():
   parser = ArgumentParser()
@@ -21,29 +23,62 @@ def parse_args():
     default=False,
     help="Should the output be csv formatted? If True then the return will "+
       "be: <IRB timing>,<SVD timing>,<absolue error>")
+  parser.add_argument("-s", "--sparsity", default=0., type=float,
+    help="The approximate sparsity of the matrix. Default is 0 (dense).")
+  parser.add_argument("-l", "--irlb-only", action="store_const", const=True,
+    default=False)
+  parser.add_argument("-p", "--scipy-only", action="store_const", const=True,
+    default=False)
+
   return parser.parse_args()
 
 def main():
   args = parse_args()
-  A=np.random.rand(args.num_row*args.num_col).reshape(args.num_row,args.num_col)
+  if not args.scipy_only:
+    if args.sparsity > 0.:
+      if args.sparsity >= 1.:
+        raise(RuntimeError("Sparsity must be less than 1."))
+      nnz = (1.-args.sparsity) * args.num_row * args.num_col
+      v = np.random.rand(nnz)
+      i = np.random.randint(0, high=args.num_row, size=nnz)
+      j = np.random.randint(0, high=args.num_col, size=nnz)
+      A = sp.csc_matrix((v,(i,j)), shape=(args.num_row, args.num_col))
+    else:
+      A=np.random.rand(args.num_row*args.num_col).reshape(
+        args.num_row,args.num_col)
 
-  start_irlb = timeit.timeit()
-  X = irlb(A, args.nu, args.tol)
-  end_irlb = timeit.timeit()
+    start_irlb = time.time()
+    X = irlb(A, args.nu, tol=args.tol, maxit=args.max_it)
+    end_irlb = time.time()
 
   # Compare estimated values with np.linalg.svd:
-  start_svd = timeit.timeit()
-  S = np.linalg.svd(A, 0)
-  end_svd = timeit.timeit()
+  if not args.irlb_only:
+    if args.sparsity > 0.:
+      start_svd = time.time()
+      S = splinalg.svds(A=A, k=args.nu, tol=args.tol, maxiter=args.max_it,
+        return_singular_vectors=True)
+      end_svd = time.time()
+    else:
+      start_svd = time.time()
+      S = np.linalg.svd(A, 0)
+      end_svd = time.time()
 
-  abserr = np.max(np.abs(X[1]-S[1][0:args.nu]))
+    abserr = np.max(np.abs(X[1]-S[1][0:args.nu]))
 
   if (args.csv):
-    print("{},{},{}".format(end_irlb-start_irlb, end_svd-start_svd, abserr))
+    if args.irlb_only:
+      print("{}".format(end_irlb-start_irlb))
+    elif args.scipy_only:
+      print("{}".format(end_svd-start_svd))
+    else:
+      print("{},{},{}".format(end_irlb-start_irlb, end_svd-start_svd, abserr))
   else:
-    print("IRLB timing: {}".format(end_irlb-start_irlb))
-    print("SVD timing: {}".format(end_svd-start_svd))
-    print("Absolute error: {}".format(abserr))
+    if not args.scipy_only:
+      print("IRLB timing: {}".format(end_irlb-start_irlb))
+    if not args.irlb_only:
+      print("SVD timing: {}".format(end_svd-start_svd))
+    if not args.irlb_only and not args.scipy_only:
+      print("Absolute error: {}".format(abserr))
 
 if __name__=="__main__":
   main()

tests/test_sparse.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 import scipy.sparse as sp
+from irlb import *
 
 def main():
   m   = 20000   # Number of rows
@@ -13,8 +14,11 @@ def main():
 
   v   = np.random.rand(nnz)
   i   = np.random.randint(0,high=m,size=nnz)
-  j   = np.random.randint(0,high=m,size=nnz)
+  j   = np.random.randint(0,high=n,size=nnz)
   A   = sp.csc_matrix((v,(i,j)),shape=(m,n))
   X = irlb(A,nu,tol=0.000001)
   # A gut check measurement of absolute decomposition error:
   print("TSVD: ||AV - US||_F = %f" % np.linalg.norm(A.dot(sp.csr_matrix(X[2])).todense().A - X[0]*X[1],"fro"))
+
+  S = sp.linalg.svds(A, nu)
+  abserr = np.max(np.abs(X[1]-S[1][0:nu]))