Fast Hadamard Transform

Project.toml

@@ -6,9 +6,11 @@ version = "0.1.9"
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
+Hadamard = "4a05ff16-5f95-55f4-bb53-bb3f467c689a"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

docs/src/api/helper_functions.md

@@ -0,0 +1,10 @@
+## Internal Helper Functions
+
+```@contents
+Pages = ["helper_functions.md"]
+```
+## Fast Hadamard transform
+
+```@docs
+RLinearAlgebra.fwht!
+```

docs/src/api/linear_samplers.md

@@ -49,6 +49,13 @@ LinSysVecRowMaxDistance
 
 LinSysVecRowDistCyclic
 ```
+## Block Vector Row Samplers
+
+```@docs
+LinSysBlockRowFJLT
+
+LinSysBlockRowSRHT
+```
 
 ## Vector Column Samplers
 
@@ -58,7 +65,7 @@ LinSysVecColDetermCyclic
 LinSysVecColOneRandCyclic
 ```
 
-## Block Vector Row Samplers
+## Block Row Samplers
 
 ```@docs
 LinSysBlkRowGaussSampler
@@ -68,14 +75,19 @@ LinSysBlkRowReplace
 LinSysBlkRowRandCyclic
 ```
 
-## Block Vector Col Samplers
+## Block Col Samplers
 
 ```@docs
 LinSysBlkColGaussSampler
 
 LinSysBlkColReplace
 
 LinSysBlkColRandCyclic
+
+LinSysBlockColFJLT
+
+LinSysBlockColSRHT
+
 ```
 
 ## Sample Function

docs/src/man/cls_overview.md

@@ -170,3 +170,8 @@ solver = RLSSolver(
 # Solve the system
 sol = rsolve(solver, A, b)
 ```
+
+## Block Methods for Linear Systems
+
+Because of the way that computers operate, it is often more efficient to work using 
+blocks of data rather than single vectors to generate updates to solutions. 

src/RLinearAlgebra.jl

@@ -23,6 +23,9 @@ module RLinearAlgebra
 
 using LinearAlgebra, Random, Distributions
 
+import SparseArrays: sprandn, SparseMatrixCSC
+
+import Hadamard: hadamard
 ###########################################
 # Exports
 ###########################################
@@ -43,14 +46,16 @@ export LinSysVecRowDetermCyclic, LinSysVecRowHopRandCyclic, LinSysVecRowOneRandC
     LinSysVecRowPropToNormSampler, LinSysVecRowSVSampler, LinSysVecRowRandCyclic,
     LinSysVecRowUnidSampler, LinSysVecRowUnifSampler, LinSysVecRowGaussSampler,
     LinSysVecRowSparseUnifSampler, LinSysVecRowSparseGaussSampler, LinSysVecRowMaxResidual,
-    LinSysVecRowMaxDistance, LinSysVecRowResidCyclic, LinSysVecRowDistCyclic
+    LinSysVecRowMaxDistance, LinSysVecRowResidCyclic, LinSysVecRowDistCyclic, LinSysBlockRowSRHT,
+    LinSysBlockRowFJLT
 
 # Vector Column Samplers
 export LinSysVecColDetermCyclic, LinSysVecColOneRandCyclic
 #Vector Block Row Samplers
 export LinSysBlkRowGaussSampler, LinSysBlkRowRandCyclic, LinSysBlkRowReplace
 #Vector Block Column Samplers
-export LinSysBlkColRandCyclic, LinSysBlkColGaussSampler, LinSysBlkColReplace
+export LinSysBlkColRandCyclic, LinSysBlkColGaussSampler, LinSysBlkColReplace, LinSysBlockColSRHT, LinSysBlockColFJLT
+
 #*****************************************#
 # Linear Solver Routine Exports
 #*****************************************#

src/linear_samplers.jl

@@ -141,6 +141,8 @@ include("linear_samplers/vec_row_uniform.jl")
 include("linear_samplers/vec_row_gaussian.jl")
 include("linear_samplers/vec_row_uniform_sparse.jl")
 include("linear_samplers/vec_row_gaussian_sparse.jl")
+include("linear_samplers/block_row_SRHT.jl")
+include("linear_samplers/block_row_FJLT.jl")
 #include("linear_samplers/vec_row_uniform_sym_sparse.jl")
 #include("linear_samplers/vec_row_uniform_sym.jl")
 
@@ -161,6 +163,13 @@ include("linear_samplers/vec_col_one_rand_cyclic.jl")
 # Non-adaptive Sampling (with replacement)
 #Leventhal-Lewis (non-symmetric)
 
+# Non-adaptive Sketching
+include("linear_samplers/block_col_SRHT.jl")
+include("linear_samplers/block_col_FJLT.jl")
+#############################################
+# Fast Transforms
+#############################################
+include("linear_samplers/sampler_helpers/hadamard.jl")
 #############################################
 # Block Row Sampler/Sketch/Selector
 #############################################

src/linear_samplers/block_col_FJLT.jl

@@ -0,0 +1,83 @@
+"""
+    LinSysBlockColFJLT <: LinSysBlkColSampler
+
+A mutable structure with fields to handle FJLT row sketching. For this procedure,
+the hadamard transform and random sign swaps are applied once, then that matrix is repeatably
+sampled.
+
+# Fields
+- `blockSize::Int64`, the size of the sketching dimension
+- `sparsity::Float64`, the sparsity of the sampling matrix
+- `paddedSize::Int64`, the size of the matrix when padded
+- `Sketch::Union{SparseMatrixCSC, Nothing}`, storage for sparse sketching matrix 
+- `Ap::Union{AbstractMatrix, Nothing}`, storage for padded matrix
+- `bp::Union{AbstractMatrix, Nothing}`, storage for padded vector
+- `signs::Union{Vector{Bool}, Nothing}`, storage for random sign flips.
+- `scaling::Float64`, storage for the scaling of the sketches.
+
+Calling `LinSysBlockColFJLT()` defaults to setting `sparsity` to .3 and the blocksize to 2.
+
+Nir Ailon and Bernard Chazelle. 2006. Approximate nearest neighbors and the fast Johnson-Lindenstrauss transform. In Proceedings of the thirty-eighth annual ACM symposium on Theory of Computing (STOC '06). Association for Computing Machinery, New York, NY, USA, 557–563. https://doi.org/10.1145/1132516.1132597
+"""
+mutable struct LinSysBlockColFJLT <: LinSysBlkColSampler
+    blockSize::Int64
+    sparsity::Float64 
+    paddedSize::Int64
+    Sketch::Union{SparseMatrixCSC, Nothing}
+    Ap::Union{AbstractMatrix, Nothing}
+    bp::Union{AbstractVector, Nothing}
+    signs::Union{Vector{Bool}, Nothing}
+    scaling::Float64
+end
+
+LinSysBlockColFJLT(;blocksize = 2, sparsity = .3) = LinSysBlockColFJLT(
+                                                   blocksize,
+                                                   sparsity, 
+                                                   0, 
+                                                   nothing, 
+                                                   nothing,
+                                                   nothing,
+                                                   nothing,
+                                                   0.0
+                                                  )
+
+# Common sample interface for linear systems
+function sample(
+    type::LinSysBlockColFJLT,
+    A::AbstractArray,
+    b::AbstractVector,
+    x::AbstractVector,
+    iter::Int64
+)
+    if iter == 1
+        m, n = size(A)
+        # If matrix is not a power of 2 then pad the rows
+        if rem(log(2, n), 1) != 0
+            type.paddedSize = Int64(2^(div(log(2, n), 1) + 1)) 
+            # Find nearest power 2 and allocate
+            type.Ap = zeros(m, type.paddedSize)
+            # Pad matrix and constant vector
+            @views type.Ap[:, 1:n] .= A
+        else
+            type.paddedSize = n
+            type.Ap = A
+        end
+        # Compute scaling and sign flips
+        type.scaling = sqrt(type.blockSize / (type.paddedSize * type.sparsity))
+        type.signs = bitrand(type.paddedSize)
+        # Apply FWHT to padded matrix and vector
+        for i = 1:m
+            @views fwht!(type.Ap[i, :], signs = type.signs, scaling = type.scaling)
+        end
+
+    end
+
+    type.Sketch = sprandn(type.paddedSize, type.blockSize, type.sparsity) 
+    AS = type.Ap * type.Sketch
+    # Residual of the linear system
+    res = A * x - b
+    grad = AS'res
+    H = hadamard(type.paddedSize)
+    sgn = [type.signs[i] ? 1 : -1 for i in 1:type.paddedSize]
+    return [sgn, type.Sketch .* type.scaling], AS, res, grad
+end

src/linear_samplers/block_col_SRHT.jl

@@ -0,0 +1,79 @@
+
+"""
+    LinSysBlockColSRHT <: LinSysBlkColSampler
+
+A mutable structure with fields to handle SRHT column sketching. For this procedure,
+the hadamard transform and random sign swaps are applied once, then that matrix is repeatably
+sampled.
+
+# Fields
+- `blockSize::Int64`, the size of blocks being chosen
+- `paddedSize::Int64`, the size of the matrix when padded
+- `block::Union{Vector{Int64}, Nothing}`, storage for block indices
+- `Ap::Union{AbstractMatrix, Nothing}`, storage for padded matrix
+- `signs::Union{Vector{Bool}, Nothing}`, storage for random sign flips.
+- `scaling::Float64`, storage for the scaling of the sketches.
+
+Calling `LinSysBlockColSRHT()` defaults to setting `blockSize` to 2.
+
+Nir Ailon and Bernard Chazelle. 2006. Approximate nearest neighbors and the fast Johnson-Lindenstrauss transform. In Proceedings of the thirty-eighth annual ACM symposium on Theory of Computing (STOC '06). Association for Computing Machinery, New York, NY, USA, 557–563. https://doi.org/10.1145/1132516.1132597
+"""
+mutable struct LinSysBlockColSRHT <: LinSysBlkColSampler
+    blockSize::Int64
+    paddedSize::Int64
+    block::Union{Vector{Int64}, Nothing}
+    Ap::Union{AbstractMatrix, Nothing}
+    signs::Union{Vector{Bool}, Nothing}
+    scaling::Float64
+end
+
+LinSysBlockColSRHT(blockSize) = LinSysBlockColSRHT(
+                                                   blockSize, 
+                                                   0, 
+                                                   nothing, 
+                                                   nothing,
+                                                   nothing,
+                                                   0.0
+                                                  )
+LinSysBlockColSRHT() = LinSysBlockColSRHT(2, 0, nothing, nothing, nothing, 0.0)
+
+# Common sample interface for linear systems
+function sample(
+    type::LinSysBlockColSRHT,
+    A::AbstractArray,
+    b::AbstractVector,
+    x::AbstractVector,
+    iter::Int64
+)
+    if iter == 1
+        m, n = size(A)
+        # If matrix is not a power of 2 then pad the rows
+        if rem(log(2, n), 1) != 0
+            type.paddedSize = Int64(2^(div(log(2, n), 1) + 1)) 
+            # Find nearest power 2 and allocate
+            type.Ap = zeros(m, type.paddedSize)
+            # Pad matrix and constant vector
+            type.Ap[:, 1:n] .= A
+        else
+            type.paddedSize = n
+            type.Ap = A
+        end
+        # Compute scaling and sign flips
+        type.scaling = sqrt(type.blockSize / type.paddedSize)
+        type.signs = bitrand(type.paddedSize)
+        for i = 1:m
+            @views fwht!(type.Ap[i, :], signs = type.signs, scaling = type.scaling)
+        end
+
+        type.block = zeros(Int64, type.blockSize) 
+    end
+
+    type.block .= randperm(type.paddedSize)[1:type.blockSize] 
+    AS = type.Ap[:, type.block]
+    # Residual of the linear system
+    res = A * x - b
+    grad = AS'res
+    H = hadamard(type.paddedSize)
+    sgn = [type.signs[i] ? 1 : -1 for i in 1:type.paddedSize]
+    return [sgn .* type.scaling, type.block], AS, res, grad
+end

src/linear_samplers/block_row_FJLT.jl

@@ -0,0 +1,85 @@
+"""
+    LinSysBlockRowFJLT <: LinSysBlkRowSampler
+
+A mutable structure with fields to handle FJLT row sketching. For this procedure,
+the hadamard transform and random sign swaps are applied once, then that matrix is repeatably
+sampled.
+
+# Fields
+- `blockSize::Int64`, the size of the sketching dimension
+- `sparsity::Float64`, the sparsity of the sampling matrix
+- `paddedSize::Int64`, the size of the matrix when padded
+- `Sketch::Union{SparseMatrixCSC, Nothing}`, storage for sparse sketching matrix 
+- `Ap::Union{AbstractMatrix, Nothing}`, storage for padded matrix
+- `bp::Union{AbstractMatrix, Nothing}`, storage for padded vector
+- `signs::Union{Vector{Bool}, Nothing}`, storage for random sign flips.
+- `scaling::Float64`, storage for the scaling of the sketches.
+
+Calling `LinSysBlockRowFJLT()` defaults to setting `sparsity` to .3 and the blocksize to 2.
+
+Nir Ailon and Bernard Chazelle. 2006. Approximate nearest neighbors and the fast Johnson-Lindenstrauss transform. In Proceedings of the thirty-eighth annual ACM symposium on Theory of Computing (STOC '06). Association for Computing Machinery, New York, NY, USA, 557–563. https://doi.org/10.1145/1132516.1132597
+"""
+mutable struct LinSysBlockRowFJLT <: LinSysBlkRowSampler
+    blockSize::Int64
+    sparsity::Float64 
+    paddedSize::Int64
+    Sketch::Union{SparseMatrixCSC, Nothing}
+    Ap::Union{AbstractMatrix, Nothing}
+    bp::Union{AbstractVector, Nothing}
+    signs::Union{Vector{Bool}, Nothing}
+    scaling::Float64
+end
+
+LinSysBlockRowFJLT(;blocksize = 2, sparsity = .3) = LinSysBlockRowFJLT(
+                                                   blocksize,
+                                                   sparsity, 
+                                                   0, 
+                                                   nothing, 
+                                                   nothing,
+                                                   nothing,
+                                                   nothing,
+                                                   0.0
+                                                  )
+
+# Common sample interface for linear systems
+function sample(
+    type::LinSysBlockRowFJLT,
+    A::AbstractArray,
+    b::AbstractVector,
+    x::AbstractVector,
+    iter::Int64
+)
+    if iter == 1
+        m, n = size(A)
+        # If matrix is not a power of 2 then pad the rows
+        if rem(log(2, m), 1) != 0
+            type.paddedSize = Int64(2^(div(log(2, m), 1) + 1)) 
+            # Find nearest power 2 and allocate
+            type.Ap = zeros(type.paddedSize, n)
+            type.bp = zeros(type.paddedSize)
+            # Pad matrix and constant vector
+            @views type.Ap[1:m, :] .= A
+            @views type.bp[1:m] .= b
+        else
+            type.paddedSize = m
+            type.Ap = A
+            type.bp = b
+        end
+        # Compute scaling and sign flips
+        type.scaling = sqrt(type.blockSize / (type.paddedSize * type.sparsity))
+        type.signs = bitrand(type.paddedSize)
+        # Apply FWHT to padded matrix and vector
+        fwht!(type.bp, signs = type.signs, scaling = type.scaling)
+        for i = 1:n
+            @views fwht!(type.Ap[:, i], signs = type.signs, scaling = type.scaling)
+        end
+
+    end
+
+    type.Sketch = sprandn(type.blockSize, type.paddedSize, type.sparsity) 
+    SA = type.Sketch * type.Ap
+    Sb = type.Sketch * type.bp
+    # Residual of the linear system
+    res = SA * x - Sb
+    return type.Sketch, SA, res
+end