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 a66a6b8..6a96513 100644 --- a/src/linear_solver_logs/solve_log_ma.jl +++ b/src/linear_solver_logs/solve_log_ma.jl @@ -8,22 +8,22 @@ 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 - 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. +- `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. 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 @@ -37,23 +37,30 @@ 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`. +- `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 `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 `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 + 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} @@ -62,39 +69,33 @@ mutable struct SEDistInfo sigma2::Union{Float64, Nothing} omega::Union{Float64, Nothing} eta::Float64 + scaling::Float64 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. -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`). - `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. -- `iterations::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. +- `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 `false`. - `dist_info::SEDistInfo`, [`SEDistInfo`](@ref) # Constructors @@ -106,7 +107,14 @@ Solvers through Uncertainty Quantification.” SIAM/ASA J. Uncertain. Quantifica sigma2 = nothing, 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 @@ -115,9 +123,8 @@ 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 end @@ -128,17 +135,15 @@ LSLogMA(; sigma2 = nothing, omega = nothing, eta = 1, - true_res = false ) = LSLogMA( collection_rate, - MAInfo(lambda1, lambda2, 1, false, 1, zeros(lambda2)), + MAInfo(lambda1, lambda2, lambda1, false, 1, zeros(lambda2)), Float64[], Float64[], Int64[], norm, -1, false, - true_res, - SEDistInfo(nothing, 0, 0, sigma2, omega, eta) + SEDistInfo(nothing, 0, 0, sigma2, omega, eta, 0) ) #Function to update the moving average @@ -149,56 +154,50 @@ 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 + 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 - 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 - get_SE_constants!(log, log.dist_info.sampler) + # 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, typeof(sampler)) end - + end ma_info = log.ma_info - log.iterations = 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 = 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 + log.iteration = iter + # 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 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] - update_ma!(log, res, ma_info.lambda1, iter) - else - update_ma!(log, res, ma_info.lambda1, iter) - ma_info.flag = true - end + # 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 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 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,13 +224,15 @@ 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 ? 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 == 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] @@ -246,15 +246,19 @@ function update_ma!(log::LSLogMA, res::Union{AbstractVector, Real}, lambda_base: 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 @@ -274,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 @@ -283,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 @@ -292,7 +297,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 @@ -311,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 @@ -328,14 +333,17 @@ 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 -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. 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, @@ -347,6 +355,7 @@ for type in (LinSysVecRowDetermCyclic,LinSysVecRowHopRandCyclic, @eval begin function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) 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 @@ -359,6 +368,7 @@ for type in (LinSysVecColOneRandCyclic, LinSysVecColDetermCyclic) @eval begin function get_SE_constants!(log::LSLogMA, sampler::Type{$type}) 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 @@ -371,6 +381,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/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 9f70a49..f031bb6 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 @@ -23,19 +23,28 @@ Random.seed!(1010) z = rand(2) sampler = LinSysVecRowOneRandCyclic() - log = LSLogMA() + 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!(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 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) + @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 moving average + # Verify late moving average let A = rand(2,2) x = rand(2) @@ -43,33 +52,70 @@ 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(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(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.iteration == 10 + @test logger.converged == false + + #Test uncertainty set + 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] - 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 + A = rand(2,2) + x = rand(2) + b = A * x + z = rand(2) + + sampler = LinSysVecRowOneRandCyclic() + logger = LSLogMA(lambda1 = 2, + lambda2 = 10) + samp = (A[1,:], b[1]) + + 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!(logger, sampler, x + .3^(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] - @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 log.iterations == 10 - @test log.converged == false + 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(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.iteration == 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 + + + end 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