Renamed eigSVD, fixed types, exact_svd parameter

Author: Viktor Petukhov

Project.toml

@@ -4,14 +4,14 @@ authors = ["Viktor Petukhov <viktor.petukhov+gh@pm.me> and contributors"]
 version = "0.1.0"
 
 [deps]
-IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-IterativeSolvers = "0.9"
 julia = "1"
+KrylovKit = "0.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -4,7 +4,7 @@
 
 Fast randomized PCA algorithms for Sparse Data. It's a julia re-implementation of Matlab [frPCA_sparse](https://github.com/XuFengthucs/frPCA_sparse).
 
-The package includes an implementation of two functions `eigSVD` and `pca`, which are designed to work with sparse matrices (`SparseMatrixCSC`), though can also be applied with normal ones (`Matrix`).
+The package includes an implementation of two functions `eig_svd` and `pca`, which are designed to work with sparse matrices (`SparseMatrixCSC`), though can also be applied with normal ones (`Matrix`).
 
 ## References
 

src/FastRandPCA.jl

@@ -1,32 +1,39 @@
 module FastRandPCA
 
-import IterativeSolvers: lobpcg
+using KrylovKit: svdsolve, KrylovDefaults
 
 using SparseArrays
 using LinearAlgebra
 using Random
 
-export eigSVD, pca
+export eig_svd, pca
 
 """
-Perform truncated singular value decomposition of a matrix `A` to return `k` first values.
+`eig_svd(A, [k]; [kwargs...])`
 
-Designed to work with sparse matrices. If `k` is not specified, all values are returned.
+Perform singular value decomposition of a matrix `A` by estimating `svd(A' * A)`. If `k` is provided,
+Krylov iterative solver is used to return only `k` largest singular values and corresponding singular vectors.
+
+Designed to work with sparse matrices. Keyword arguments are forwarded to `KrylovKit.svdsolve`.
 
 # Examples
+
 ```julia-repl
 using SparseArrays
-eigSVD(sprand(100, 100, 0.2), 10)
+eig_svd(sprand(100, 100, 0.2), 10)
 ```
 """
-function eigSVD(A::AbstractMatrix{<:Number}, k::Int=-1)
+function eig_svd(A::AbstractMatrix{<:Number}, k::Int=-1; krylovdim::Int=-1, kwargs...)
     B = A' * A;
-    if k == -1
+    if (k == -1) || (k == size(B, 1))
         res = eigen(Matrix(B))
         V, λ = res.vectors, res.values
     else
-        r = lobpcg(B, true, k)
-        V, λ = r.X, r.λ
+        if krylovdim <= 0
+            krylovdim = max(KrylovDefaults.krylovdim, k + 2)
+        end
+        r = svdsolve(B, size(B, 1), k, krylovdim=krylovdim, kwargs...)
+        V, λ = hcat(r[2]...), r[1]
     end
 
     λ[λ .< 0] .= 0 # Some of zero values become negative due to numerical errors
@@ -36,49 +43,59 @@ function eigSVD(A::AbstractMatrix{<:Number}, k::Int=-1)
     return SVD(U, d, V)
 end
 
+svd_wrap(A::AbstractMatrix{<:Number}, exact_svd::Bool) = (exact_svd ? svd(A) : eig_svd(A))
+
 """
+`pca(A, k; [q=10], [exact_svd=false], [s=5])`
+
 Perform fast randomized PCA of a matrix `A` to return `k` first components.
 
 Designed to work with sparse matrices.
 
-Parameter `q` is the number of pass over `A`, `q` should larger than 1 and `q` times pass eqauls to `(q-2)/2` times power iteration.
+# Keyword arguments
+
+- `exact_svd` - if `true`, use exact SVD on random projections, otherwise use `eig_svd` approximation.
+  Setting it to `true` would improve precision of the algorithm on small eigenvalues, but would drastically decrease performance
+  on large matrices (>1M rows).
+- `q` is the number of pass over `A`, `q` should larger than 1 and `q` times pass eqauls to `(q-2)/2` times power iteration
+- `s` is the excess dimension for the random matrix. This parameter is not supposed to be changed normally.
 
 # Examples
+
 ```julia-repl
 using SparseArrays
 pca(sprand(100, 100, 0.2), 10; q=3)
 ```
 """
-function pca(A::AbstractMatrix{<:Number}, k::Int; q::Int=10)
+function pca(A::AbstractMatrix{<:Number}, k::Int; q::Int=10, exact_svd::Bool=false, s=5)
     q >= 2 || error("Pass parameter q must be larger than 1 !");
 
-    s = 5;
     if size(A, 1) > size(A, 2)
         r = pca(A', k; q=q)
         return SVD(Matrix(r.Vt), r.S, Matrix(r.U))
     end
 
     m,n = size(A);
     if q % 2 == 0
-        Q = randn(n, k+s);
+        Q = randn(promote_type(Float32, eltype(A)), n, k+s);
         Q = A*Q;
         if q == 2
-            Q = eigSVD(Q).U;
+            Q = svd_wrap(Q, exact_svd).U;
         else
             Q = lu(Q).L;
         end
     else
-        Q = randn(m, k+s);
+        Q = randn(promote_type(Float32, eltype(A)), m, k+s);
     end
     upper = floor((q-1)/2);
     for i = 1:upper
         if i == upper
-            Q = eigSVD(A*(A'*Q)).U;
+            Q = svd_wrap(A*(A'*Q), exact_svd).U;
         else
             Q = lu(A*(A'*Q)).L;
         end
     end
-    V,S,U = eigSVD(A'*Q);
+    V,S,U = svd_wrap(A'*Q, exact_svd);
     ind = s+1:k+s;
     U = (Q * U[:, ind])';
     V = V[:, ind]';

test/runtests.jl

@@ -6,18 +6,20 @@ using SparseArrays
     for (n,m) in [[100, 20], [20, 100], [1000, 1000]]
         mat = sprand(n, m, 0.2)
 
-        eigs = eigSVD(mat)
+        eigs = eig_svd(mat)
         @test all(size(eigs.V) .== size(mat, 2))
         @test length(eigs.S) == size(mat, 2)
         @test all(size(eigs.U) .== size(mat))
 
         k = 10
-        pcs = pca(mat, k; q=3)
-        @test size(pcs.Vt, 2) == size(mat, 2)
-        @test size(pcs.Vt, 1) == k
-        @test length(pcs.S) == k
-        @test size(pcs.U, 2) == size(mat, 1)
-        @test size(pcs.U, 1) == k
-        @test all(size(pcs.U * mat) .== size(pcs.Vt))
+        for exact in [false, true]
+            pcs = pca(mat, k; q=3, exact_svd=exact)
+            @test size(pcs.Vt, 2) == size(mat, 2)
+            @test size(pcs.Vt, 1) == k
+            @test length(pcs.S) == k
+            @test size(pcs.U, 2) == size(mat, 1)
+            @test size(pcs.U, 1) == k
+            @test all(size(pcs.U * mat) .== size(pcs.Vt))
+        end
     end
 end