From 44a288bf3e9a18aaee2e96f1f486993c10e003f0 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Fri, 10 May 2024 09:15:20 -0500 Subject: [PATCH 01/14] Corrected Scaling issue so now sketched block expectation is same as true block --- docs/src/man/cls_overview.md | 5 +++++ src/linear_solver_logs/solve_log_ma.jl | 17 ++++++++++++----- test/linear_solver_logs/proc_solve_log_ma.jl | 4 ++-- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/docs/src/man/cls_overview.md b/docs/src/man/cls_overview.md index c8d5a51..d4abd90 100644 --- a/docs/src/man/cls_overview.md +++ b/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. diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 4f54376..75f94d3 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -38,6 +38,7 @@ A mutable structure that stores information about the sub-Exponential distributi - `omega::Union{Float64, Nothing}`, The exponential distrbution parameter, if not given is determined for sampling methods. - `eta::Float64`, A parameter for adjusting the conservativeness of the distribution, higher value means a less conservative estimate. By default, this is set to one. +- `scaling::Float64`, constant multiplied by norm to ensure expectation of block norms is the same as the full norm. """ mutable struct DistInfo sampler::Union{DataType, Nothing} @@ -46,6 +47,7 @@ mutable struct DistInfo sigma2::Union{Float64, Nothing} omega::Union{Float64, Nothing} eta::Float64 + scaling::Float64 end """ @@ -114,7 +116,7 @@ LSLogMA(; -1, false, true_res, - DistInfo(nothing, 0, 0, sigma2, omega, eta) + DistInfo(nothing, 0, 0, sigma2, omega, eta, 0) ) #Function to update the moving average @@ -147,7 +149,7 @@ function log_update!( log.dist_info.sampler = typeof(sampler) # If the constants for the sub-Exponential distribution are not defined then define them - if typeof(log.dist_info.sigma2) <: Nothing + if typeof(log.dist_info.sigma2) <: Nothing || log.dist_info.sigma2 == 0 get_SE_constants!(log, log.dist_info.sampler) end @@ -158,8 +160,10 @@ function log_update!( #Check if we want exact residuals computed if !log.true_res && iter > 0 # Compute the current residual to second power to align with theory - res::Float64 = eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? - log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2 + res::Float64 = log.dist_info.scaling * + (eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? + log.resid_norm(samp[3])^2 : + log.resid_norm(dot(samp[1], x) - samp[2])^2) else res = log.resid_norm(A * x - b)^2 end @@ -305,7 +309,8 @@ A function that returns a default set of sub-Exponential constants for each samp - `sampler::Type{LinSysSampler}`, the type of sampler being used. # Outputs -Performs an inplace update of the sub-Exponential constants for the log. +Performs an inplace update of the sub-Exponential constants for the log. Additionally, updates the scaling constant to ensure expectation of +block norms is equal to true norm. """ function get_SE_constants!(log::LSLogMA, sampler::Type{T}) where T<:LinSysSampler return nothing @@ -320,6 +325,7 @@ for type in (LinSysVecRowDetermCyclic,LinSysVecRowHopRandCyclic, @eval begin function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) log.dist_info.sigma2 = log.dist_info.max_dimension^2 / (4 * log.dist_info.block_dimension^2 * log.dist_info.eta) + log.dist_info.scaling = log.dist_info.max_dimension / log.dist_info.block_dimension end end @@ -332,6 +338,7 @@ for type in (LinSysVecColOneRandCyclic, LinSysVecColDetermCyclic) @eval begin function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) log.dist_info.sigma2 = log.dist_info.max_dimension^2 / (4 * log.dist_info.block_dimension^2 * log.dist_info.eta) + log.dist_info.scaling = log.dist_info.max_dimension / log.dist_info.block_dimension end end diff --git a/test/linear_solver_logs/proc_solve_log_ma.jl b/test/linear_solver_logs/proc_solve_log_ma.jl index 4c4b446..830ff05 100644 --- a/test/linear_solver_logs/proc_solve_log_ma.jl +++ b/test/linear_solver_logs/proc_solve_log_ma.jl @@ -53,8 +53,8 @@ Random.seed!(1010) RLinearAlgebra.log_update!(log, sampler, x + (i+1)*(z-x), samp, i, A, b) end #compute sampled residuals - obs_res = [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^2 for i = 0:10] - obs_res2 = [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^4 for i = 0:10] + obs_res = 2 .* [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^2 for i = 0:10] + obs_res2 = 4 .* [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^4 for i = 0:10] @test length(log.resid_hist) == 11 @test norm(log.resid_hist[2:11] - vcat(obs_res[2], [(obs_res[i] + obs_res[i-1])/2 for i = 3:11])) < 1e2 * eps() From 893a850220817359e86889f58416c059d9cf973e Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Sat, 11 May 2024 22:06:40 -0500 Subject: [PATCH 02/14] Fixed bug with lambda1>1 and added test to prevent issues --- src/linear_solver_logs/solve_log_ma.jl | 33 +++++++++++++-- test/linear_solver_logs/proc_solve_log_ma.jl | 43 +++++++++++++++++--- 2 files changed, 68 insertions(+), 8 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 75f94d3..7ed83a3 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -108,7 +108,7 @@ LSLogMA(; eta = 1, true_res = false ) = LSLogMA( cr, - MAInfo(1, lambda2, 1, false, 1, zeros(lambda2)), + MAInfo(lambda1, lambda2, lambda1, false, 1, zeros(lambda2)), Float64[], Float64[], Int64[], @@ -172,7 +172,7 @@ function log_update!( update_ma!(log, res, ma_info.lambda2, iter) else #Check if we can switch between lambda1 and lambda2 regime - if res < ma_info.res_window[ma_info.idx] + if iter == 0 || res <= ma_info.res_window[ma_info.idx] update_ma!(log, res, ma_info.lambda1, iter) else update_ma!(log, res, ma_info.lambda1, iter) @@ -209,7 +209,7 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: ma_info.idx = ma_info.idx < ma_info.lambda2 ? ma_info.idx + 1 : 1 ma_info.res_window[ma_info.idx] = res #Check if entire storage buffer can be used - if ma_info.lambda == lambda_base + if ma_info.lambda == ma_info.lambda2 # Compute the moving average for i in 1:ma_info.lambda2 accum += ma_info.res_window[i] @@ -221,6 +221,32 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: push!(log.resid_hist, accum / ma_info.lambda) push!(log.iota_hist, accum2 / ma_info.lambda) end + + elseif ma_info.lambda == ma_info.lambda1 && !ma_info.flag + diff = ma_info.idx - ma_info.lambda + # Determine start point for first loop + startp1 = diff < 0 ? 1 : (diff + 1) + + # Determine start and endpoints for second loop + startp2 = diff < 0 ? ma_info.lambda2 + diff + 1 : 2 + endp2 = diff < 0 ? ma_info.lambda2 : 1 + # Compute the moving average two loop setup required when lambda < lambda2 + for i in startp1:ma_info.idx + accum += ma_info.res_window[i] + accum2 += ma_info.res_window[i]^2 + end + + for i in startp2:endp2 + accum += ma_info.res_window[i] + accum2 += ma_info.res_window[i]^2 + end + + #Update the log variable with the information for this update + if mod(iter, log.collection_rate) == 0 || iter == 0 + push!(log.width_hist, ma_info.lambda) + push!(log.resid_hist, accum / ma_info.lambda) + push!(log.iota_hist, accum2 / ma_info.lambda) + end else # Get the difference between the start and current lambda @@ -351,6 +377,7 @@ for type in (LinSysVecRowGaussSampler, LinSysVecRowSparseGaussSampler) function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) log.dist_info.sigma2 = log.dist_info.block_dimension / (0.2345 * log.dist_info.eta) log.dist_info.omega = .1127 + log.dist_info.scaling = 1. end end diff --git a/test/linear_solver_logs/proc_solve_log_ma.jl b/test/linear_solver_logs/proc_solve_log_ma.jl index 830ff05..3d63b49 100644 --- a/test/linear_solver_logs/proc_solve_log_ma.jl +++ b/test/linear_solver_logs/proc_solve_log_ma.jl @@ -35,7 +35,7 @@ Random.seed!(1010) @test log.converged == false end - # Verify moving average + # Verify late moving average let A = rand(2,2) x = rand(2) @@ -56,10 +56,10 @@ Random.seed!(1010) obs_res = 2 .* [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^2 for i = 0:10] obs_res2 = 4 .* [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^4 for i = 0:10] @test length(log.resid_hist) == 11 - @test norm(log.resid_hist[2:11] - vcat(obs_res[2], - [(obs_res[i] + obs_res[i-1])/2 for i = 3:11])) < 1e2 * eps() - @test norm(log.iota_hist[2:11] - vcat(obs_res2[2], - [(obs_res2[i] + obs_res2[i-1])/2 for i = 3:11])) < 1e2 * eps() + @test norm(log.resid_hist[3:11] - vcat(obs_res[3], + [(obs_res[i] + obs_res[i-1])/2 for i = 4:11])) < 1e2 * eps() + @test norm(log.iota_hist[3:11] - vcat(obs_res2[3], + [(obs_res2[i] + obs_res2[i-1])/2 for i = 4:11])) < 1e2 * eps() @test log.iterations == 10 @test log.converged == false @@ -70,7 +70,40 @@ Random.seed!(1010) @test norm((Uncertainty_set[2] - log.resid_hist) ./ sqrt.(2 * Base.log(2/.05) * log.iota_hist * log.dist_info.sigma2 .* (1 .+ Base.log.(log.width_hist)) ./ log.width_hist) .- 1) < 1e2 * eps() end + # Verify early moving average + let + A = rand(2,2) + x = rand(2) + b = A * x + z = rand(2) + + sampler = LinSysVecRowOneRandCyclic() + log = LSLogMA(lambda1 = 2, + lambda2 = 10) + samp = (A[1,:], b[1]) + RLinearAlgebra.log_update!(log, sampler, z, (A[1,:],b[1]), 0, A, b) + # Test moving average of log when the residual only decreases + for i = 1:10 + samp = (A[1,:], b[1]) + RLinearAlgebra.log_update!(log, sampler, x + .3^(i+1)*(z-x), samp, i, A, b) + end + #compute sampled residuals + obs_res = 2 .* [abs(dot(A[1,:],x + .3^(i+1)*(z-x)) - b[1])^2 for i = 0:10] + obs_res2 = 4 .* [abs(dot(A[1,:],x + .3^(i+1)*(z-x)) - b[1])^4 for i = 0:10] + @test length(log.resid_hist) == 11 + @test norm(log.resid_hist[3:11] - vcat( [(obs_res[i] + obs_res[i-1])/2 for i = 3:11])) < 1e2 * eps() + @test norm(log.iota_hist[3:11] - vcat( [(obs_res2[i] + obs_res2[i-1])/2 for i = 3:11])) < 1e2 * eps() + @test log.iterations == 10 + @test log.converged == false + + #Test uncertainty set + Uncertainty_set = get_uncertainty(log) + @test length(Uncertainty_set[1]) == 11 + #If you undo the steps of the interval calculation should be 1 + @test norm((Uncertainty_set[2] - log.resid_hist) ./ sqrt.(2 * Base.log(2/.05) * log.iota_hist * + log.dist_info.sigma2 .* (1 .+ Base.log.(log.width_hist)) ./ log.width_hist) .- 1) < 1e2 * eps() + end end end # End Module From ebb9ceec81d260a8ca4dafdc603d2740caf0eb8f Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Fri, 9 Aug 2024 13:12:01 -0500 Subject: [PATCH 03/14] fixed tests to align with updates to master --- src/linear_solver_logs/solve_log_ma.jl | 20 +++---- test/linear_solver_logs/proc_solve_log_ma.jl | 60 ++++++++++---------- 2 files changed, 40 insertions(+), 40 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 51ef86a..b1ded39 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -130,7 +130,7 @@ LSLogMA(; omega = nothing, eta = 1, true_res = false - ) = LSLogMA( cr, + ) = LSLogMA( collection_rate, MAInfo(lambda1, lambda2, lambda1, false, 1, zeros(lambda2)), Float64[], Float64[], @@ -139,7 +139,7 @@ LSLogMA(; -1, false, true_res, - SEDistInfo(nothing, 0, 0, sigma2, omega, eta, 0 + SEDistInfo(nothing, 0, 0, sigma2, omega, eta, 0) ) #Function to update the moving average @@ -150,7 +150,7 @@ function log_update!( samp::Tuple, iter::Int64, A::AbstractArray, - b::AbstractVector + b::AbstractVector, ) if iter == 0 # Check if it is a row or column method and record dimensions @@ -184,8 +184,8 @@ function log_update!( if !log.true_res && iter > 0 # Compute the current residual to second power to align with theory # Check if it is one dimensional or block sampling method - res::Float64 = eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? - log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2 + res::Float64 = log.dist_info.scaling * (eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? + log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2) else res = log.resid_norm(A * x - b)^2 end @@ -268,7 +268,7 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: #Update the log variable with the information for this update if mod(iter, log.collection_rate) == 0 || iter == 0 - push!(log.width_hist, ma_info.lambda) + push!(log.lambda_hist, ma_info.lambda) push!(log.resid_hist, accum / ma_info.lambda) push!(log.iota_hist, accum2 / ma_info.lambda) end @@ -375,8 +375,8 @@ for type in (LinSysVecRowDetermCyclic,LinSysVecRowHopRandCyclic, LinSysVecRowMaxDistance,) @eval begin function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) - log.dist_info.sigma2 = log.dist_info.max_dimension^2 / (4 * log.dist_info.block_dimension^2 * log.dist_info.eta) - log.dist_info.scaling = log.dist_info.max_dimension / log.dist_info.block_dimension + log.dist_info.sigma2 = log.dist_info.dimension^2 / (4 * log.dist_info.block_dimension^2 * log.dist_info.eta) + log.dist_info.scaling = log.dist_info.dimension / log.dist_info.block_dimension end end @@ -388,8 +388,8 @@ end for type in (LinSysVecColOneRandCyclic, LinSysVecColDetermCyclic) @eval begin function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) - log.dist_info.sigma2 = log.dist_info.max_dimension^2 / (4 * log.dist_info.block_dimension^2 * log.dist_info.eta) - log.dist_info.scaling = log.dist_info.max_dimension / log.dist_info.block_dimension + log.dist_info.sigma2 = log.dist_info.dimension^2 / (4 * log.dist_info.block_dimension^2 * log.dist_info.eta) + log.dist_info.scaling = log.dist_info.dimension / log.dist_info.block_dimension end end diff --git a/test/linear_solver_logs/proc_solve_log_ma.jl b/test/linear_solver_logs/proc_solve_log_ma.jl index ff30296..c308756 100644 --- a/test/linear_solver_logs/proc_solve_log_ma.jl +++ b/test/linear_solver_logs/proc_solve_log_ma.jl @@ -23,16 +23,16 @@ Random.seed!(1010) z = rand(2) sampler = LinSysVecRowOneRandCyclic() - log = LSLogMA() + logger = LSLogMA() - RLinearAlgebra.log_update!(log, sampler, z, (A[1,:],b[1]), 0, A, b) + RLinearAlgebra.log_update!(logger, sampler, z, (A[1,:],b[1]), 0, A, b) - @test length(log.resid_hist) == 1 - @test log.resid_hist[1] == norm(A * z - b)^2 - @test norm(log.iota_hist[1] - norm(A * z - b)^4) < 1e2 * eps() - @test log.iterations == 0 - @test log.converged == false + @test length(logger.resid_hist) == 1 + @test logger.resid_hist[1] == norm(A * z - b)^2 + @test norm(logger.iota_hist[1] - norm(A * z - b)^4) < 1e2 * eps() + @test logger.iterations == 0 + @test logger.converged == false end # Verify late moving average @@ -43,32 +43,32 @@ Random.seed!(1010) z = rand(2) sampler = LinSysVecRowOneRandCyclic() - log = LSLogMA(lambda2 = 2) + logger = LSLogMA(lambda2 = 2) samp = (A[1,:], b[1]) - RLinearAlgebra.log_update!(log, sampler, z, (A[1,:],b[1]), 0, A, b) + RLinearAlgebra.log_update!(logger, sampler, z, (A[1,:],b[1]), 0, A, b) # Test moving average of log for i = 1:10 samp = (A[1,:], b[1]) - RLinearAlgebra.log_update!(log, sampler, x + (i+1)*(z-x), samp, i, A, b) + RLinearAlgebra.log_update!(logger, sampler, x + (i+1)*(z-x), samp, i, A, b) end #compute sampled residuals obs_res = 2 .* [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^2 for i = 0:10] obs_res2 = 4 .* [abs(dot(A[1,:],x + (i+1)*(z-x)) - b[1])^4 for i = 0:10] - @test length(log.resid_hist) == 11 - @test norm(log.resid_hist[3:11] - vcat(obs_res[3], + @test length(logger.resid_hist) == 11 + @test norm(logger.resid_hist[3:11] - vcat(obs_res[3], [(obs_res[i] + obs_res[i-1])/2 for i = 4:11])) < 1e2 * eps() - @test norm(log.iota_hist[3:11] - vcat(obs_res2[3], + @test norm(logger.iota_hist[3:11] - vcat(obs_res2[3], [(obs_res2[i] + obs_res2[i-1])/2 for i = 4:11])) < 1e2 * eps() - @test log.iterations == 10 - @test log.converged == false + @test logger.iterations == 10 + @test logger.converged == false #Test uncertainty set - Uncertainty_set = get_uncertainty(log) + Uncertainty_set = get_uncertainty(logger) @test length(Uncertainty_set[1]) == 11 #If you undo the steps of the interval calculation should be 1 - @test norm((Uncertainty_set[2] - log.resid_hist) ./ sqrt.(2 * Base.log(2/.05) * log.iota_hist * - log.dist_info.sigma2 .* (1 .+ Base.log.(log.lambda_hist)) ./ log.lambda_hist) .- 1) < 1e2 * eps() + @test norm((Uncertainty_set[2] - logger.resid_hist) ./ sqrt.(2 * log(2/.05) * logger.iota_hist * + logger.dist_info.sigma2 .* (1 .+ log.(logger.lambda_hist)) ./ logger.lambda_hist) .- 1) < 1e2 * eps() end # Verify early moving average let @@ -78,31 +78,31 @@ Random.seed!(1010) z = rand(2) sampler = LinSysVecRowOneRandCyclic() - log = LSLogMA(lambda1 = 2, + logger = LSLogMA(lambda1 = 2, lambda2 = 10) samp = (A[1,:], b[1]) - RLinearAlgebra.log_update!(log, sampler, z, (A[1,:],b[1]), 0, A, b) - # Test moving average of log when the residual only decreases + RLinearAlgebra.log_update!(logger, sampler, z, (A[1,:],b[1]), 0, A, b) + # Test moving average of log when the residual only decreases to not trigger switch for i = 1:10 samp = (A[1,:], b[1]) - RLinearAlgebra.log_update!(log, sampler, x + .3^(i+1)*(z-x), samp, i, A, b) + RLinearAlgebra.log_update!(logger, sampler, x + .3^(i+1)*(z-x), samp, i, A, b) end #compute sampled residuals obs_res = 2 .* [abs(dot(A[1,:],x + .3^(i+1)*(z-x)) - b[1])^2 for i = 0:10] obs_res2 = 4 .* [abs(dot(A[1,:],x + .3^(i+1)*(z-x)) - b[1])^4 for i = 0:10] - @test length(log.resid_hist) == 11 - @test norm(log.resid_hist[3:11] - vcat( [(obs_res[i] + obs_res[i-1])/2 for i = 3:11])) < 1e2 * eps() - @test norm(log.iota_hist[3:11] - vcat( [(obs_res2[i] + obs_res2[i-1])/2 for i = 3:11])) < 1e2 * eps() - @test log.iterations == 10 - @test log.converged == false + @test length(logger.resid_hist) == 11 + @test norm(logger.resid_hist[3:11] - vcat( [(obs_res[i] + obs_res[i-1])/2 for i = 3:11])) < 1e2 * eps() + @test norm(logger.iota_hist[3:11] - vcat( [(obs_res2[i] + obs_res2[i-1])/2 for i = 3:11])) < 1e2 * eps() + @test logger.iterations == 10 + @test logger.converged == false #Test uncertainty set - Uncertainty_set = get_uncertainty(log) + Uncertainty_set = get_uncertainty(logger) @test length(Uncertainty_set[1]) == 11 #If you undo the steps of the interval calculation should be 1 - @test norm((Uncertainty_set[2] - log.resid_hist) ./ sqrt.(2 * Base.log(2/.05) * log.iota_hist * - log.dist_info.sigma2 .* (1 .+ Base.log.(log.width_hist)) ./ log.width_hist) .- 1) < 1e2 * eps() + @test norm((Uncertainty_set[2] - logger.resid_hist) ./ sqrt.(2 * log(2/.05) * logger.iota_hist * + logger.dist_info.sigma2 .* (1 .+ log.(logger.lambda_hist)) ./ logger.lambda_hist) .- 1) < 1e2 * eps() end end From dd2e5457894a40e985037df6088f1d085db8b11d Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Mon, 19 Aug 2024 17:13:24 -0400 Subject: [PATCH 04/14] Updated comments to clarify two phase processing --- src/linear_solver_logs/solve_log_ma.jl | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index b1ded39..2f2db53 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -47,8 +47,8 @@ A mutable structure that stores information about the sub-Exponential family. if not specified by the user it is provided based on the sampling method. - `omega::Union{Float64, Nothing}`, the exponential distrbution parameter, if not specified by the user, it is provided based on the sampling method. -- `eta::Float64`, a parameter for adjusting the conservativeness of the distribution, higher value means a less conservative - estimate. By default, this is set to `1`. +- `eta::Float64`, a parameter for adjusting the conservativeness of the distribution, higher value means a less conservative estimate. By default, this is set to `1`. +- `scaling::Float64`, a scaling parameter for the norm-squared of the sketched residual to ensure its expectation is the norm-squared of the residual. For more information see: - Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. @@ -185,7 +185,7 @@ function log_update!( # Compute the current residual to second power to align with theory # Check if it is one dimensional or block sampling method res::Float64 = log.dist_info.scaling * (eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? - log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2) + log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2) else res = log.resid_norm(A * x - b)^2 end @@ -194,7 +194,9 @@ function log_update!( if ma_info.flag update_ma!(log, res, ma_info.lambda2, iter) else - #Check if we can switch between lambda1 and lambda2 regime + # Check if we can switch between lambda1 and lambda2 regime + # If it is in the monotonic decreasing of the sketched residual then we are in a lambda1 regime + # otherwise we switch to the lambda2 regime which is indicated by the changing of the flag if iter == 0 || res <= ma_info.res_window[ma_info.idx] update_ma!(log, res, ma_info.lambda1, iter) else @@ -246,13 +248,19 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: push!(log.resid_hist, accum / ma_info.lambda) push!(log.iota_hist, accum2 / ma_info.lambda) end - + + # In the lambda1 regime we are not summing all the entries in the ma and must be careful elseif ma_info.lambda == ma_info.lambda1 && !ma_info.flag diff = ma_info.idx - ma_info.lambda - # Determine start point for first loop + # Because lambda1 is less than lambda2 and the index storing the new norm of the residual + # moves sequentially there could be an issue of summing at the end and the begining of the buffer. + # When lambda is greater than the index we will have to go to end of the buffer to get the remaining + # terms in the sum and add it to the terms summed from index 1 to current index. Otherwise we can + # just start the sum lambda terms in front of the index startp1 = diff < 0 ? 1 : (diff + 1) - - # Determine start and endpoints for second loop + # Since the index is less than the width, we must add the remianing terms to the sum summing the + # terms starting at the end of the buffer-remaining terms to the length of the buffer. If this is not + # necessary we skip the second loop startp2 = diff < 0 ? ma_info.lambda2 + diff + 1 : 2 endp2 = diff < 0 ? ma_info.lambda2 : 1 # Compute the moving average two loop setup required when lambda < lambda2 From 93056cfa2b3491f86d2377bb416472e91450a2f7 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Wed, 21 Aug 2024 15:05:07 -0400 Subject: [PATCH 05/14] removed sampler lines from code --- src/linear_solver_logs/solve_log_ma.jl | 24 ++++++++------------ src/linear_solver_stops/stop_ma_thres.jl | 6 ++--- test/linear_solver_logs/proc_solve_log_ma.jl | 10 ++++---- 3 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 2f2db53..f514917 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -8,7 +8,6 @@ A mutable structure that stores information relevant to the moving average of th - `lambda1::Int64`, the width of the moving average during the fast convergence phase of the algorithm. During this fast convergence phase, the majority of variation of the sketched estimator comes from improvement in the solution and thus wide moving average windows inaccurately represent progress. - By default this parameter is set to `1`. - `lambda2::Int64`, the width of the moving average in the slower convergence phase. In the slow convergence phase, each iterate differs from the previous one by a small amount and thus most of the observed variation arises from the randomness of the sketched progress estimator, which is best smoothed by a wide moving @@ -92,7 +91,7 @@ Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertain. Quantifica calculation, it is collected at a rate specified by `collection_rate`. - `resid_norm::Function`, a function that accepts a single vector argument and returns a scalar. Used to compute the residual size. -- `iterations::Int64`, the current iteration of the solver. +- `iteration::Int64`, the current iteration of the solver. - `converged::Bool`, a flag to indicate whether the system has converged by some measure. By default this is set to false. - `true_res::Bool`, a boolean indicating if we want the true residual computed instead of approximate. @@ -116,7 +115,7 @@ mutable struct LSLogMA <: LinSysSolverLog iota_hist::Vector{Float64} lambda_hist::Vector{Int64} resid_norm::Function - iterations::Int64 + iteration::Int64 converged::Bool true_res::Bool dist_info::SEDistInfo @@ -170,16 +169,15 @@ function log_update!( log.dist_info.block_dimension = sampler.block_size end - log.dist_info.sampler = typeof(sampler) - # If the constants for the sub-Exponential distribution are not defined then define them + # If the constants for the sub-Exponential distribution are not defined then define them if typeof(log.dist_info.sigma2) <: Nothing || log.dist_info.sigma2 == 0 - get_SE_constants!(log, log.dist_info.sampler) + get_SE_constants!(log, typeof(sampler)) end end ma_info = log.ma_info - log.iterations = iter + log.iteration = iter #Check if we want exact residuals computed if !log.true_res && iter > 0 # Compute the current residual to second power to align with theory @@ -197,12 +195,10 @@ function log_update!( # Check if we can switch between lambda1 and lambda2 regime # If it is in the monotonic decreasing of the sketched residual then we are in a lambda1 regime # otherwise we switch to the lambda2 regime which is indicated by the changing of the flag - if iter == 0 || res <= ma_info.res_window[ma_info.idx] - update_ma!(log, res, ma_info.lambda1, iter) - else - update_ma!(log, res, ma_info.lambda1, iter) - ma_info.flag = true - end + # because update_ma changes res_window and ma_info.idx must check condition first + flag_cond = iter == 0 || res <= ma_info.res_window[ma_info.idx] + update_ma!(log, res, ma_info.lambda1, iter) + ma_info.flag = flag_cond ? false : true end @@ -328,7 +324,7 @@ function get_uncertainty(hist::LSLogMA; alpha::Float64 = 0.05) lower = zeros(l) # If the constants for the sub-Exponential distribution are not defined then define them if typeof(hist.dist_info.sigma2) <: Nothing - get_SE_constants!(hist, hist.dist_info.sampler) + throw(ArgumentError("The SE constants are empty, please set them in dist_info field of LSLogMA first.")) end for i in 1:l diff --git a/src/linear_solver_stops/stop_ma_thres.jl b/src/linear_solver_stops/stop_ma_thres.jl index 4cd6a11..c37e862 100644 --- a/src/linear_solver_stops/stop_ma_thres.jl +++ b/src/linear_solver_stops/stop_ma_thres.jl @@ -39,7 +39,7 @@ function check_stop_criterion( log::LSLogMA, stop::LSStopMA ) - its = log.iterations + its = log.iteration if its > 0 I_threshold = iota_threshold(log, stop) thresholdChecks = sqrt(log.iota_hist[its]) <= I_threshold && @@ -81,10 +81,10 @@ function iota_threshold(hist::LSLogMA, stop::LSStopMA) # Compute the threshold bound in the case where there is no omega c = min((1 - delta1)^2 * threshold^2 / (2 * log(1/chi1)), (delta2 - 1)^2 * threshold^2 / (2 * log(1/chi2))) - c /= (hist.dist_info.sigma2 * sqrt(hist.iota_hist[hist.iterations])) * (1 + log(lambda)) / lambda + c /= (hist.dist_info.sigma2 * sqrt(hist.iota_hist[hist.iteration])) * (1 + log(lambda)) / lambda else #compute error bound when there is an omega - siota = (hist.dist_info.sigma2 * sqrt(hist.iota_hist[hist.iterations])) * (1 + log(lambda)) / lambda + siota = (hist.dist_info.sigma2 * sqrt(hist.iota_hist[hist.iteration])) * (1 + log(lambda)) / lambda min1 = min((1 - delta1)^2 * threshold^2 / (2 * log(1/chi1) * siota), lambda * (1 - delta1) * threshold / (2 * log(1/chi1) * hist.dist_info.omega)) min2 = min((delta2 - 1)^2 * threshold^2 / (2 * log(1/chi2) * siota), diff --git a/test/linear_solver_logs/proc_solve_log_ma.jl b/test/linear_solver_logs/proc_solve_log_ma.jl index c308756..2543c3a 100644 --- a/test/linear_solver_logs/proc_solve_log_ma.jl +++ b/test/linear_solver_logs/proc_solve_log_ma.jl @@ -12,7 +12,7 @@ Random.seed!(1010) @test supertype(LSLogMA) == LinSysSolverLog # Verify Required Fields - @test :iterations in fieldnames(LSLogMA) + @test :iteration in fieldnames(LSLogMA) @test :converged in fieldnames(LSLogMA) # Verify log_update initialization @@ -24,14 +24,14 @@ Random.seed!(1010) sampler = LinSysVecRowOneRandCyclic() logger = LSLogMA() - + @test_throws ArgumentError("The SE constants are empty, please set them in dist_info field of LSLogMA first.") get_uncertainty(logger) RLinearAlgebra.log_update!(logger, sampler, z, (A[1,:],b[1]), 0, A, b) @test length(logger.resid_hist) == 1 @test logger.resid_hist[1] == norm(A * z - b)^2 @test norm(logger.iota_hist[1] - norm(A * z - b)^4) < 1e2 * eps() - @test logger.iterations == 0 + @test logger.iteration == 0 @test logger.converged == false end @@ -60,7 +60,7 @@ Random.seed!(1010) [(obs_res[i] + obs_res[i-1])/2 for i = 4:11])) < 1e2 * eps() @test norm(logger.iota_hist[3:11] - vcat(obs_res2[3], [(obs_res2[i] + obs_res2[i-1])/2 for i = 4:11])) < 1e2 * eps() - @test logger.iterations == 10 + @test logger.iteration == 10 @test logger.converged == false #Test uncertainty set @@ -94,7 +94,7 @@ Random.seed!(1010) @test length(logger.resid_hist) == 11 @test norm(logger.resid_hist[3:11] - vcat( [(obs_res[i] + obs_res[i-1])/2 for i = 3:11])) < 1e2 * eps() @test norm(logger.iota_hist[3:11] - vcat( [(obs_res2[i] + obs_res2[i-1])/2 for i = 3:11])) < 1e2 * eps() - @test logger.iterations == 10 + @test logger.iteration == 10 @test logger.converged == false #Test uncertainty set From 774f92c34db85fa15180574e90363b930c3ef2c0 Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Wed, 21 Aug 2024 15:08:06 -0400 Subject: [PATCH 06/14] Corrected stopping test issue relating to change in iterations terminology --- test/linear_solver_stops/proc_stop_ma.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/linear_solver_stops/proc_stop_ma.jl b/test/linear_solver_stops/proc_stop_ma.jl index f981d6e..dc577aa 100644 --- a/test/linear_solver_stops/proc_stop_ma.jl +++ b/test/linear_solver_stops/proc_stop_ma.jl @@ -17,13 +17,13 @@ using Test, RLinearAlgebra log.dist_info.dimension = 100 log.dist_info.sampler = LinSysVecRowDetermCyclic - log.iterations = 0 + log.iteration = 0 @test RLinearAlgebra.check_stop_criterion(log, stop) == false - log.iterations = 2 + log.iteration = 2 @test RLinearAlgebra.check_stop_criterion(log, stop) == true - log.iterations = 3 + log.iteration = 3 @test RLinearAlgebra.check_stop_criterion(log, stop) == false #Verify threshold stopping @@ -33,13 +33,13 @@ using Test, RLinearAlgebra log.ma_info.lambda = 15 stop = LSStopMA(4, 1e-10, 1.1, .9, .01, .01) - log.iterations = 1 + log.iteration = 1 @test RLinearAlgebra.check_stop_criterion(log, stop) == false - log.iterations = 2 + log.iteration = 2 @test RLinearAlgebra.check_stop_criterion(log, stop) == false - log.iterations = 3 + log.iteration = 3 @test RLinearAlgebra.check_stop_criterion(log, stop) == true end From ce01ef866ccee40c3762609cfb5bc14700d8cfc0 Mon Sep 17 00:00:00 2001 From: Vivak Date: Wed, 28 Aug 2024 11:41:42 -0400 Subject: [PATCH 07/14] Removed default for `lambda2` --- src/linear_solver_logs/solve_log_ma.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index f514917..6dfd91a 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -11,7 +11,7 @@ A mutable structure that stores information relevant to the moving average of th - `lambda2::Int64`, the width of the moving average in the slower convergence phase. In the slow convergence phase, each iterate differs from the previous one by a small amount and thus most of the observed variation arises from the randomness of the sketched progress estimator, which is best smoothed by a wide moving - average width. By default this parameter is set to `30`. + average width. - `lambda::Int64`, the width of the moving average at the current iteration. This value is not controlled by the user. - `flag::Bool`, a boolean indicating which phase we are in, a value of `true` indicates second phase. From cf85d7bb4e1e4181fd0909d07cf0f1d65ad9aa1e Mon Sep 17 00:00:00 2001 From: Vivak Date: Wed, 28 Aug 2024 11:43:15 -0400 Subject: [PATCH 08/14] Changed `flag` description to slow convergence phase --- src/linear_solver_logs/solve_log_ma.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 6dfd91a..3ce4626 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -14,7 +14,7 @@ A mutable structure that stores information relevant to the moving average of th average width. - `lambda::Int64`, the width of the moving average at the current iteration. This value is not controlled by the user. -- `flag::Bool`, a boolean indicating which phase we are in, a value of `true` indicates second phase. +- `flag::Bool`, a boolean indicating which phase we are in, a value of `true` indicates slow convergence phase. - `idx::Int64`, the index indcating what value should be replaced in the moving average buffer. - `res_window::Vector{Float64}`, the moving average buffer. From 6ab723b21cbbbb6f4ee3a6c17454ebcc79495229 Mon Sep 17 00:00:00 2001 From: Vivak Date: Wed, 28 Aug 2024 11:46:12 -0400 Subject: [PATCH 09/14] Fixed reference formats for MAInfo and SEDistInfo --- src/linear_solver_logs/solve_log_ma.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 3ce4626..954b23c 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -19,10 +19,11 @@ A mutable structure that stores information relevant to the moving average of th - `res_window::Vector{Float64}`, the moving average buffer. For more information see: - - Pritchard, Nathaniel, and Patel, Vivak. Solving, tracking and stopping streaming linear inverse problems. -Inverse Problems 40.8 Web. doi:10.1088/1361-6420/ad5583. - - Pritchard, Nathaniel and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative Least Squares -Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertain. Quantification 11 (2022): 996-1024. +- Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear +inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. +- Pritchard, Nathaniel, and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative +Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty +Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 """ mutable struct MAInfo lambda1::Int64 @@ -50,9 +51,11 @@ A mutable structure that stores information about the sub-Exponential family. - `scaling::Float64`, a scaling parameter for the norm-squared of the sketched residual to ensure its expectation is the norm-squared of the residual. For more information see: -- Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. -- Pritchard, Nathaniel and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative Least Squares -Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertain. Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 +- Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear +inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. +- Pritchard, Nathaniel, and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative +Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty +Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 """ mutable struct SEDistInfo sampler::Union{DataType, Nothing} From 32c4369751ab19cb31584af2a092d35a07be8754 Mon Sep 17 00:00:00 2001 From: Vivak Date: Wed, 28 Aug 2024 11:58:40 -0400 Subject: [PATCH 10/14] Edits to docstring for SEDistInfo --- src/linear_solver_logs/solve_log_ma.jl | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 954b23c..a7844d4 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -37,18 +37,23 @@ end """ SEDistInfo -A mutable structure that stores information about the sub-Exponential family. +A mutable structure that stores information about a distribution (i.e., sampling method) +in the sub-Exponential family. # Fields - `sampler::Union{DataType, Nothing}`, the type of sampling method. - `dimension::Int64`, the dimension that of the space that is being sampled. - `block_dimension::Int64`, the dimension of the sample. -- `sigma2::Union{Float64, Nothing}`, the variance parameter in the sub-Exponential family, - if not specified by the user it is provided based on the sampling method. -- `omega::Union{Float64, Nothing}`, the exponential distrbution parameter, if not specified by the user, - it is provided based on the sampling method. -- `eta::Float64`, a parameter for adjusting the conservativeness of the distribution, higher value means a less conservative estimate. By default, this is set to `1`. -- `scaling::Float64`, a scaling parameter for the norm-squared of the sketched residual to ensure its expectation is the norm-squared of the residual. +- `sigma2::Union{Float64, Nothing}`, the variance parameter in the sub-Exponential family. + If not specified by the user, a value is selected from a table based on the `sampler`. + If the `sampler` is not in the table, then ?. +- `omega::Union{Float64, Nothing}`, the exponential distrbution parameter. If not specified + by the user, a value is selected from a table based on the `sampler`. + If the `sampler` is not in the table, then ?. +- `eta::Float64`, a parameter for adjusting the conservativeness of the distribution, higher + value means a less conservative estimate. A recommended value is `1`. +- `scaling::Float64`, a scaling parameter for the norm-squared of the sketched residual to + ensure its expectation is the norm-squared of the residual. For more information see: - Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear From 9f976f0a91112839f78b646a5ae339ade85c94ad Mon Sep 17 00:00:00 2001 From: Vivak Date: Wed, 28 Aug 2024 12:19:13 -0400 Subject: [PATCH 11/14] Fixed references to ensure formatted correctly in docs --- src/linear_solver_logs/solve_log_ma.jl | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index a7844d4..06c21d2 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -20,10 +20,10 @@ A mutable structure that stores information relevant to the moving average of th For more information see: - Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear -inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. + inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. - Pritchard, Nathaniel, and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative -Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty -Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 + Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty + Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 """ mutable struct MAInfo lambda1::Int64 @@ -57,10 +57,10 @@ in the sub-Exponential family. For more information see: - Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear -inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. + inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. - Pritchard, Nathaniel, and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative -Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty -Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 + Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty + Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 """ mutable struct SEDistInfo sampler::Union{DataType, Nothing} @@ -80,11 +80,6 @@ The log assumes that only a random block of full linear system is available for The goal of this log is usually for all use cases as it acts as a cheap approximation of the full residual. -For more information see: -- Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. -- Pritchard, Nathaniel and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative Least Squares -Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertain. Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 - # Fields - `collection_rate::Int64`, the frequency with which to record information about progress estimators to append to the remaining fields, starting with the initialization (i.e., iterate `0`). @@ -115,6 +110,13 @@ Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertain. Quantifica omega = nothing, eta = 1, true_res = false)` + +For more information see: +- Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear + inverse problems." Inverse Problems (2024). doi:10.1088/1361-6420/ad5583. +- Pritchard, Nathaniel, and Vivak Patel. “Towards Practical Large-Scale Randomized Iterative + Least Squares Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertainty + Quantification 11 (2022): 996-1024. doi.org/10.1137/22M1515057 """ mutable struct LSLogMA <: LinSysSolverLog collection_rate::Int64 From 1bb8056db285ea85084024be7cfc172727f1d2b7 Mon Sep 17 00:00:00 2001 From: Vivak Date: Wed, 28 Aug 2024 12:34:48 -0400 Subject: [PATCH 12/14] Updates to docstring for `LSLogMa` --- src/linear_solver_logs/solve_log_ma.jl | 31 +++++++++++++------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 06c21d2..a10726d 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -75,28 +75,27 @@ end """ LSLogMA <: LinSysSolverLog -A mutable structure that stores information about a randomized linear solver's behavior. -The log assumes that only a random block of full linear system is available for processing. -The goal of this log is usually for all use cases as it acts as a cheap approximation of the -full residual. +A mutable structure that stores information for tracking a randomized linear solver's + progress. + The log assumes that the entire linear system is **not** available for processing. + # Fields - `collection_rate::Int64`, the frequency with which to record information about progress estimators to append to the remaining fields, starting with the initialization (i.e., iterate `0`). - `ma_info::MAInfo`, [`MAInfo`](@ref) -- `resid_hist::Vector{Float64}`, a structure that contains the moving average of the error proxy squared, - typically the norm of residual or gradient of normal equations, it is collected at a rate - specified by `collection_rate`. -- `iota_hist::Vector{Float64}`, a structure that contains the moving average of the error proxy - to the fourth power, typically the norm of residual or gradient of normal equations. This is used - in part to approximate the variance of the estimator, it is collected at a rate specified by `collection_rate`. -- `lambda_hist::Vector{Int64}`, data structure that contains the lambdas, widths of the moving average, - calculation, it is collected at a rate specified by `collection_rate`. -- `resid_norm::Function`, a function that accepts a single vector argument and returns a - scalar. Used to compute the residual size. +- `resid_hist::Vector{Float64}`, contains an estimate of the progress of the randomized + solver. These values are stored at iterates specified by `collection_rate`. +- `iota_hist::Vector{Float64}`, contains an estimate used for calculating the variability + of the progress estimators. These values are stored at iterates specified by + `collection_rate`. +- `lambda_hist::Vector{Int64}`, contains the widths of the moving average. + These values are stored at iterates specified by `collection_rate`. +- `resid_norm::Function`, the desired `norm` used. The default constructor sets this to the + Euclidean norm. - `iteration::Int64`, the current iteration of the solver. -- `converged::Bool`, a flag to indicate whether the system has converged by some measure. By default this - is set to false. +- `converged::Bool`, a flag to indicate whether the system has converged by some measure. + By default this is `false`. - `true_res::Bool`, a boolean indicating if we want the true residual computed instead of approximate. - `dist_info::SEDistInfo`, [`SEDistInfo`](@ref) From f9088aa7bc2f5eb6a64686564cb44e4bbec2787e Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Wed, 28 Aug 2024 16:55:10 -0400 Subject: [PATCH 13/14] Updated default metric for get uncertainty --- src/linear_solver_logs/solve_log_ma.jl | 29 ++++++++------------ test/linear_solver_logs/proc_solve_log_ma.jl | 18 ++++++++++-- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index a10726d..63e1a8c 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -96,7 +96,6 @@ A mutable structure that stores information for tracking a randomized linear sol - `iteration::Int64`, the current iteration of the solver. - `converged::Bool`, a flag to indicate whether the system has converged by some measure. By default this is `false`. -- `true_res::Bool`, a boolean indicating if we want the true residual computed instead of approximate. - `dist_info::SEDistInfo`, [`SEDistInfo`](@ref) # Constructors @@ -108,7 +107,7 @@ A mutable structure that stores information for tracking a randomized linear sol sigma2 = nothing, omega = nothing, eta = 1, - true_res = false)` + )` For more information see: - Pritchard, Nathaniel, and Vivak Patel. "Solving, tracking and stopping streaming linear @@ -126,7 +125,6 @@ mutable struct LSLogMA <: LinSysSolverLog resid_norm::Function iteration::Int64 converged::Bool - true_res::Bool dist_info::SEDistInfo end @@ -137,7 +135,6 @@ LSLogMA(; sigma2 = nothing, omega = nothing, eta = 1, - true_res = false ) = LSLogMA( collection_rate, MAInfo(lambda1, lambda2, lambda1, false, 1, zeros(lambda2)), Float64[], @@ -146,7 +143,6 @@ LSLogMA(; norm, -1, false, - true_res, SEDistInfo(nothing, 0, 0, sigma2, omega, eta, 0) ) @@ -182,20 +178,15 @@ function log_update!( if typeof(log.dist_info.sigma2) <: Nothing || log.dist_info.sigma2 == 0 get_SE_constants!(log, typeof(sampler)) end - + end ma_info = log.ma_info log.iteration = iter - #Check if we want exact residuals computed - if !log.true_res && iter > 0 - # Compute the current residual to second power to align with theory - # Check if it is one dimensional or block sampling method - res::Float64 = log.dist_info.scaling * (eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? - log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2) - else - res = log.resid_norm(A * x - b)^2 - end + # Compute the current residual to second power to align with theory + # Check if it is one dimensional or block sampling method + res::Float64 = log.dist_info.scaling * (eltype(samp[1]) <: Int64 || size(samp[1],2) != 1 ? + log.resid_norm(samp[3])^2 : log.resid_norm(dot(samp[1], x) - samp[2])^2) # Check if MA is in lambda1 or lambda2 regime if ma_info.flag @@ -238,7 +229,7 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: accum = 0 accum2 = 0 ma_info = log.ma_info - ma_info.idx = ma_info.idx < ma_info.lambda2 ? ma_info.idx + 1 : 1 + ma_info.idx = ma_info.idx < ma_info.lambda2 && iter != 0 ? ma_info.idx + 1 : 1 ma_info.res_window[ma_info.idx] = res #Check if entire storage buffer can be used if ma_info.lambda == ma_info.lambda2 @@ -374,10 +365,12 @@ This function is not exported and thus the user does not have direct access to i # Outputs Performs an inplace update of the sub-Exponential constants for the log. Additionally, updates the scaling constant to ensure expectation of -block norms is equal to true norm. +block norms is equal to true norm. If default is not a defined a warning is returned that sigma2 is set 1 and scaling is set to 1. """ function get_SE_constants!(log::LSLogMA, sampler::Type{T}) where T<:LinSysSampler - return nothing + @warn "No constants defined for method of type $sampler. By default we set sigma2 to 1 and scaling to 1." + log.dist_info.sigma2 = 1 + log.dist_info.scaling = 1 end for type in (LinSysVecRowDetermCyclic,LinSysVecRowHopRandCyclic, diff --git a/test/linear_solver_logs/proc_solve_log_ma.jl b/test/linear_solver_logs/proc_solve_log_ma.jl index 2543c3a..8e56ef3 100644 --- a/test/linear_solver_logs/proc_solve_log_ma.jl +++ b/test/linear_solver_logs/proc_solve_log_ma.jl @@ -24,15 +24,23 @@ Random.seed!(1010) sampler = LinSysVecRowOneRandCyclic() logger = LSLogMA() + + # Test the error message on the get_uncertainty function @test_throws ArgumentError("The SE constants are empty, please set them in dist_info field of LSLogMA first.") get_uncertainty(logger) + # Test warning message in default SEconstant look up - RLinearAlgebra.log_update!(logger, sampler, z, (A[1,:],b[1]), 0, A, b) + RLinearAlgebra.log_update!(logger, sampler, z, (A[1, :], b[1]), 0, A, b) @test length(logger.resid_hist) == 1 - @test logger.resid_hist[1] == norm(A * z - b)^2 - @test norm(logger.iota_hist[1] - norm(A * z - b)^4) < 1e2 * eps() + @test norm(logger.resid_hist[1] - 2 * norm(dot(A[1, :], z) - b[1])^2) < 1e2 * eps() + @test norm(logger.iota_hist[1] - 4 * norm(dot(A[1, :], z) - b[1])^4) < 1e2 * eps() @test logger.iteration == 0 @test logger.converged == false + + struct MadeUpSampler <: LinSysSampler + end + @test_logs (:warn, "No constants defined for method of type Main.ProceduralTestLSLogMA.MadeUpSampler. By default we set sigma2 to 1 and scaling to 1.") RLinearAlgebra.get_SE_constants!(logger, MadeUpSampler) + end # Verify late moving average @@ -104,6 +112,10 @@ Random.seed!(1010) @test norm((Uncertainty_set[2] - logger.resid_hist) ./ sqrt.(2 * log(2/.05) * logger.iota_hist * logger.dist_info.sigma2 .* (1 .+ log.(logger.lambda_hist)) ./ logger.lambda_hist) .- 1) < 1e2 * eps() end + + + + end end # End Module From 7c5e3d3ea0971b06133956334792f7e16598959e Mon Sep 17 00:00:00 2001 From: Nathaniel pritchard Date: Thu, 5 Sep 2024 10:15:05 -0400 Subject: [PATCH 14/14] Removed unnecessary condition added error for sampler type and added to test --- src/linear_solver_logs/solve_log_ma.jl | 83 +++++++------------- test/linear_solver_logs/proc_solve_log_ma.jl | 3 +- 2 files changed, 30 insertions(+), 56 deletions(-) diff --git a/src/linear_solver_logs/solve_log_ma.jl b/src/linear_solver_logs/solve_log_ma.jl index 63e1a8c..6a96513 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -46,10 +46,10 @@ in the sub-Exponential family. - `block_dimension::Int64`, the dimension of the sample. - `sigma2::Union{Float64, Nothing}`, the variance parameter in the sub-Exponential family. If not specified by the user, a value is selected from a table based on the `sampler`. - If the `sampler` is not in the table, then ?. + If the `sampler` is not in the table, then `sigma2` is set to `1`. - `omega::Union{Float64, Nothing}`, the exponential distrbution parameter. If not specified by the user, a value is selected from a table based on the `sampler`. - If the `sampler` is not in the table, then ?. + If the `sampler` is not in the table, then `omega` is set to `1`. - `eta::Float64`, a parameter for adjusting the conservativeness of the distribution, higher value means a less conservative estimate. A recommended value is `1`. - `scaling::Float64`, a scaling parameter for the norm-squared of the sketched residual to @@ -158,20 +158,19 @@ function log_update!( ) if iter == 0 # Check if it is a row or column method and record dimensions + log.dist_info.dimension = size(A,1) if supertype(typeof(sampler)) <: LinSysVecRowSampler - log.dist_info.dimension = size(A,1) log.dist_info.block_dimension = 1 elseif supertype(typeof(sampler)) <: LinSysBlkRowSampler - log.dist_info.dimension = size(A,1) # For the block methods samp[3] is always the sketched residual, its length is block size log.dist_info.block_dimension = sampler.block_size elseif supertype(typeof(sampler)) <: LinSysVecColSampler - log.dist_info.dimension = size(A,1) log.dist_info.block_dimension = 1 - else - log.dist_info.dimension = size(A,1) + elseif supertype(typeof(sampler)) <: LinSysBlkColSampler # For the block methods samp[3] is always the sketched residual, its length is block size log.dist_info.block_dimension = sampler.block_size + else + throw(ArgumentError("`sampler` is not of type `LinSysBlkColSampler`, `LinSysVecColSampler`, `LinSysBlkRowSampler`, or `LinSysVecRowSampler`")) end # If the constants for the sub-Exponential distribution are not defined then define them @@ -195,10 +194,10 @@ function log_update!( # Check if we can switch between lambda1 and lambda2 regime # If it is in the monotonic decreasing of the sketched residual then we are in a lambda1 regime # otherwise we switch to the lambda2 regime which is indicated by the changing of the flag - # because update_ma changes res_window and ma_info.idx must check condition first + # because update_ma changes res_window and ma_info.idx we must check condition first flag_cond = iter == 0 || res <= ma_info.res_window[ma_info.idx] update_ma!(log, res, ma_info.lambda1, iter) - ma_info.flag = flag_cond ? false : true + ma_info.flag = !flag_cond end @@ -217,8 +216,7 @@ Function that updates the moving average tracking statistic. # Inputs - `log::LSLogMA`, the moving average log structure. -- `res::Union{AbstractVector, Real}, the residual for the current iteration. - This could be sketeched or full residual depending on the inputs when creating the log structor. +- `res::Union{AbstractVector, Real}`, the sketched residual for the current iteration. - `lambda_base::Int64`, which lambda, between lambda1 and lambda2, is currently being used. - `iter::Int64`, the current iteration. @@ -226,7 +224,9 @@ Function that updates the moving average tracking statistic. Updates the log datatype and does not explicitly return anything. """ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base::Int64, iter::Int64) + # Variable to store the sum of the terms for rho accum = 0 + # Variable to store the sum of the terms for iota accum2 = 0 ma_info = log.ma_info ma_info.idx = ma_info.idx < ma_info.lambda2 && iter != 0 ? ma_info.idx + 1 : 1 @@ -245,48 +245,20 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: push!(log.iota_hist, accum2 / ma_info.lambda) end - # In the lambda1 regime we are not summing all the entries in the ma and must be careful - elseif ma_info.lambda == ma_info.lambda1 && !ma_info.flag - diff = ma_info.idx - ma_info.lambda - # Because lambda1 is less than lambda2 and the index storing the new norm of the residual - # moves sequentially there could be an issue of summing at the end and the begining of the buffer. - # When lambda is greater than the index we will have to go to end of the buffer to get the remaining - # terms in the sum and add it to the terms summed from index 1 to current index. Otherwise we can - # just start the sum lambda terms in front of the index - startp1 = diff < 0 ? 1 : (diff + 1) - # Since the index is less than the width, we must add the remianing terms to the sum summing the - # terms starting at the end of the buffer-remaining terms to the length of the buffer. If this is not - # necessary we skip the second loop - startp2 = diff < 0 ? ma_info.lambda2 + diff + 1 : 2 - endp2 = diff < 0 ? ma_info.lambda2 : 1 - # Compute the moving average two loop setup required when lambda < lambda2 - for i in startp1:ma_info.idx - accum += ma_info.res_window[i] - accum2 += ma_info.res_window[i]^2 - end - - for i in startp2:endp2 - accum += ma_info.res_window[i] - accum2 += ma_info.res_window[i]^2 - end - - #Update the log variable with the information for this update - if mod(iter, log.collection_rate) == 0 || iter == 0 - push!(log.lambda_hist, ma_info.lambda) - push!(log.resid_hist, accum / ma_info.lambda) - push!(log.iota_hist, accum2 / ma_info.lambda) - end - else - # Get the difference between the start and current lambda + # Consider the case when lambda <= lambda1 or lambda1 < lambda < lambda2 diff = ma_info.idx - ma_info.lambda - - # Determine start point for first loop + # Because the storage of the residual is based dependent on lambda2 and + # we want to sum only the previous lamdda terms we could have a situation + # where we want the first `idx` terms of the buffer and the last `diff` + # terms of the buffer. Doing this requires two loops + # If `diff` is negative there idx is not far enough into the buffer and + # two sums will be needed startp1 = diff < 0 ? 1 : (diff + 1) - # Determine start and endpoints for second loop - startp2 = diff > 0 ? 2 : lambda_base + diff + 1 - endp2 = diff > 0 ? 1 : lambda_base + # Assuming that the width of the buffer is lambda2 + startp2 = diff < 0 ? ma_info.lambda2 + diff + 1 : 2 + endp2 = diff < 0 ? ma_info.lambda2 : 1 # Compute the moving average two loop setup required when lambda < lambda2 for i in startp1:ma_info.idx @@ -306,7 +278,7 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: push!(log.iota_hist, accum2 / ma_info.lambda) end - ma_info.lambda += 1 + ma_info.lambda += ma_info.lambda < lambda_base ? 1 : 0 end end @@ -315,7 +287,8 @@ end """ get_uncertainty(log::LSLogMA; alpha = 0.05) -A function that takes a LSLogMA type and a confidence level, `alpha`, and returns a `(1-alpha)`-credible intervals for for every rho in the log, specifically it returns a tuple with (rho, Upper bound, Lower bound). +A function that takes a LSLogMA type and a confidence level, `alpha`, and returns a `(1-alpha)`-credible intervals +for every `rho` in the `log`, specifically it returns a tuple with (rho, Upper bound, Lower bound). """ function get_uncertainty(hist::LSLogMA; alpha::Float64 = 0.05) lambda = hist.ma_info.lambda @@ -343,9 +316,9 @@ function get_uncertainty(hist::LSLogMA; alpha::Float64 = 0.05) #compute error bound when there is an omega diffG = sqrt(cG * 2 * log(2/(alpha))) diffO = sqrt(iota) * 2 * log(2/(alpha)) * hist.dist_info.omega / (hist.dist_info.eta * width) - diffM = min(diffG, diffO) - upper[i] = rho + diffG - lower[i] = rho - diffG + diffM = max(diffG, diffO) + upper[i] = rho + diffM + lower[i] = rho - diffM end end @@ -360,7 +333,7 @@ A function that returns a default set of sub-Exponential constants for each samp This function is not exported and thus the user does not have direct access to it. # Inputs -- `log::LSLogMA`, the log containing al the tracking information. +- `log::LSLogMA`, the log containing all the tracking information. - `sampler::Type{LinSysSampler}`, the type of sampler being used. # Outputs diff --git a/test/linear_solver_logs/proc_solve_log_ma.jl b/test/linear_solver_logs/proc_solve_log_ma.jl index 8e56ef3..f031bb6 100644 --- a/test/linear_solver_logs/proc_solve_log_ma.jl +++ b/test/linear_solver_logs/proc_solve_log_ma.jl @@ -39,8 +39,9 @@ Random.seed!(1010) struct MadeUpSampler <: LinSysSampler end - @test_logs (:warn, "No constants defined for method of type Main.ProceduralTestLSLogMA.MadeUpSampler. By default we set sigma2 to 1 and scaling to 1.") RLinearAlgebra.get_SE_constants!(logger, MadeUpSampler) + @test_logs (:warn, "No constants defined for method of type Main.ProceduralTestLSLogMA.MadeUpSampler. By default we set sigma2 to 1 and scaling to 1.") RLinearAlgebra.get_SE_constants!(logger, MadeUpSampler) + @test_throws ArgumentError("`sampler` is not of type `LinSysBlkColSampler`, `LinSysVecColSampler`, `LinSysBlkRowSampler`, or `LinSysVecRowSampler`") RLinearAlgebra.log_update!(logger, MadeUpSampler(), z, (A[1, :], b[1]), 0, A, b) end # Verify late moving average