From 8e18d8bafc7e14b25548d34c37942744a428da12 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Mon, 7 Jul 2025 21:00:44 +0100 Subject: [PATCH 1/8] initialized IHS --- src/Solvers/ihs.jl | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 src/Solvers/ihs.jl diff --git a/src/Solvers/ihs.jl b/src/Solvers/ihs.jl new file mode 100644 index 00000000..2845854d --- /dev/null +++ b/src/Solvers/ihs.jl @@ -0,0 +1,45 @@ +mutable struct IHSRecipe{M<:AbstractArray, MV<:SubArray, V<:AbstractVector} + logger::LoggerRecipe + compressor::CompressorRecipe + sub_solver::SubSolverRecipe + compressed_mat::M + mat_view::MV + res::V + grad::V + update_vec::V + solution_vec::V +end + +function complete_solver() + +end + +function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::AbstractVector) + reset_logger!(solver.logger) + solver.solution_vec = x + err = 0.0 + for i in 1:solver.log.max_it + mul!(solver.res, A, solver.solution_vec, -1.0, 1.0) + mul!(solver.grad, A', solver.res) + err = compute_error(solver.error, solver, A, b) + update_logger!(solver.logger, err, i) + if solver.log.converged + return solver.solution_vec, solver.log + end + + # generate a new compressor + update_compressor!(solver.compressor, x, A, b) + # Based on the size of the compressor update views of the matrix + rows_s, cols_s = size(solver.compressor) + solver.mat_view = view(solver.compressed_mat, 1:rows_s, :) + # Compress the matrix + mul!(solver.mat_view, solver.compressor, A) + # Update the subsolver + update_sub_solver!(solver.sub_solver, solver.mat_view) + ldiv!(solver.update_vec, solver.sub_solver, solver.grad) + axpy!(1.0, solver.update_vec, 1.0, solver.solution_vec) + end + + return solver.solution_vec, solver + +end From 082eb34e98bd596816c357bf41ed40e163965301 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Tue, 8 Jul 2025 11:46:09 +0100 Subject: [PATCH 2/8] finished complete compressor function --- src/Solvers/ihs.jl | 100 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 88 insertions(+), 12 deletions(-) diff --git a/src/Solvers/ihs.jl b/src/Solvers/ihs.jl index 2845854d..299dd640 100644 --- a/src/Solvers/ihs.jl +++ b/src/Solvers/ihs.jl @@ -1,17 +1,89 @@ -mutable struct IHSRecipe{M<:AbstractArray, MV<:SubArray, V<:AbstractVector} - logger::LoggerRecipe - compressor::CompressorRecipe - sub_solver::SubSolverRecipe +mutable struct IHSRecipe{ + Type<:Number, + LR<:LoggerRecipe, + CR<:CompressorRecipe, + ER<:ErrorRecipe, + M<:AbstractArray, + MV<:SubArray, + V<:AbstractVector +} + logger::LR + compressor::CR + error::ER compressed_mat::M mat_view::MV - res::V - grad::V + residual_vec::V + gradient_vec::V update_vec::V solution_vec::V + R::UpperTriangular{T, M} end -function complete_solver() +function complete_solver( + ingredients::IHS, + x::AbstractVector, + A::AbstractMatrix, + b::AbstractVector +) + compressor = complete_compressor(ingredients.compressor, x, A, b) + logger = complete_logger(ingredients.log) + error = complete_error(ingredients.error) + sample_size::Int64 = compressor.n_rows + rows_a, cols_a = size(A) + # Check that required fields are in the types + if !isdefined(error, :residual) + throw( + ArgumentError( + "ErrorRecipe $(typeof(error)) does not contain the \ + field 'residual' and is not valid for an IHS solver." + ) + ) + end + + if !isdefined(logger, :converged) + throw( + ArgumentError( + "LoggerRecipe $(typeof(logger)) does not contain \ + the field 'converged' and is not valid for an IHS solver." + ) + ) + end + + # Check that the sketch size is larger than the column dimension and return a warning + # otherwise + if cols_a > sample_size + @warn "Compression dimension not larger than column dimension this will lead to \ + singular QR decompositions, which cannot be inverted" + end + + compressed_mat = zeros(ingredients.type, sample_size, cols_a) + res = zeros(ingredients.type, rows_a) + grad = zeros(ingredients.type, cols_a) + update_vec = zeros(ingredients.type, cols_a) + solution_vec = x + mat_view = view(compressed_mat, 1:sample_size, :) + R = UpperTriangular(mat_view[1:cols_a]) + return IHSRecipe{ + ingredients.type, + typeof(logger), + typeof(compressor), + typeof(error), + typeof(compressed_mat), + typeof(mat_view), + typeof(update_vec) + }( + logger, + compressor, + error, + compressed_mat, + mat_view, + res, + grad, + update_vec, + solution_vec, + R + ) end function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::AbstractVector) @@ -19,8 +91,8 @@ function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::Abs solver.solution_vec = x err = 0.0 for i in 1:solver.log.max_it - mul!(solver.res, A, solver.solution_vec, -1.0, 1.0) - mul!(solver.grad, A', solver.res) + mul!(solver.resdual_vec, A, solver.solution_vec, -1.0, 1.0) + mul!(solver.gradient_vec, A', solver.res) err = compute_error(solver.error, solver, A, b) update_logger!(solver.logger, err, i) if solver.log.converged @@ -35,9 +107,13 @@ function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::Abs # Compress the matrix mul!(solver.mat_view, solver.compressor, A) # Update the subsolver - update_sub_solver!(solver.sub_solver, solver.mat_view) - ldiv!(solver.update_vec, solver.sub_solver, solver.grad) - axpy!(1.0, solver.update_vec, 1.0, solver.solution_vec) + solver.R = UpperTriangular(qr!(solver.ma_view).R) + # Compute first R' solver + ldiv!(solver.update_vec, solver.R', solver.gradient_vec) + # Compute second R Solve + ldiv!(solver.gradient_vec, solver.R, solver.update_vec) + # update the solution + axpy!(1.0, solver.gradient_vec, 1.0, solver.solution_vec) end return solver.solution_vec, solver From 7be38a3e0ceaffe2404aea81f00c5cf6a2f21842 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Tue, 8 Jul 2025 20:58:20 +0100 Subject: [PATCH 3/8] added documentation to IHS --- docs/src/refs.bib | 62 ++++++++++++++++++++++++++++++ src/Solvers/ihs.jl | 95 +++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 156 insertions(+), 1 deletion(-) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index b7297ec8..a873f869 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -9,7 +9,9 @@ @article{ailon2009fast pages = {302--322}, issn = {0097-5397, 1095-7111}, doi = {10.1137/060673096}, + urldate = {2025-04-26}, langid = {english}, + file = {/Users/nathanielpritchard/Zotero/storage/WKB6MTMQ/Ailon and Chazelle - 2009 - The Fast Johnson–Lindenstrauss Transform and Approximate Nearest Neighbors.pdf} } @article{halko2011finding, @@ -24,6 +26,7 @@ @article{halko2011finding pages = {217--288}, issn = {0036-1445, 1095-7200}, doi = {10.1137/090771806}, + urldate = {2025-01-24}, langid = {english} } @@ -34,10 +37,41 @@ @misc{martinsson2020randomized year = {2020}, publisher = {arXiv}, doi = {10.48550/ARXIV.2002.01387}, + urldate = {2025-02-10}, + abstract = {This survey describes probabilistic algorithms for linear algebra computations, such as factorizing matrices and solving linear systems. It focuses on techniques that have a proven track record for real-world problem instances. The paper treats both the theoretical foundations of the subject and the practical computational issues. Topics covered include norm estimation; matrix approximation by sampling; structured and unstructured random embeddings; linear regression problems; low-rank approximation; subspace iteration and Krylov methods; error estimation and adaptivity; interpolatory and CUR factorizations; Nystr{\"o}m approximation of positive-semidefinite matrices; single view ("streaming") algorithms; full rank-revealing factorizations; solvers for linear systems; and approximation of kernel matrices that arise in machine learning and in scientific computing.}, copyright = {arXiv.org perpetual, non-exclusive license}, keywords = {FOS: Mathematics,Numerical Analysis (math.NA)} } +@article{needell2014paved, + title = {Paved with Good Intentions: {{Analysis}} of a Randomized Block {{Kaczmarz}} Method}, + shorttitle = {Paved with Good Intentions}, + author = {Needell, Deanna and Tropp, Joel A.}, + year = {2014}, + month = jan, + journal = {Linear Algebra and its Applications}, + volume = {441}, + pages = {199--221}, + issn = {00243795}, + doi = {10.1016/j.laa.2012.12.022}, + urldate = {2025-06-06}, + langid = {english}, + file = {/Users/nathanielpritchard/Zotero/storage/FPYFBHCI/Needell and Tropp - 2014 - Paved with good intentions Analysis of a randomized block Kaczmarz method.pdf} +} + +@misc{pilanci2014iterative, + title = {Iterative {{Hessian}} Sketch: {{Fast}} and Accurate Solution Approximation for Constrained Least-Squares}, + shorttitle = {Iterative {{Hessian}} Sketch}, + author = {Pilanci, Mert and Wainwright, Martin J.}, + year = {2014}, + publisher = {arXiv}, + doi = {10.48550/ARXIV.1411.0347}, + urldate = {2025-07-08}, + abstract = {We study randomized sketching methods for approximately solving least-squares problem with a general convex constraint. The quality of a least-squares approximation can be assessed in different ways: either in terms of the value of the quadratic objective function (cost approximation), or in terms of some distance measure between the approximate minimizer and the true minimizer (solution approximation). Focusing on the latter criterion, our first main result provides a general lower bound on any randomized method that sketches both the data matrix and vector in a least-squares problem; as a surprising consequence, the most widely used least-squares sketch is sub-optimal for solution approximation. We then present a new method known as the iterative Hessian sketch, and show that it can be used to obtain approximations to the original least-squares problem using a projection dimension proportional to the statistical complexity of the least-squares minimizer, and a logarithmic number of iterations. We illustrate our general theory with simulations for both unconstrained and constrained versions of least-squares, including \${\textbackslash}ell\_1\$-regularization and nuclear norm constraints. We also numerically demonstrate the practicality of our approach in a real face expression classification experiment.}, + copyright = {arXiv.org perpetual, non-exclusive license}, + keywords = {FOS: Computer and information sciences,FOS: Mathematics,Information Theory (cs.IT),Machine Learning (cs.LG),Machine Learning (stat.ML),Optimization and Control (math.OC)} +} + @article{pritchard2023practical, title = {Towards {{Practical Large-Scale Randomized Iterative Least Squares Solvers}} through {{Uncertainty Quantification}}}, author = {Pritchard, Nathaniel and Patel, Vivak}, @@ -49,6 +83,7 @@ @article{pritchard2023practical pages = {996--1024}, issn = {2166-2525}, doi = {10.1137/22M1515057}, + urldate = {2025-01-24}, langid = {english} } @@ -63,4 +98,31 @@ @article{pritchard2024solving pages = {085003}, issn = {0266-5611, 1361-6420}, doi = {10.1088/1361-6420/ad5583}, + urldate = {2025-01-24}, + abstract = {Abstract In large-scale applications including medical imaging, collocation differential equation solvers, and estimation with differential privacy, the underlying linear inverse problem can be reformulated as a streaming problem. In theory, the streaming problem can be effectively solved using memory-efficient, exponentially-converging streaming solvers. In special cases when the underlying linear inverse problem is finite-dimensional, streaming solvers can periodically evaluate the residual norm at a substantial computational cost. When the underlying system is infinite dimensional, streaming solver can only access noisy estimates of the residual. While such noisy estimates are computationally efficient, they are useful only when their accuracy is known. In this work, we rigorously develop a general family of computationally-practical residual estimators and their uncertainty sets for streaming solvers, and we demonstrate the accuracy of our methods on a number of large-scale linear problems. Thus, we further enable the practical use of streaming solvers for important classes of linear inverse problems.} } + +@article{strohmer2009randomized, + title = {A {{Randomized Kaczmarz Algorithm}} with {{Exponential Convergence}}}, + author = {Strohmer, Thomas and Vershynin, Roman}, + year = {2009}, + month = apr, + journal = {Journal of Fourier Analysis and Applications}, + volume = {15}, + number = {2}, + pages = {262--278}, + issn = {1069-5869, 1531-5851}, + doi = {10.1007/s00041-008-9030-4}, + urldate = {2025-06-06}, + copyright = {http://www.springer.com/tdm}, + langid = {english} +} +@misc{pilanci2014iterative, + title = {Iterative {{Hessian}} Sketch: {{Fast}} and Accurate Solution Approximation for Constrained Least-Squares}, + shorttitle = {Iterative {{Hessian}} Sketch}, + author = {Pilanci, Mert and Wainwright, Martin J.}, + year = {2014}, + publisher = {arXiv}, + doi = {10.48550/ARXIV.1411.0347}, +} + diff --git a/src/Solvers/ihs.jl b/src/Solvers/ihs.jl index 299dd640..6a5177fd 100644 --- a/src/Solvers/ihs.jl +++ b/src/Solvers/ihs.jl @@ -1,3 +1,96 @@ +""" + IHS <: Solver + +An implementation of the Iterative Hessian Sketch solver for solving over determined +least squares problems (@cite)[pilanci2014iterative]. + +# Mathematical Description +Let ``A`` be an ``m \\times n`` matrix and consider the least square problem ``\\min_x +\\|Ax - b \\|_2^2. Iterative Hessian Sketch + +# Fields +- `alpha::Float64`, a step size parameter. +- `compressor::Compressor`, a technique for forming the compressed linear system. +- `log::Logger`, a technique for logging the progress of the solver. +- `error::SolverError', a method for estimating the progress of the solver. + +# Constructor + function IHS(; + compressor::Compressor = SparseSign(cardinality = Left()), + log::Logger = BasicLogger(), + error::SolverError = FullResidual(), + alpha::Float64 = 1.0 + ) +## Keywords +- `compressor::Compressor`, a technique for forming the compressed linear system. +- `log::Logger`, a technique for logging the progress of the solver. +- `error::SolverError', a method for estimating the progress of the solver. +- `alpha::Float64`, a step size parameter. + +# Returns +- A `IHS` object. +""" +mutable struct IHS <: Solver + alpha::Float64 + log::Logger + compressor::Compressor + error::Error + function IHS(alpha, compressor, log, error) + if typeof(compressor.cardinality) != Left + @warn "Compressor has cardinality `Right` but IHS compresses from the `Left`." + end + + new(alpha, compressor, log, error, subsolver) + end + +end + +function IHS(; + compressor::Compressor = SparseSign(cardinality = Left()), + log::Logger = BasicLogger(), + error::SolverError = FullResidual(), + alpha::Float64 = 1.0 +) + return IHS( + alpha, + compressor, + log, + error + ) +end + +""" + IHSRecipe{ + Type<:Number, + LR<:LoggerRecipe, + CR<:CompressorRecipe, + ER<:ErrorRecipe, + M<:AbstractArray, + MV<:SubArray, + V<:AbstractVector + } <: SolverRecip + +A mutable structure containing all information relevant to the Iterative Hessian Sketch +solver. It is formed by calling the function `complete_solver` on a `IHS` solver, which +includes all the user controlled parameters, the linear system matrix `A` and constant +vector `b`. + +# Fields +- `compressor::CompressorRecipe`, a technique for compressing the matrix ``A``. +- `logger::LoggerRecipe`, a technique for logging the progress of the solver. +- `error::SolverErrorRecipe`, a technique for estimating the progress of the solver. +- `compresed_mat::AbstractMatrix`, a buffer for storing the compressed matrix. +- `mat_view::SubArray`, a container for storing a view of the compressed matrix buffer. +- `residual_vec::AbstractVector`, a vector that contains the residual of the linear system + ``Ax-b``. +- `gradient_vec::AbstractVector`, a vector that contains the gradient of the least squares + problem, ``A^\\top(b-Ax)``. +- `update_vec::AbstractVector`, a vector that contains the update by solving the IHS + subproblem. +- `solution_vec::AbstractVector`, a vector storing the current IHS solution. +- `R::UpperTriangular`, a containter for storing the upper triangular portion of the R + factor from a QR factorization of `mat_view`. This is used to solve the IHS sub-problem. +""" mutable struct IHSRecipe{ Type<:Number, LR<:LoggerRecipe, @@ -6,7 +99,7 @@ mutable struct IHSRecipe{ M<:AbstractArray, MV<:SubArray, V<:AbstractVector -} +} <: SolverRecipe logger::LR compressor::CR error::ER From 009c6b5b0083cd73198a431625154b77492159ed Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Wed, 9 Jul 2025 15:22:38 +0100 Subject: [PATCH 4/8] addded tests for IHS (not converging) --- docs/src/api/solvers.md | 2 + src/RLinearAlgebra.jl | 4 +- src/Solvers.jl | 1 + src/Solvers/ihs.jl | 51 +++-- test/Solvers/ihs.jl | 444 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 479 insertions(+), 23 deletions(-) create mode 100644 test/Solvers/ihs.jl diff --git a/docs/src/api/solvers.md b/docs/src/api/solvers.md index 05f9c1a1..ca98abc3 100644 --- a/docs/src/api/solvers.md +++ b/docs/src/api/solvers.md @@ -12,7 +12,9 @@ SolverRecipe ## Solver Structures ```@docs +IHS +IHSRecipe ``` ## Exported Functions diff --git a/src/RLinearAlgebra.jl b/src/RLinearAlgebra.jl index cfd4324f..8816e23b 100644 --- a/src/RLinearAlgebra.jl +++ b/src/RLinearAlgebra.jl @@ -1,7 +1,8 @@ module RLinearAlgebra import Base.:* import Base: transpose, adjoint -import LinearAlgebra: Adjoint, axpby!, dot, ldiv!, lmul!, lq!, lq, LQ, mul!, norm, qr! +import LinearAlgebra: Adjoint, axpby!, axpy!, dot, ldiv!, lmul!, lq!, lq, LQ, mul!, norm +import LinearAlgebra: qr!, UpperTriangular import StatsBase: sample!, ProbabilityWeights, wsample! import Random: bitrand, rand!, randn! import SparseArrays: SparseMatrixCSC, sprandn @@ -33,6 +34,7 @@ export Uniform, UniformRecipe # Export Solver types and functions export Solver, SolverRecipe export complete_solver, update_solver!, rsolve! +export IHS, IHSRecipe # Export Logger types and functions export Logger, LoggerRecipe diff --git a/src/Solvers.jl b/src/Solvers.jl index ffe0aaa1..4242c0e4 100644 --- a/src/Solvers.jl +++ b/src/Solvers.jl @@ -150,3 +150,4 @@ include("Solvers/ErrorMethods.jl") ############################# # The Solver Routine Files ############################ +include("Solvers/ihs.jl") diff --git a/src/Solvers/ihs.jl b/src/Solvers/ihs.jl index 6a5177fd..2deac486 100644 --- a/src/Solvers/ihs.jl +++ b/src/Solvers/ihs.jl @@ -34,13 +34,13 @@ mutable struct IHS <: Solver alpha::Float64 log::Logger compressor::Compressor - error::Error - function IHS(alpha, compressor, log, error) + error::SolverError + function IHS(alpha, log, compressor, error) if typeof(compressor.cardinality) != Left @warn "Compressor has cardinality `Right` but IHS compresses from the `Left`." end - new(alpha, compressor, log, error, subsolver) + new(alpha, log, compressor, error) end end @@ -53,8 +53,8 @@ function IHS(; ) return IHS( alpha, - compressor, log, + compressor, error ) end @@ -95,21 +95,22 @@ mutable struct IHSRecipe{ Type<:Number, LR<:LoggerRecipe, CR<:CompressorRecipe, - ER<:ErrorRecipe, + ER<:SolverErrorRecipe, M<:AbstractArray, MV<:SubArray, V<:AbstractVector } <: SolverRecipe - logger::LR + log::LR compressor::CR error::ER + alpha::Float64 compressed_mat::M mat_view::MV residual_vec::V gradient_vec::V update_vec::V solution_vec::V - R::UpperTriangular{T, M} + R::UpperTriangular{Type, M} end function complete_solver( @@ -120,7 +121,7 @@ function complete_solver( ) compressor = complete_compressor(ingredients.compressor, x, A, b) logger = complete_logger(ingredients.log) - error = complete_error(ingredients.error) + error = complete_error(ingredients.error, ingredients, A, b) sample_size::Int64 = compressor.n_rows rows_a, cols_a = size(A) # Check that required fields are in the types @@ -145,20 +146,24 @@ function complete_solver( # Check that the sketch size is larger than the column dimension and return a warning # otherwise if cols_a > sample_size - @warn "Compression dimension not larger than column dimension this will lead to \ - singular QR decompositions, which cannot be inverted" + throw( + ArgumentError( + "Compression dimension not larger than column dimension this will lead to \ + singular QR decompositions, which cannot be inverted" + ) + ) end - compressed_mat = zeros(ingredients.type, sample_size, cols_a) - res = zeros(ingredients.type, rows_a) - grad = zeros(ingredients.type, cols_a) - update_vec = zeros(ingredients.type, cols_a) + compressed_mat = zeros(eltype(A), sample_size, cols_a) + res = zeros(eltype(A), rows_a) + grad = zeros(eltype(A), cols_a) + update_vec = zeros(eltype(A), cols_a) solution_vec = x mat_view = view(compressed_mat, 1:sample_size, :) - R = UpperTriangular(mat_view[1:cols_a]) + R = UpperTriangular(mat_view[1:cols_a, :]) return IHSRecipe{ - ingredients.type, + eltype(compressed_mat), typeof(logger), typeof(compressor), typeof(error), @@ -169,6 +174,7 @@ function complete_solver( logger, compressor, error, + ingredients.alpha, compressed_mat, mat_view, res, @@ -180,14 +186,15 @@ function complete_solver( end function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::AbstractVector) - reset_logger!(solver.logger) + reset_logger!(solver.log) solver.solution_vec = x err = 0.0 + copyto!(solver.residual_vec, b) for i in 1:solver.log.max_it - mul!(solver.resdual_vec, A, solver.solution_vec, -1.0, 1.0) - mul!(solver.gradient_vec, A', solver.res) + mul!(solver.residual_vec, A, solver.solution_vec, -1.0, 1.0) + mul!(solver.gradient_vec, A', solver.residual_vec) err = compute_error(solver.error, solver, A, b) - update_logger!(solver.logger, err, i) + update_logger!(solver.log, err, i) if solver.log.converged return solver.solution_vec, solver.log end @@ -200,13 +207,13 @@ function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::Abs # Compress the matrix mul!(solver.mat_view, solver.compressor, A) # Update the subsolver - solver.R = UpperTriangular(qr!(solver.ma_view).R) + solver.R = UpperTriangular(qr!(solver.mat_view).R) # Compute first R' solver ldiv!(solver.update_vec, solver.R', solver.gradient_vec) # Compute second R Solve ldiv!(solver.gradient_vec, solver.R, solver.update_vec) # update the solution - axpy!(1.0, solver.gradient_vec, 1.0, solver.solution_vec) + axpy!(solver.alpha, solver.gradient_vec, solver.solution_vec) end return solver.solution_vec, solver diff --git a/test/Solvers/ihs.jl b/test/Solvers/ihs.jl new file mode 100644 index 00000000..1b4cbc3e --- /dev/null +++ b/test/Solvers/ihs.jl @@ -0,0 +1,444 @@ +module IHSTest +using Test, RLinearAlgebra, LinearAlgebra +import RLinearAlgebra: complete_compressor +import LinearAlgebra: mul!, norm +import Random: randn! +using ..FieldTest +using ..ApproxTol + +# Define the test structures +########################## +# Compressors +########################## +mutable struct ITestCompressor <: Compressor + cardinality::Cardinality + compression_dim::Int64 +end + +ITestCompressor() = ITestCompressor(Left(), 5) + +mutable struct ITestCompressorRecipe <: CompressorRecipe + cardinality::Cardinality + n_rows::Int64 + n_cols::Int64 + op::AbstractMatrix +end + +function RLinearAlgebra.complete_compressor( + comp::ITestCompressor, + A::AbstractMatrix, + b::AbstractVector +) + n_rows = comp.compression_dim + n_cols = size(A, 1) + # Make a gaussian compressor + op = randn(n_rows, n_cols) ./ sqrt(n_cols) + return ITestCompressorRecipe(comp.cardinality, n_rows, n_cols, op) +end + +function RLinearAlgebra.update_compressor!( + comp::ITestCompressorRecipe, + x::AbstractVector, + A::AbstractMatrix, + b::AbstractVector +) + randn!(comp.op) + comp.op ./= sqrt(comp.n_cols) +end + +# Define a mul function for the test compressor +function RLinearAlgebra.mul!( + C::AbstractArray, + S::Main.IHSTest.ITestCompressorRecipe, + A::AbstractArray, + alpha::Number, + beta::Number +) + mul!(C, S.op, A, alpha, beta) +end + +########################## +# Error Method +########################## +mutable struct ITestError <: RLinearAlgebra.SolverError + g::Real +end + +mutable struct ITestErrorRecipe <: RLinearAlgebra.SolverErrorRecipe + residual::Vector{Number} +end + +ITestError() = ITestError(1.0) + +function RLinearAlgebra.complete_error( + error::ITestError, + solver::IHS, + A::AbstractMatrix, + b::AbstractVector +) + return ITestErrorRecipe(zeros(typeof(error.g), size(A, 1))) +end + +function RLinearAlgebra.compute_error(error::ITestErrorRecipe, solver, A, b) + error.residual = A * solver.solution_vec - b + return norm(error.residual) +end + +############################## +# Residual-less Error Recipe +############################## +mutable struct ITestErrorNoRes <: RLinearAlgebra.SolverError end + +mutable struct ITestErrorRecipeNoRes <: RLinearAlgebra.SolverErrorRecipe end + +function RLinearAlgebra.complete_error( + error::ITestErrorNoRes, + solver::IHS, + A::AbstractMatrix, + b::AbstractVector +) + return ITestErrorRecipeNoRes() +end + +############################ +# Loggers +############################ +mutable struct ITestLog <: Logger + max_it::Int64 + g::Number +end + +ITestLog() = ITestLog(5, 1.0) + +ITestLog(max_it) = ITestLog(max_it, 1.0) + +mutable struct ITestLogRecipe <: LoggerRecipe + max_it::Int64 + hist::Vector{Real} + thresh::Float64 + converged::Bool +end + +function RLinearAlgebra.complete_logger(logger::ITestLog) + return ITestLogRecipe( + logger.max_it, + zeros(typeof(logger.g), logger.max_it), + logger.g, + false + ) +end + +function RLinearAlgebra.update_logger!(logger::ITestLogRecipe, err::Real, i::Int64) + logger.hist[i] = err + logger.converged = err < logger.thresh ? true : false +end + +function RLinearAlgebra.reset_logger!(logger::ITestLogRecipe) + fill!(logger.hist, 0.0) +end + +############################## +# Converged-less Logger +############################## +mutable struct ITestLogNoCov <: Logger end + +mutable struct ITestLogRecipeNoCov <: LoggerRecipe end + +function RLinearAlgebra.complete_logger(logger::ITestLogNoCov) + return ITestLogRecipeNoCov() +end + + +@testset "IHS" begin + n_rows = 40 + n_cols = 2 + A = rand(n_rows, n_cols) + xsol = rand(n_cols) + b = A * xsol + + @testset "IHS" begin + @test supertype(IHS) == Solver + + # test fieldnames and types + @test fieldnames(IHS) == (:alpha, :log, :compressor, :error) + @test fieldtypes(IHS) == (Float64, Logger, Compressor, SolverError) + + # test default constructor + + let solver = IHS() + @test solver.alpha == 1.0 + @test typeof(solver.compressor) == SparseSign + @test typeof(solver.compressor.cardinality) == Left + @test typeof(solver.log) == BasicLogger + @test typeof(solver.error) == FullResidual + end + + # test constructor + let solver = IHS( + alpha = 2.0, + compressor = ITestCompressor(), + log = ITestLog(), + error = ITestError(), + ) + + @test solver.alpha == 2.0 + @test typeof(solver.compressor) == ITestCompressor + @test typeof(solver.log) == ITestLog + @test typeof(solver.error) == ITestError + end + + # Test that error gets returned with right compressor + @test_logs (:warn, + "Compressor has cardinality `Right` but IHS compresses from the `Left`." + ) IHS( + alpha = 2.0, + compressor = ITestCompressor(Right(), 5), + log = ITestLog(), + error = ITestError(), + ) + end + + @testset "IHSRecipe" begin + @test supertype(IHSRecipe) == SolverRecipe + + # test fieldnames and types + @test fieldnames(IHSRecipe) == ( + :log, + :compressor, + :error, + :alpha, + :compressed_mat, + :mat_view, + :residual_vec, + :gradient_vec, + :update_vec, + :solution_vec, + :R + ) + @test fieldtypes(IHSRecipe) == ( + LoggerRecipe, + CompressorRecipe, + SolverErrorRecipe, + Float64, + AbstractArray, + SubArray, + AbstractVector, + AbstractVector, + AbstractVector, + AbstractVector, + UpperTriangular{Type, M} where {Type<:Number, M<:AbstractArray} + ) + end + + @testset "IHS: Complete Solver" begin + # test error method with no residual error + let A = A, + xsol = xsol, + b = b, + comp_dim = 2 * size(A, 2), + alpha = 1.0, + n_rows = size(A, 1), + n_cols = size(A, 2), + x = zeros(n_cols) + + comp = ITestCompressor(Left(), comp_dim) + log = ITestLog() + err = ITestErrorNoRes() + solver = IHS( + log = log, + compressor = comp, + error = err, + alpha = alpha + ) + + @test_throws ArgumentError( + "ErrorRecipe $(typeof(ITestErrorRecipeNoRes())) does not contain the \ + field 'residual' and is not valid for an IHS solver." + ) complete_solver(solver, x, A, b) + end + + # test logger method with no converged field + let A = A, + xsol = xsol, + b = b, + comp_dim = 2 * size(A, 2), + alpha = 1.0, + n_rows = size(A, 1), + n_cols = size(A, 2), + x = zeros(n_cols) + + comp = ITestCompressor(Left(), comp_dim) + log = ITestLogNoCov() + err = ITestError() + solver = IHS( + compressor = comp, + log = log, + error = err, + alpha = alpha + ) + + @test_throws ArgumentError( + "LoggerRecipe $(typeof(ITestLogRecipeNoCov())) does not contain \ + the field 'converged' and is not valid for an IHS solver." + ) complete_solver(solver, x, A, b) + end + + # Test the error message about not large enough compression_dim + let A = A, + xsol = xsol, + b = b, + comp_dim = 1, + alpha = 1.0, + n_rows = size(A, 1), + n_cols = size(A, 2), + x = zeros(n_cols) + + comp = ITestCompressor(Left(), comp_dim) + log = ITestLog() + err = ITestError() + solver = IHS( + compressor = comp, + log = log, + error = err, + alpha = alpha + ) + + @test_throws ArgumentError( + "Compression dimension not larger than column dimension this will lead to \ + singular QR decompositions, which cannot be inverted" + ) complete_solver(solver, x, A, b) + end + + let A = A, + xsol = xsol, + b = b, + comp_dim = 8 * size(A, 2), + alpha = 1.0, + n_rows = size(A, 1), + n_cols = size(A, 2), + x = zeros(n_cols) + + comp = ITestCompressor(Left(), comp_dim) + log = ITestLog() + err = ITestError() + solver = IHS( + compressor = comp, + log = log, + error = err, + alpha = alpha + ) + + solver_rec = complete_solver(solver, x, A, b) + + # test types of the contents of the solver + @test typeof(solver_rec) == IHSRecipe{ + Float64, + Main.IHSTest.ITestLogRecipe, + Main.IHSTest.ITestCompressorRecipe, + Main.IHSTest.ITestErrorRecipe, + Matrix{Float64}, + SubArray{ + Float64, + 2, + Matrix{Float64}, + Tuple{UnitRange{Int64}, Base.Slice{Base.OneTo{Int64}}}, + false + }, + Vector{Float64}, + } + @test typeof(solver_rec.compressor) == ITestCompressorRecipe + @test typeof(solver_rec.log) == ITestLogRecipe + @test typeof(solver_rec.error) == ITestErrorRecipe + @test typeof(solver_rec.alpha) == Float64 + @test typeof(solver_rec.compressed_mat) == Matrix{Float64} + @test typeof(solver_rec.solution_vec) == Vector{Float64} + @test typeof(solver_rec.update_vec) == Vector{Float64} + @test typeof(solver_rec.gradient_vec) == Vector{Float64} + @test typeof(solver_rec.residual_vec) == Vector{Float64} + @test typeof(solver_rec.mat_view) <: SubArray + @test typeof(solver_rec.R) <: UpperTriangular + # Test sizes of vectors and matrices + @test size(solver_rec.compressor) == (comp_dim, n_rows) + @test size(solver_rec.compressed_mat) == (comp_dim, n_cols) + @test size(solver_rec.update_vec) == (n_cols,) + @test size(solver_rec.solution_vec) == (n_cols,) + @test size(solver_rec.gradient_vec) == (n_cols,) + @test size(solver_rec.residual_vec) == (n_rows,) + + # test values of entries + solver_rec.alpha == alpha + solver_rec.solution_vec == x + solver_rec.update_vec == zeros(n_cols) + solver_rec.gradient_vec == zeros(n_cols) + solver_rec.residual_vec == zeros(n_rows) + end + + end + + @testset "IHS: rsolve!" begin + for type in [Float16, Float32, Float64, ComplexF32, ComplexF64] + # test maxit stop + let A = rand(type, n_rows, n_cols), + xsol = ones(type, n_cols), + b = A * xsol, + comp_dim = 10 * size(A, 2), + alpha = 1.0, + n_rows = size(A, 1), + n_cols = size(A, 2), + x = zeros(type, n_cols) + x_st = deepcopy(x) + + comp = ITestCompressor(Left(), comp_dim) + #check 10 iterations + log = ITestLog(10, 0.0) + err = ITestError() + solver = IHS( + compressor = comp, + log = log, + error = err, + alpha = alpha + ) + + solver_rec = complete_solver(solver, x, A, b) + + result = rsolve!(solver_rec, x, A, b) + #test that the error decreases + @test norm(A * x_st - b) > norm(A * x - b) + end + + # using threshold stop + let A = rand(type, n_rows, n_cols), + xsol = ones(type, n_cols), + b = A * xsol, + comp_dim = 10 * size(A, 2), + alpha = 1.0, + n_rows = size(A, 1), + n_cols = size(A, 2), + x = zeros(type, n_cols) + x_st = deepcopy(x) + + comp = ITestCompressor(Left(), comp_dim) + #check 20 iterations + log = ITestLog(20, 0.05) + err = ITestError() + solver = IHS( + compressor = comp, + log = log, + error = err, + alpha = alpha + ) + + solver_rec = complete_solver(solver, x, A, b) + + result = rsolve!(solver_rec, x, A, b) + #test that the error decreases + @test norm(A * x_st - b) > norm(A * x - b) + end + + end + + end + +end + +end From db48731651b212250908a62d1359a953dea70695 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Fri, 11 Jul 2025 17:51:42 +0100 Subject: [PATCH 5/8] corrected test issues and cleaned documentation --- src/Solvers/ihs.jl | 55 ++++++++++++++++++++++++++++----------------- test/Solvers/ihs.jl | 31 +++++++++++++------------ 2 files changed, 51 insertions(+), 35 deletions(-) diff --git a/src/Solvers/ihs.jl b/src/Solvers/ihs.jl index 2deac486..1616edb2 100644 --- a/src/Solvers/ihs.jl +++ b/src/Solvers/ihs.jl @@ -5,14 +5,21 @@ An implementation of the Iterative Hessian Sketch solver for solving over determ least squares problems (@cite)[pilanci2014iterative]. # Mathematical Description -Let ``A`` be an ``m \\times n`` matrix and consider the least square problem ``\\min_x -\\|Ax - b \\|_2^2. Iterative Hessian Sketch +Let ``A \\in \\mathbb{R}^{m \\times n}`` and consider the least square problem ``\\min_x +\\|Ax - b \\|_2^2``. If we let ``S \\in \\mathbb{R}^{s \\times m}`` be a compression matrix, then +Iterative Hessian Sketch iteratively finds a solution to this problem +by repeatedly updating ``x_{k+1} = x_k + \\alpha u_k``where ``u_k`` is the solution to the +convex optimization problem, +``u_k = \\min_u \\{\\|S_k Au\\|_2^2 - \\langle A, b - Ax_k \\rangle \\}.`` This method +has been to shown to converge geometrically at a rate ``\\rho \\in (0, 1/2]``, typically the +required compression dimension needs to be 4-8 times the size of n for the algorithm to +perform successfully. # Fields - `alpha::Float64`, a step size parameter. - `compressor::Compressor`, a technique for forming the compressed linear system. - `log::Logger`, a technique for logging the progress of the solver. -- `error::SolverError', a method for estimating the progress of the solver. +- `error::SolverError`, a method for estimating the progress of the solver. # Constructor function IHS(; @@ -37,7 +44,7 @@ mutable struct IHS <: Solver error::SolverError function IHS(alpha, log, compressor, error) if typeof(compressor.cardinality) != Left - @warn "Compressor has cardinality `Right` but IHS compresses from the `Left`." + @warn "Compressor has cardinality `Right` but IHS compresses from the `Left`." end new(alpha, log, compressor, error) @@ -72,23 +79,22 @@ end A mutable structure containing all information relevant to the Iterative Hessian Sketch solver. It is formed by calling the function `complete_solver` on a `IHS` solver, which -includes all the user controlled parameters, the linear system matrix `A` and constant +includes all the user controlled parameters, the linear system `A`, and the constant vector `b`. # Fields - `compressor::CompressorRecipe`, a technique for compressing the matrix ``A``. - `logger::LoggerRecipe`, a technique for logging the progress of the solver. - `error::SolverErrorRecipe`, a technique for estimating the progress of the solver. -- `compresed_mat::AbstractMatrix`, a buffer for storing the compressed matrix. +- `compressed_mat::AbstractMatrix`, a buffer for storing the compressed matrix. - `mat_view::SubArray`, a container for storing a view of the compressed matrix buffer. - `residual_vec::AbstractVector`, a vector that contains the residual of the linear system ``Ax-b``. - `gradient_vec::AbstractVector`, a vector that contains the gradient of the least squares problem, ``A^\\top(b-Ax)``. -- `update_vec::AbstractVector`, a vector that contains the update by solving the IHS - subproblem. +- `buffer_vec::AbstractVector`, a buffer vector for storing intermediate linear system solves. - `solution_vec::AbstractVector`, a vector storing the current IHS solution. -- `R::UpperTriangular`, a containter for storing the upper triangular portion of the R +- `R::UpperTriangular`, a container for storing the upper triangular portion of the R factor from a QR factorization of `mat_view`. This is used to solve the IHS sub-problem. """ mutable struct IHSRecipe{ @@ -108,7 +114,7 @@ mutable struct IHSRecipe{ mat_view::MV residual_vec::V gradient_vec::V - update_vec::V + buffer_vec::V solution_vec::V R::UpperTriangular{Type, M} end @@ -149,7 +155,7 @@ function complete_solver( throw( ArgumentError( "Compression dimension not larger than column dimension this will lead to \ - singular QR decompositions, which cannot be inverted" + singular QR decompositions, which cannot be inverted." ) ) end @@ -157,7 +163,7 @@ function complete_solver( compressed_mat = zeros(eltype(A), sample_size, cols_a) res = zeros(eltype(A), rows_a) grad = zeros(eltype(A), cols_a) - update_vec = zeros(eltype(A), cols_a) + buffer_vec = zeros(eltype(A), cols_a) solution_vec = x mat_view = view(compressed_mat, 1:sample_size, :) R = UpperTriangular(mat_view[1:cols_a, :]) @@ -169,7 +175,7 @@ function complete_solver( typeof(error), typeof(compressed_mat), typeof(mat_view), - typeof(update_vec) + typeof(buffer_vec) }( logger, compressor, @@ -179,7 +185,7 @@ function complete_solver( mat_view, res, grad, - update_vec, + buffer_vec, solution_vec, R ) @@ -190,13 +196,15 @@ function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::Abs solver.solution_vec = x err = 0.0 copyto!(solver.residual_vec, b) + # compute the initial residual r = b - Ax + mul!(solver.residual_vec, A, solver.solution_vec, -1.0, 1.0) for i in 1:solver.log.max_it - mul!(solver.residual_vec, A, solver.solution_vec, -1.0, 1.0) + # compute the gradient A'r mul!(solver.gradient_vec, A', solver.residual_vec) err = compute_error(solver.error, solver, A, b) update_logger!(solver.log, err, i) if solver.log.converged - return solver.solution_vec, solver.log + return nothing end # generate a new compressor @@ -207,15 +215,20 @@ function rsolve!(solver::IHSRecipe, x::AbstractVector, A::AbstractMatrix, b::Abs # Compress the matrix mul!(solver.mat_view, solver.compressor, A) # Update the subsolver + # This is the only piece of allocating code solver.R = UpperTriangular(qr!(solver.mat_view).R) - # Compute first R' solver - ldiv!(solver.update_vec, solver.R', solver.gradient_vec) - # Compute second R Solve - ldiv!(solver.gradient_vec, solver.R, solver.update_vec) + # Compute first R' solver R'R x = g + ldiv!(solver.buffer_vec, solver.R', solver.gradient_vec) + # Compute second R Solve Rx = (R')^(-1)g will be stored in gradient_vec + ldiv!(solver.gradient_vec, solver.R, solver.buffer_vec) # update the solution + # solver.solution_vec = solver.solution_vec + alpha * solver.gradient_vec axpy!(solver.alpha, solver.gradient_vec, solver.solution_vec) + # compute the fast update of r = r - A * gradient_vec + # note: in this case gradient vec stores the update + mul!(solver.residual_vec, A, solver.gradient_vec, -1.0, 1.0) end - return solver.solution_vec, solver + return nothing end diff --git a/test/Solvers/ihs.jl b/test/Solvers/ihs.jl index 1b4cbc3e..c806f4c8 100644 --- a/test/Solvers/ihs.jl +++ b/test/Solvers/ihs.jl @@ -150,7 +150,7 @@ end @testset "IHS" begin - n_rows = 40 + n_rows = 20 n_cols = 2 A = rand(n_rows, n_cols) xsol = rand(n_cols) @@ -189,7 +189,7 @@ end # Test that error gets returned with right compressor @test_logs (:warn, - "Compressor has cardinality `Right` but IHS compresses from the `Left`." + "Compressor has cardinality `Right` but IHS compresses from the `Left`." ) IHS( alpha = 2.0, compressor = ITestCompressor(Right(), 5), @@ -211,7 +211,7 @@ end :mat_view, :residual_vec, :gradient_vec, - :update_vec, + :buffer_vec, :solution_vec, :R ) @@ -305,7 +305,7 @@ end @test_throws ArgumentError( "Compression dimension not larger than column dimension this will lead to \ - singular QR decompositions, which cannot be inverted" + singular QR decompositions, which cannot be inverted." ) complete_solver(solver, x, A, b) end @@ -352,7 +352,7 @@ end @test typeof(solver_rec.alpha) == Float64 @test typeof(solver_rec.compressed_mat) == Matrix{Float64} @test typeof(solver_rec.solution_vec) == Vector{Float64} - @test typeof(solver_rec.update_vec) == Vector{Float64} + @test typeof(solver_rec.buffer_vec) == Vector{Float64} @test typeof(solver_rec.gradient_vec) == Vector{Float64} @test typeof(solver_rec.residual_vec) == Vector{Float64} @test typeof(solver_rec.mat_view) <: SubArray @@ -360,7 +360,7 @@ end # Test sizes of vectors and matrices @test size(solver_rec.compressor) == (comp_dim, n_rows) @test size(solver_rec.compressed_mat) == (comp_dim, n_cols) - @test size(solver_rec.update_vec) == (n_cols,) + @test size(solver_rec.buffer_vec) == (n_cols,) @test size(solver_rec.solution_vec) == (n_cols,) @test size(solver_rec.gradient_vec) == (n_cols,) @test size(solver_rec.residual_vec) == (n_rows,) @@ -368,7 +368,7 @@ end # test values of entries solver_rec.alpha == alpha solver_rec.solution_vec == x - solver_rec.update_vec == zeros(n_cols) + solver_rec.buffer_vec == zeros(n_cols) solver_rec.gradient_vec == zeros(n_cols) solver_rec.residual_vec == zeros(n_rows) end @@ -378,10 +378,11 @@ end @testset "IHS: rsolve!" begin for type in [Float16, Float32, Float64, ComplexF32, ComplexF64] # test maxit stop - let A = rand(type, n_rows, n_cols), + let A = Array(qr(rand(type, n_rows, n_cols)).Q), xsol = ones(type, n_cols), b = A * xsol, - comp_dim = 10 * size(A, 2), + # need to choose compression dim to be large enough + comp_dim = 7 * size(A, 2), alpha = 1.0, n_rows = size(A, 1), n_cols = size(A, 2), @@ -402,15 +403,16 @@ end solver_rec = complete_solver(solver, x, A, b) result = rsolve!(solver_rec, x, A, b) - #test that the error decreases + #test that the residual decreases @test norm(A * x_st - b) > norm(A * x - b) end # using threshold stop - let A = rand(type, n_rows, n_cols), + let A = Array(qr(rand(type, n_rows, n_cols)).Q), xsol = ones(type, n_cols), b = A * xsol, - comp_dim = 10 * size(A, 2), + # need to choose compression dim to be large enough + comp_dim = 7 * size(A, 2), alpha = 1.0, n_rows = size(A, 1), n_cols = size(A, 2), @@ -418,8 +420,8 @@ end x_st = deepcopy(x) comp = ITestCompressor(Left(), comp_dim) - #check 20 iterations - log = ITestLog(20, 0.05) + #check 40 iterations we make a 1% improvement + log = ITestLog(40, norm(b) * .99) err = ITestError() solver = IHS( compressor = comp, @@ -433,6 +435,7 @@ end result = rsolve!(solver_rec, x, A, b) #test that the error decreases @test norm(A * x_st - b) > norm(A * x - b) + @test solver_rec.log.converged end end From c48420ff5b73e73b88a21cc06c4c85af989cd1f2 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Fri, 11 Jul 2025 19:56:30 +0100 Subject: [PATCH 6/8] more refs clean ups --- docs/src/refs.bib | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index a873f869..e6be4075 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -9,9 +9,7 @@ @article{ailon2009fast pages = {302--322}, issn = {0097-5397, 1095-7111}, doi = {10.1137/060673096}, - urldate = {2025-04-26}, langid = {english}, - file = {/Users/nathanielpritchard/Zotero/storage/WKB6MTMQ/Ailon and Chazelle - 2009 - The Fast Johnson–Lindenstrauss Transform and Approximate Nearest Neighbors.pdf} } @article{halko2011finding, @@ -26,7 +24,6 @@ @article{halko2011finding pages = {217--288}, issn = {0036-1445, 1095-7200}, doi = {10.1137/090771806}, - urldate = {2025-01-24}, langid = {english} } @@ -37,10 +34,6 @@ @misc{martinsson2020randomized year = {2020}, publisher = {arXiv}, doi = {10.48550/ARXIV.2002.01387}, - urldate = {2025-02-10}, - abstract = {This survey describes probabilistic algorithms for linear algebra computations, such as factorizing matrices and solving linear systems. It focuses on techniques that have a proven track record for real-world problem instances. The paper treats both the theoretical foundations of the subject and the practical computational issues. Topics covered include norm estimation; matrix approximation by sampling; structured and unstructured random embeddings; linear regression problems; low-rank approximation; subspace iteration and Krylov methods; error estimation and adaptivity; interpolatory and CUR factorizations; Nystr{\"o}m approximation of positive-semidefinite matrices; single view ("streaming") algorithms; full rank-revealing factorizations; solvers for linear systems; and approximation of kernel matrices that arise in machine learning and in scientific computing.}, - copyright = {arXiv.org perpetual, non-exclusive license}, - keywords = {FOS: Mathematics,Numerical Analysis (math.NA)} } @article{needell2014paved, @@ -54,9 +47,7 @@ @article{needell2014paved pages = {199--221}, issn = {00243795}, doi = {10.1016/j.laa.2012.12.022}, - urldate = {2025-06-06}, langid = {english}, - file = {/Users/nathanielpritchard/Zotero/storage/FPYFBHCI/Needell and Tropp - 2014 - Paved with good intentions Analysis of a randomized block Kaczmarz method.pdf} } @misc{pilanci2014iterative, @@ -66,9 +57,6 @@ @misc{pilanci2014iterative year = {2014}, publisher = {arXiv}, doi = {10.48550/ARXIV.1411.0347}, - urldate = {2025-07-08}, - abstract = {We study randomized sketching methods for approximately solving least-squares problem with a general convex constraint. The quality of a least-squares approximation can be assessed in different ways: either in terms of the value of the quadratic objective function (cost approximation), or in terms of some distance measure between the approximate minimizer and the true minimizer (solution approximation). Focusing on the latter criterion, our first main result provides a general lower bound on any randomized method that sketches both the data matrix and vector in a least-squares problem; as a surprising consequence, the most widely used least-squares sketch is sub-optimal for solution approximation. We then present a new method known as the iterative Hessian sketch, and show that it can be used to obtain approximations to the original least-squares problem using a projection dimension proportional to the statistical complexity of the least-squares minimizer, and a logarithmic number of iterations. We illustrate our general theory with simulations for both unconstrained and constrained versions of least-squares, including \${\textbackslash}ell\_1\$-regularization and nuclear norm constraints. We also numerically demonstrate the practicality of our approach in a real face expression classification experiment.}, - copyright = {arXiv.org perpetual, non-exclusive license}, keywords = {FOS: Computer and information sciences,FOS: Mathematics,Information Theory (cs.IT),Machine Learning (cs.LG),Machine Learning (stat.ML),Optimization and Control (math.OC)} } @@ -83,7 +71,6 @@ @article{pritchard2023practical pages = {996--1024}, issn = {2166-2525}, doi = {10.1137/22M1515057}, - urldate = {2025-01-24}, langid = {english} } @@ -98,8 +85,6 @@ @article{pritchard2024solving pages = {085003}, issn = {0266-5611, 1361-6420}, doi = {10.1088/1361-6420/ad5583}, - urldate = {2025-01-24}, - abstract = {Abstract In large-scale applications including medical imaging, collocation differential equation solvers, and estimation with differential privacy, the underlying linear inverse problem can be reformulated as a streaming problem. In theory, the streaming problem can be effectively solved using memory-efficient, exponentially-converging streaming solvers. In special cases when the underlying linear inverse problem is finite-dimensional, streaming solvers can periodically evaluate the residual norm at a substantial computational cost. When the underlying system is infinite dimensional, streaming solver can only access noisy estimates of the residual. While such noisy estimates are computationally efficient, they are useful only when their accuracy is known. In this work, we rigorously develop a general family of computationally-practical residual estimators and their uncertainty sets for streaming solvers, and we demonstrate the accuracy of our methods on a number of large-scale linear problems. Thus, we further enable the practical use of streaming solvers for important classes of linear inverse problems.} } @article{strohmer2009randomized, @@ -113,8 +98,6 @@ @article{strohmer2009randomized pages = {262--278}, issn = {1069-5869, 1531-5851}, doi = {10.1007/s00041-008-9030-4}, - urldate = {2025-06-06}, - copyright = {http://www.springer.com/tdm}, langid = {english} } @misc{pilanci2014iterative, From 4bb21ee4adb00bc898813846b63647433bdbcb60 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Fri, 11 Jul 2025 20:31:17 +0100 Subject: [PATCH 7/8] corrected test failure issue --- test/Solvers/ihs.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/Solvers/ihs.jl b/test/Solvers/ihs.jl index c806f4c8..202e2936 100644 --- a/test/Solvers/ihs.jl +++ b/test/Solvers/ihs.jl @@ -150,6 +150,7 @@ end @testset "IHS" begin + Random.seed!(12312) n_rows = 20 n_cols = 2 A = rand(n_rows, n_cols) @@ -382,7 +383,7 @@ end xsol = ones(type, n_cols), b = A * xsol, # need to choose compression dim to be large enough - comp_dim = 7 * size(A, 2), + comp_dim = 9 * size(A, 2), alpha = 1.0, n_rows = size(A, 1), n_cols = size(A, 2), @@ -412,7 +413,7 @@ end xsol = ones(type, n_cols), b = A * xsol, # need to choose compression dim to be large enough - comp_dim = 7 * size(A, 2), + comp_dim = 9 * size(A, 2), alpha = 1.0, n_rows = size(A, 1), n_cols = size(A, 2), From 7cc79b08de6ac344951a28233e3a082988064bf2 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Fri, 11 Jul 2025 20:36:01 +0100 Subject: [PATCH 8/8] added seed to test --- test/Solvers/ihs.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Solvers/ihs.jl b/test/Solvers/ihs.jl index 202e2936..1eee8d20 100644 --- a/test/Solvers/ihs.jl +++ b/test/Solvers/ihs.jl @@ -2,7 +2,7 @@ module IHSTest using Test, RLinearAlgebra, LinearAlgebra import RLinearAlgebra: complete_compressor import LinearAlgebra: mul!, norm -import Random: randn! +import Random: randn!, seed! using ..FieldTest using ..ApproxTol @@ -150,7 +150,7 @@ end @testset "IHS" begin - Random.seed!(12312) + seed!(12312) n_rows = 20 n_cols = 2 A = rand(n_rows, n_cols)