Merge pull request #1 from haampie/add-factorization

Add interpolative decomposition

docs/src/index.md

@@ -28,6 +28,12 @@ rnorm
 rnorms
 ```
 
+## Interpolative Decomposition
+
+```@docs
+idfact
+```
+
 [^Halko2011]: Halko, Nathan, Per-Gunnar Martinsson, and Joel A. Tropp. "Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions." SIAM review 53.2 (2011): 217-288.
 
 [^Dixon1983]: Dixon, John D. "Estimating extremal eigenvalues and condition numbers of matrices." SIAM Journal on Numerical Analysis 20.4 (1983): 812-814.

src/RandomizedLinAlg.jl

@@ -5,6 +5,7 @@ for randomized methods in numerical linear algebra.
 """
 module RandomizedLinAlg
 
+include("factorization.jl")
 include("rlinalg.jl")
 include("rsvd.jl")
 include("rsvd_fnkz.jl")

src/factorization.jl

@@ -0,0 +1,90 @@
+############################
+# Specialized factorizations
+############################
+
+export idfact
+
+"""
+An interpolative decomposition.
+
+For a matrix `A`, the interpolative decomposition `F` contains the matrices `B`
+and `P` computed by `idfact()`. See the documentation of `idfact()` for details.
+
+# References
+
+\\cite{Cheng2005, Liberty2007}
+"""
+struct Interpolative{T} <: Factorization{T}
+    B :: AbstractMatrix{T}
+    P :: AbstractMatrix{T}
+end
+
+"""
+    idfact(A, k, l)
+
+Compute and return the interpolative decomposition of `A`: A ≈ B * P
+
+Where:
+* `B`'s columns are a subset of the columns of `A`
+* some subset of `P`'s columns are the `k x k` identity, no entry of `P` exceeds magnitude 2, and
+* ||B * P - A|| ≲ σ(A, k+1), the (`k+1`)st singular value of `A`.
+
+# Arguments
+
+`A`: Matrix to factorize
+
+`k::Int`: Number of columns of A to return in B
+
+`l::Int`: Length of random vectors to project onto
+
+# Output
+
+`(::Interpolative)`: interpolative decomposition.
+
+# Implementation note
+
+This is a hacky version of the algorithms described in \\cite{Liberty2007}
+and \\cite{Cheng2005}. The former refers to the factorization (3.1) of the
+latter.  However, it is not actually necessary to compute this
+factorization in its entirely to compute an interpolative decomposition.
+
+Instead, it suffices to find some permutation of the first k columns of Y =
+R * A, extract the subset of A into B, then compute the P matrix as B\\A
+which will automatically compute P using a suitable least-squares
+algorithm.
+
+The approximation we use here is to compute the column pivots of Y,
+rather then use the true column pivots as would be computed by a column-
+pivoted QR process.
+
+# References
+
+\\cite[Algorithm I]{Liberty2007}
+
+```bibtex
+@article{Cheng2005,
+    author = {Cheng, H and Gimbutas, Z and Martinsson, P G and Rokhlin, V},
+    doi = {10.1137/030602678},
+    issn = {1064-8275},
+    journal = {SIAM Journal on Scientific Computing},
+    month = jan,
+    number = {4},
+    pages = {1389--1404},
+    title = {On the Compression of Low Rank Matrices},
+    volume = {26},
+    year = {2005}
+}
+```
+"""
+function idfact(A, k::Int, l::Int)
+    m, n = size(A)
+    R = randn(l, m)
+    Y = R * A #size l x n
+
+    #Compute column pivots of first k columns of Y
+    maxvals = map(j->maximum(abs.(view(Y, :, j))), 1:n)
+    piv = sortperm(maxvals, rev=true)[1:k]
+
+    B = A[:, piv]
+    Interpolative(B, B\A)
+end

test/factorization.jl

@@ -0,0 +1,11 @@
+using RandomizedLinAlg
+using Base.Test
+
+@testset "IDfact" begin
+    srand(1)
+    M = randn(4,5)
+    k = 3
+
+    F = idfact(M, k, 3)
+    @test vecnorm(F.B * F.P - M) ≤ 2svdvals(M)[k + 1]
+end

test/runtests.jl

@@ -1,5 +1,6 @@
 using RandomizedLinAlg
 
+include("factorization.jl")
 include("rlinalg.jl")
 include("rsvd.jl")
 include("rsvd_fnkz.jl")