From a8ab564b27faf8abea9c1583d73c97b9b36e5d09 Mon Sep 17 00:00:00 2001
From: Jiahao Chen <jiahao@mit.edu>
Date: Fri, 16 Jan 2015 21:39:31 -0500
Subject: [PATCH 1/9] Provide Interpolative factorization for low rank
 approximation of matrices

---
 src/factorization.jl  | 70 +++++++++++++++++++++++++++++++++++++++++++
 test/factorization.jl |  9 ++++++
 2 files changed, 79 insertions(+)
 create mode 100644 src/factorization.jl
 create mode 100644 test/factorization.jl

diff --git a/src/factorization.jl b/src/factorization.jl
new file mode 100644
index 0000000..7d87890
--- /dev/null
+++ b/src/factorization.jl
@@ -0,0 +1,70 @@
+############################
+# Specialized factorizations
+############################
+
+export idfact
+
+immutable Interpolative{T} <: Factorization{T}
+    B :: AbstractMatrix{T}
+    P :: AbstractMatrix{T}
+end
+
+"""
+Compute the interpolative decomposition of A
+
+    A ≈ B * P
+
+    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)
+
+Inputs:
+    A: Matrix to factorize
+    k: Number of columns of A to return in B
+    l: Length of random vectors to project onto
+
+Reference:
+
+    Algorithm I of Liberty2007
+
+Implementation note:
+
+    This is a hacky version of the algorithms described in Liberty2007 and
+    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.
+
+    @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(sub(Y, :, j))), 1:n)
+    piv = sortperm(maxvals, rev=true)[1:k]
+
+    B = A[:, piv]
+    Interpolative(B, B\A)
+end
+
diff --git a/test/factorization.jl b/test/factorization.jl
new file mode 100644
index 0000000..b911be8
--- /dev/null
+++ b/test/factorization.jl
@@ -0,0 +1,9 @@
+using Base.Test
+using IterativeSolvers
+srand(1)
+M = randn(4,5)
+k = 3
+
+F = idfact(M, k, 3)
+vecnorm(F.B*F.P-M) <= 2svdvals(M)[k+1]
+

From 61f918605e1e6a8ed4f1ca4b17e6656e3c02b52b Mon Sep 17 00:00:00 2001
From: Jiahao Chen <jiahao@mit.edu>
Date: Mon, 19 Jan 2015 21:26:53 -0500
Subject: [PATCH 2/9] Update documentation of idfact

---
 src/factorization.jl | 55 ++++++++++++++++++++++++++++----------------
 1 file changed, 35 insertions(+), 20 deletions(-)

diff --git a/src/factorization.jl b/src/factorization.jl
index 7d87890..ad6c4d6 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -4,44 +4,59 @@
 
 export idfact
 
+@doc doc"""
+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}
+""" ->
 immutable Interpolative{T} <: Factorization{T}
     B :: AbstractMatrix{T}
     P :: AbstractMatrix{T}
 end
 
-"""
-Compute the interpolative decomposition of A
+@doc doc"""
+Compute the interpolative decomposition of `A`
 
+    $$
     A ≈ B * P
+    $$
 
-    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)
+    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`.
 
 Inputs:
-    A: Matrix to factorize
-    k: Number of columns of A to return in B
-    l: Length of random vectors to project onto
-
-Reference:
 
-    Algorithm I of Liberty2007
+    `A`: Matrix to factorize
+    `k`: Number of columns of A to return in B
+    `l`: Length of random vectors to project onto
 
 Implementation note:
 
-    This is a hacky version of the algorithms described in Liberty2007 and
-    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.
+    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.
 
+Reference:
+
+    \cite[Algorithm I]{Liberty2007}
+
     @article{Cheng2005,
         author = {Cheng, H and Gimbutas, Z and Martinsson, P G and Rokhlin, V},
         doi = {10.1137/030602678},
@@ -54,7 +69,7 @@ Implementation note:
         volume = {26},
         year = {2005}
     }
-"""
+""" ->
 function idfact(A, k::Int, l::Int)
     m, n = size(A)
     R = randn(l, m)

From ea2a1179129710fb3370ac6d11bd118a6eea30b5 Mon Sep 17 00:00:00 2001
From: Jiahao Chen <jiahao@mit.edu>
Date: Sat, 22 Aug 2015 03:34:27 -0400
Subject: [PATCH 3/9] Use FactCheck for testing framework

---
 test/factorization.jl | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/test/factorization.jl b/test/factorization.jl
index b911be8..ee69006 100644
--- a/test/factorization.jl
+++ b/test/factorization.jl
@@ -1,9 +1,11 @@
-using Base.Test
 using IterativeSolvers
-srand(1)
-M = randn(4,5)
-k = 3
+using FactCheck
 
-F = idfact(M, k, 3)
-vecnorm(F.B*F.P-M) <= 2svdvals(M)[k+1]
+facts("idfact") do
+    srand(1)
+    M = randn(4,5)
+    k = 3
 
+    F = idfact(M, k, 3)
+    @fact vecnorm(F.B*F.P-M) --> less_than_or_equal(2svdvals(M)[k+1])
+end

From ccb86934e82462254631f439cf59b291c98417bc Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juan=20Antonio=20L=C3=B3pez=20Mendoza?=
 <juanlopezm94@gmail.com>
Date: Thu, 18 Aug 2016 22:25:57 +0200
Subject: [PATCH 4/9] Add view instead slice and sub

---
 src/factorization.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/factorization.jl b/src/factorization.jl
index ad6c4d6..024f5ca 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -48,7 +48,7 @@ Implementation note:
     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.
@@ -76,10 +76,9 @@ function idfact(A, k::Int, l::Int)
     Y = R * A #size l x n
 
     #Compute column pivots of first k columns of Y
-    maxvals = map(j->maximum(abs(sub(Y, :, j))), 1:n)
+    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
-

From 4c8ab867a7cc8a32301b7420881cfd09bdfeb347 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Juan=20Antonio=20L=C3=B3pez=20Mendoza?=
 <juanlopezm94@gmail.com>
Date: Fri, 16 Dec 2016 17:19:17 +0100
Subject: [PATCH 5/9] Add Manual (#101)

---
 src/factorization.jl | 116 +++++++++++++++++++++++--------------------
 1 file changed, 61 insertions(+), 55 deletions(-)

diff --git a/src/factorization.jl b/src/factorization.jl
index 024f5ca..4162746 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -4,72 +4,78 @@
 
 export idfact
 
-@doc doc"""
+"""
 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:
+# References
 
-    \cite{Cheng2005, Liberty2007}
-""" ->
+\cite{Cheng2005, Liberty2007}
+"""
 immutable Interpolative{T} <: Factorization{T}
     B :: AbstractMatrix{T}
     P :: AbstractMatrix{T}
 end
 
-@doc doc"""
-Compute 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`.
-
-Inputs:
-
-    `A`: Matrix to factorize
-    `k`: Number of columns of A to return in B
-    `l`: Length of random vectors to project onto
-
-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.
-
-Reference:
-
-    \cite[Algorithm I]{Liberty2007}
-
-    @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}
-    }
-""" ->
+"""
+    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)

From 15facd0ddacf65f25e53dcaa04a2e042c34ea802 Mon Sep 17 00:00:00 2001
From: Harmen Stoppels <harmenstoppels@gmail.com>
Date: Mon, 7 Aug 2017 04:57:54 +0200
Subject: [PATCH 6/9] Remove FactCheck, use Base.test and simplify tests (#150)

* Make sure CG gets the correct types

* Use Base.Test in bicgstab

* Use Base.Test in cg

* Simplify cg

* Chebyshev to base.test and simplified tests

* Remove LinearMaps tests (basically verifying LinearMaps internals) and move to Base.Test

* Simplify GMRES tests and move partially to Base.Test

* Use new API for LinearMap

* Move factorization to Base.Test

* Use Base.Test in Hessenberg test

* Simplify tests and use Base.Test

* Indentation and move initialization to the start of the test

* Use Base.Test in lsmr

* Use Base.Test in orthogonalization tests

* Use Base.Test in rlinalg

* Use Base.Test in rsvd and use broadcasting

* Use Base.Test in rsvd_fnkz

* Use Base.Test in simple eigensolvers and update the API for LinearMaps

* Simplify stationary solvers tests and use Base.Test

* Update the Lanczos test

* Remove FastCheck related things

* Remove FactCheck dependency

* Simplify ls__ tests

* Surpress broadcast warning

* Use broadcasting properly

* Use Base.Test in minres

* Remove empty tests in history smoke test

* Move matrix stuff to benchmark folder

* Reduce dependencies in test REQUIRE file

* Test preconditioners with UmfpackLU on Julia 0.6+ only, since Julia 0.5 does not have support for A_ldiv_B!(::UmfpackLU, ::SubArray)

* Fix builds on Julia nightly by escaping strings properly

* Support 0.5

* Escape even more strings

* Remove LinearMaps where it isn't used

* Remove checking for 0.5

* Use latest API of LinearMaps

* Use I over (sp)eye

* Actually run the orthogonalization test

* Simplify orth. test

* Fix one() in CG
---
 src/factorization.jl  | 10 +++++-----
 test/factorization.jl |  6 +++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/factorization.jl b/src/factorization.jl
index 4162746..4388e67 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -12,7 +12,7 @@ and `P` computed by `idfact()`. See the documentation of `idfact()` for details.
 
 # References
 
-\cite{Cheng2005, Liberty2007}
+\\cite{Cheng2005, Liberty2007}
 """
 immutable Interpolative{T} <: Factorization{T}
     B :: AbstractMatrix{T}
@@ -43,13 +43,13 @@ Where:
 
 # 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
+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
+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.
 
@@ -59,7 +59,7 @@ pivoted QR process.
 
 # References
 
-\cite[Algorithm I]{Liberty2007}
+\\cite[Algorithm I]{Liberty2007}
 
 ```bibtex
 @article{Cheng2005,
diff --git a/test/factorization.jl b/test/factorization.jl
index ee69006..d8af590 100644
--- a/test/factorization.jl
+++ b/test/factorization.jl
@@ -1,11 +1,11 @@
 using IterativeSolvers
-using FactCheck
+using Base.Test
 
-facts("idfact") do
+@testset "IDfact" begin
     srand(1)
     M = randn(4,5)
     k = 3
 
     F = idfact(M, k, 3)
-    @fact vecnorm(F.B*F.P-M) --> less_than_or_equal(2svdvals(M)[k+1])
+    @test vecnorm(F.B * F.P - M) ≤ 2svdvals(M)[k + 1]
 end

From 0a07c531a9208061895caaa3c030226077c4015b Mon Sep 17 00:00:00 2001
From: Harmen Stoppels <harmenstoppels@gmail.com>
Date: Tue, 8 Aug 2017 01:15:06 +0200
Subject: [PATCH 7/9] Modernize syntax (#153)

* Modernize w/ struct and where

* Clean up some imports

* Update deprecated / removed exception type and add a test for it

* Remove compat

* Add row number to SingularException

* @fact_throws -> @test_throws
---
 src/factorization.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/factorization.jl b/src/factorization.jl
index 4388e67..c370915 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -14,7 +14,7 @@ and `P` computed by `idfact()`. See the documentation of `idfact()` for details.
 
 \\cite{Cheng2005, Liberty2007}
 """
-immutable Interpolative{T} <: Factorization{T}
+struct Interpolative{T} <: Factorization{T}
     B :: AbstractMatrix{T}
     P :: AbstractMatrix{T}
 end

From 4d371472391a8ffa4e7164a3b259532c6220f5cd Mon Sep 17 00:00:00 2001
From: Harmen Stoppels <harmenstoppels@gmail.com>
Date: Sat, 9 Dec 2017 12:40:23 +0100
Subject: [PATCH 8/9] Add factorization

---
 src/RandomizedLinAlg.jl | 1 +
 src/factorization.jl    | 2 +-
 test/factorization.jl   | 2 +-
 test/runtests.jl        | 1 +
 4 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/RandomizedLinAlg.jl b/src/RandomizedLinAlg.jl
index d919d19..d1daf75 100644
--- a/src/RandomizedLinAlg.jl
+++ b/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")
diff --git a/src/factorization.jl b/src/factorization.jl
index c370915..0bb29a2 100644
--- a/src/factorization.jl
+++ b/src/factorization.jl
@@ -82,7 +82,7 @@ function idfact(A, k::Int, l::Int)
     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)
+    maxvals = map(j->maximum(abs.(view(Y, :, j))), 1:n)
     piv = sortperm(maxvals, rev=true)[1:k]
 
     B = A[:, piv]
diff --git a/test/factorization.jl b/test/factorization.jl
index d8af590..3aed744 100644
--- a/test/factorization.jl
+++ b/test/factorization.jl
@@ -1,4 +1,4 @@
-using IterativeSolvers
+using RandomizedLinAlg
 using Base.Test
 
 @testset "IDfact" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 7ced9a0..0deef02 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using RandomizedLinAlg
 
+include("factorization.jl")
 include("rlinalg.jl")
 include("rsvd.jl")
 include("rsvd_fnkz.jl")
\ No newline at end of file

From f646311e5ca4d25462a370a5865a23ac4faadb75 Mon Sep 17 00:00:00 2001
From: Harmen Stoppels <harmenstoppels@gmail.com>
Date: Sat, 9 Dec 2017 12:43:53 +0100
Subject: [PATCH 9/9] Add docs for id fact

---
 docs/src/index.md | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/docs/src/index.md b/docs/src/index.md
index 078ad28..18ea391 100644
--- a/docs/src/index.md
+++ b/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.