From 927330358b573459dbf3a9b6288209a3a99487ea Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Sun, 5 Oct 2025 23:06:27 -0700 Subject: [PATCH 1/2] QR refactoring --- src/BPDual.jl | 28 ++++++++++++++++++++++------ src/helpers.jl | 2 +- src/omp.jl | 10 +++++++--- test/test_bpdn.jl | 5 +++-- test/test_omp.jl | 4 ++-- 5 files changed, 35 insertions(+), 14 deletions(-) diff --git a/src/BPDual.jl b/src/BPDual.jl index ec7452b..4816266 100644 --- a/src/BPDual.jl +++ b/src/BPDual.jl @@ -11,7 +11,15 @@ function Base.getindex(t::ASPTracer, i::Integer) end Base.lastindex(t::ASPTracer) = lastindex(t.iteration) +function Base.show(io::IO, t::ASPTracer) + nsteps = length(t.iteration) + nvars = isempty(t.solution) ? 0 : length(t.solution[end].nzind) + println(io, "ASPTracer with $nvars active variables at final step.") +end +function Base.show(io::IO, ::MIME"text/plain", t::ASPTracer) + Base.show(io, t) +end @doc raw""" ```julia function bpdual(A, b, λin, bl, bu; kwargs...)``` @@ -89,7 +97,8 @@ function bpdual( gapTol::Real = 1e-06, pivTol::Real = 1e-12, actMax::Real = Inf, - traceFlag::Bool = false + traceFlag::Bool = false, + refactor_freq::Int = 1000 ) # Start @@ -282,9 +291,8 @@ function bpdual( @info "\nOptimal solution found. Trimming multipliers..." end g = b - λin*y - trimx(x, (@view S[:, 1:cur_r_size]), (@view R[1:cur_r_size, 1:cur_r_size]), active, state, g, b, λ, feaTol, optTol, loglevel) - numtrim = nact - length(active) - nact = length(active) + x, active = trimx(x, (@view S[:,1:cur_r_size]), (@view R[1:cur_r_size, 1:cur_r_size]), active, state, g, b, λ, feaTol, optTol, loglevel) + break end # Act on any live exit conditions. @@ -349,7 +357,11 @@ function bpdual( nprodA += 1 zerovec[p] = 0 qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5) # Update R - cur_r_size +=1 + cur_r_size +=1 + if itn % refactor_freq == 0 && cur_r_size > 0 + F = qr!(S[:, 1:cur_r_size]) + @views R[1:cur_r_size, 1:cur_r_size] = F.R + end # S = [S a] push!(active, p) push!(x, 0) @@ -373,7 +385,11 @@ function bpdual( deleteat!(x, qa) # R = qrdelcol(R, qa) qrdelcol!(S, R, qa) - cur_r_size -=1 + cur_r_size -=1 + if itn % refactor_freq == 0 && cur_r_size > 0 + F = qr!(S[:, 1:cur_r_size]) + @views R[1:cur_r_size, 1:cur_r_size] = F.R + end else eFlag = :EXIT_OPTIMAL end diff --git a/src/helpers.jl b/src/helpers.jl index 84c6218..0826609 100644 --- a/src/helpers.jl +++ b/src/helpers.jl @@ -192,7 +192,7 @@ function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel) # Grab the next canddate multiplier. xmin, qa = findmin(xabs) end - return (x,S,R,active,state) + return x, active end diff --git a/src/omp.jl b/src/omp.jl index 63ba825..0216656 100644 --- a/src/omp.jl +++ b/src/omp.jl @@ -32,7 +32,8 @@ function asp_omp( gapTol::Real = 1e-06, pivTol::Real = 1e-12, actMax::Union{Real, Nothing} = nothing, - traceFlag::Bool = false) + traceFlag::Bool = false, + refactor_freq::Int = 1000) time0 = time() @@ -198,7 +199,10 @@ function asp_omp( qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5) # Update R # S = hcat(S, a) # Expand S, active cur_r_size +=1 - + if itn % refactor_freq == 0 && cur_r_size > 0 + F = qr!(S[:, 1:cur_r_size]) + @views R[1:cur_r_size, 1:cur_r_size] = F.R + end if traceFlag push!(tracer.iteration, itn) push!(tracer.lambda, zmax) @@ -218,7 +222,7 @@ function asp_omp( tottime = time() - time0 if loglevel > 0 - @info @sprintf("\nEXIT BPdual -- %s\n", EXIT_INFO[eFlag]) + @info @sprintf("\nEXIT OMP -- %s\n", EXIT_INFO[eFlag]) @info @sprintf("%-20s: %8i", "Products with A", nprodA) @info @sprintf("%-20s: %8i", "Products with At", nprodAt) @info @sprintf("%-20s: %8.1e", "Solution time (sec)", tottime) diff --git a/test/test_bpdn.jl b/test/test_bpdn.jl index 106b597..decb723 100644 --- a/test/test_bpdn.jl +++ b/test/test_bpdn.jl @@ -8,7 +8,7 @@ using Test, LinearAlgebra, Random, SparseArrays, ActiveSetPursuit, LinearOperato function test_bpdn() m = 600 n = 2560 - k = 20 + k = 100 # Generate sparse solution p = randperm(n)[1:k] # Position of nonzeros in x @@ -21,7 +21,8 @@ function test_bpdn() b = A * x # Solve the basis pursuit problem - tracer = asp_bpdn(A, b, 0.0, loglevel =0); + tracer = asp_bpdn(A, b, 0.0, loglevel =0, refactor_freq =20); + xx, λ = tracer[end] pFeas = norm(A * xx - b, Inf) / max(1, norm(xx, Inf)) diff --git a/test/test_omp.jl b/test/test_omp.jl index d581c1c..d05278a 100644 --- a/test/test_omp.jl +++ b/test/test_omp.jl @@ -8,7 +8,7 @@ using Test, LinearAlgebra, Random, SparseArrays, ActiveSetPursuit, LinearOperato function test_omp() m = 600 n = 2560 - k = 20 + k = 100 # Generate sparse solution p = randperm(n)[1:k] # Position of nonzeros in x @@ -21,7 +21,7 @@ function test_omp() b = A * x # Solve the basis pursuit problem - tracer = asp_omp(A, b, 0.0; loglevel =0); + tracer = asp_omp(A, b, 0.0; loglevel =0, refactor_freq =20); xx, _ = tracer[end] pFeas = norm(A * xx - b, Inf) / max(1, norm(xx, Inf)) From 1b94cf62637cb03aad860b73e75053861a55b5d4 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Mon, 6 Oct 2025 12:48:17 -0700 Subject: [PATCH 2/2] added exit singular --- src/omp.jl | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/omp.jl b/src/omp.jl index 0216656..bec707c 100644 --- a/src/omp.jl +++ b/src/omp.jl @@ -134,6 +134,20 @@ function asp_omp( # Main loop. while true # Compute dual obj gradient g, search direction dy, and residual r. + + + if isempty((@view R[1:cur_r_size,1:cur_r_size])) + condS = 1 + else + rmin = minimum(diag((@view R[1:cur_r_size, 1:cur_r_size]))) + rmax = maximum(diag((@view R[1:cur_r_size, 1:cur_r_size]))) + condS = rmax / rmin + end + + if condS > 1e+10 + eFlag = :EXIT_SINGULAR_LS + end + if itn == 0 x = Float64[] r = b