From 713d48e3698ef2efd732d7b37b62207ae4ac8d19 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Tue, 10 Jun 2025 15:30:32 -0700 Subject: [PATCH] update --- src/BPDual.jl | 44 ++++++++------ src/helpers.jl | 52 ++++++++++------- src/omp.jl | 134 +++++++++++++++++++++++++++---------------- test/runtests.jl | 2 + test/test_omp.jl | 49 ++++++++++++++++ test/test_triminf.jl | 39 +++++++++++++ 6 files changed, 234 insertions(+), 86 deletions(-) create mode 100644 test/test_omp.jl create mode 100644 test/test_triminf.jl diff --git a/src/BPDual.jl b/src/BPDual.jl index e8690f8..ec7452b 100644 --- a/src/BPDual.jl +++ b/src/BPDual.jl @@ -100,6 +100,13 @@ function bpdual( # ------------------------------------------------------------------ m, n = size(A) + + work = rand(n) + work2 = rand(n) + work3 = rand(n) + work4 = rand(n) + work5 = rand(m) + tracer = ASPTracer( Int[], # iteration Float64[], # lambda @@ -110,8 +117,8 @@ function bpdual( active = Vector{Int}([]) state = Vector{Int}(zeros(Int, n)) y = Vector{Float64}(zeros(Float64, m)) - S = Matrix{Float64}(zeros(Float64, m, 0)) - R = Matrix{Float64}(zeros(Float64, 0, 0)) + S = Matrix{Float64}(zeros(Float64, m, n)) + R = Matrix{Float64}(zeros(Float64, n, n)) end if homotopy @@ -171,6 +178,8 @@ function bpdual( numtrim = 0 nprodA = 0 nprodAt = 0 + cur_r_size = 0 + # ------------------------------------------------------------------ # Cold/warm-start initialization. @@ -186,7 +195,7 @@ function bpdual( end else g = b - λ*y # Compute steepest-descent dir - triminf(active, state, S, R, bl, bu, g) + active, state, S, R = triminf!(active, state, S, R, bl, bu, g) nact = length(active) x = zeros(nact) y = restorefeas(y, active, state, S, R, bl, bu) @@ -211,25 +220,25 @@ function bpdual( sL, sU = infeasibilities(bl, bu, z) g = b - λ*y # Steepest-descent direction - if isempty(R) + if isempty((@view R[1:cur_r_size,1:cur_r_size])) condS = 1 else - rmin = minimum(diag(R)) - rmax = maximum(diag(R)) + 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 # Pad x with enough zeros to make it compatible with S. - npad = size(S, 2) - size(x, 1) + npad = size((@view S[:, 1:cur_r_size]), 2) - size(x, 1) x = [x; zeros(npad)] else - dx, dy = newtonstep(S, R, g, x, λ) + dx, dy = newtonstep((@view S[:,1:cur_r_size]), (@view R[1:cur_r_size, 1:cur_r_size]), g, x, λ) x .+= dx end - r = b - S*x + r = b - S[:,1:cur_r_size]*x # Print to log. yNorm = norm(y, 2) @@ -273,7 +282,7 @@ function bpdual( @info "\nOptimal solution found. Trimming multipliers..." end g = b - λin*y - trimx(x, S, R, active, state, g, b, λ, feaTol, optTol, loglevel) + 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) end @@ -288,7 +297,7 @@ function bpdual( p = q = 0 if homotopy - x, dy, dz, step, λ, p = htpynewlam(active, state, A, R, S, x, y, sL, sU, λ, lamFinal) + x, dy, dz, step, λ, p = htpynewlam(active, state, A, (@view R[1:cur_r_size, 1:cur_r_size]), (@view S[:,1:cur_r_size]), x, y, sL, sU, λ, lamFinal) nprodAt += 1 else if norm(dy, Inf) < eps() @@ -339,8 +348,9 @@ function bpdual( a = A * zerovec nprodA += 1 zerovec[p] = 0 - R = qraddcol(S, R, a) - S = [S a] + qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5) # Update R + cur_r_size +=1 + # S = [S a] push!(active, p) push!(x, 0) else @@ -358,10 +368,12 @@ function bpdual( _, qa = findmax(abs.(x .* dropa)) q = active[qa] state[q] = 0 - S = S[:, 1:nact .!= qa] + # S = S[:, 1:nact .!= qa] deleteat!(active, qa) deleteat!(x, qa) - R = qrdelcol(R, qa) + # R = qrdelcol(R, qa) + qrdelcol!(S, R, qa) + cur_r_size -=1 else eFlag = :EXIT_OPTIMAL end @@ -390,4 +402,4 @@ function bpdual( @info "Solution time (sec): $tottime" end return tracer -end +end \ No newline at end of file diff --git a/src/helpers.jl b/src/helpers.jl index f13df18..84c6218 100644 --- a/src/helpers.jl +++ b/src/helpers.jl @@ -201,37 +201,49 @@ end # ---------------------------------------------------------------------- # Trim working constraints with "infinite" bounds. # ---------------------------------------------------------------------- - -function triminf(active::Vector, state::Vector, S::Matrix, R::Matrix, - bl::Vector, bu::Vector, g::Vector, b::Vector) - +function triminf!( + active::Vector{Int}, + state::Vector{Int}, + S::Matrix{Float64}, + R, + bl::Vector{Float64}, + bu::Vector{Float64}, +) bigbnd = 1e10 - nact = length(active) - tlistbl = find( state[active] .== -1 & bl[active] .< -bigbnd ) - tlistbu = find( state[active] .== +1 & bu[active] .> +bigbnd ) - tlist = [tlistbl; tlistbu] + tlist = [i for i in active if + (state[i] == -1 && bl[i] < -bigbnd) || + (state[i] == 1 && bu[i] > bigbnd) + ] if isempty(tlist) - return + return active, state, S, R end - for q in tlist - qa = active[q] - nact = nact - 1 - S = S[:,1:size(S,2) .!= qa] # Delete column from S - deleteat!(active, qa) # Delete index from active set - R = qrdelcol(R, qa) # Recompute new QR factorization - state[q] = 0 # Mark constraint as free - x = csne(R, S, g) # Recompute multipliers + for q in sort(tlist; rev=true) + qa = active[q] - rNorm = norm(b - S*x) - xNorm = norm(x, 1) + if qa == 1 + S = S[:, 2:end] + elseif qa == size(S, 2) + S = S[:, 1:end-1] + else + S = hcat(S[:, 1:qa-1], S[:, qa+1:end]) end + deleteat!(active, q) + + R = qrdelcol(R, qa) + + state[q] = 0 + end + + return active, state, S, R end + + ## G) find_step # ---------------------------------------------------------------------- @@ -387,4 +399,4 @@ function htpynewlam(active, state, A, R, S, x, y, s1, s2, λ, lamFinal) end return x, dy, dz, step, λ, p -end +end \ No newline at end of file diff --git a/src/omp.jl b/src/omp.jl index 14a7f32..bc6f623 100644 --- a/src/omp.jl +++ b/src/omp.jl @@ -1,16 +1,17 @@ -struct OMPTracer{T} - iteration::Vector{Int} - lambda::Vector{T} - solution::Vector{SparseVector{T}} -end +# struct OMPTracer{T} +# iteration::Vector{Int} +# lambda::Vector{T} +# active::Vector{Vector{Int}} # Store active indices instead of full sparse vectors +# values::Vector{Vector{T}} # Store corresponding values for active indices +# end -Base.length(t::OMPTracer) = length(t.iteration) +# Base.length(t::OMPTracer) = length(t.iteration) -function Base.getindex(t::OMPTracer, i::Integer) - return t.solution[i], t.lambda[i] -end +# function Base.getindex(t::OMPTracer, i::Integer) +# return t.active[i], t.values[i], t.lambda[i] +# end -Base.lastindex(t::OMPTracer) = lastindex(t.iteration) +# Base.lastindex(t::OMPTracer) = lastindex(t.iteration) @doc raw""" ```julia @@ -36,7 +37,6 @@ function asp_omp( b::Vector, λin::Real; active::Union{Nothing, Vector{Int}} = nothing, - state::Union{Nothing, Vector{Int}} = nothing, S::Matrix{Float64} = Matrix{Float64}(undef, size(A, 1), 0), R::Union{Nothing, Matrix{Float64}} = nothing, loglevel::Int = 1, @@ -46,28 +46,33 @@ function asp_omp( optTol::Real = 1e-05, gapTol::Real = 1e-06, pivTol::Real = 1e-12, - actMax::Real = Inf) + actMax::Union{Real, Nothing} = nothing, + traceFlag::Bool = false) - # Start the clock and size up the problem. time0 = time() - z = A' * b - + A_T = A' + z = A_T * b + int_ac =0 m = length(b) n = length(z) + T = eltype(A) + + work = Vector{T}(undef, size(A, 2)) + work2 = Vector{T}(undef, size(A, 2)) + work3 = Vector{T}(undef, size(A, 2)) + work4 = Vector{T}(undef, size(A, 2)) + work5 = Vector{T}(undef, size(A, 1)) + nprodA = 0 nprodAt = 1 - # Initialize the tracer - - tracer = OMPTracer( + tracer = ASPTracer( Int[], # iteration Float64[], # lambda Vector{SparseVector{Float64}}() # now stores full sparse solutions ) - # Print log header. - if loglevel > 0 @info "-"^124 @info @sprintf("%-30s : %-10d %-30s : %-10.4e", "No. rows", m, "λ", λin) @@ -76,6 +81,7 @@ function asp_omp( @info "-"^124 end + # Initialize local variables. EXIT_INFO = Dict( :EXIT_OPTIMAL => "Optimal solution found -- full Newton step", @@ -84,6 +90,7 @@ function asp_omp( :EXIT_LAMBDA => "Reached minimum value of lambda", :EXIT_RHS_ZERO => "b = 0. The solution is x = 0", :EXIT_UNCONSTRAINED => "Unconstrained solution r = b is optimal", + :EXIT_ACTMAX => "Max no. of active constraints reached", :EXIT_UNKNOWN => "unknown exit" ) @@ -92,8 +99,8 @@ function asp_omp( x = zeros(Float64, 0) zerovec = zeros(Float64, n) p = 0 + cur_r_size = 0 - # Quick exit if the RHS is zero. if norm(b, Inf) == 0 r = zeros(m) eFlag = :EXIT_RHS_ZERO @@ -106,18 +113,37 @@ function asp_omp( eFlag = :EXIT_UNCONSTRAINED end - if eFlag != :EXIT_UNKNOWN || active === nothing + if active === nothing active = Vector{Int}([]) + else + active = active + R = Matrix{Float64}(undef, size(A, 2), size(A, 2)) + S = Matrix{Float64}(undef, size(A, 1), size(A, 2)) + int_ac= length(active) + cur_r_size = length(active) + S[:, 1:cur_r_size] .= A[:, active] + _, R_ = qr(A[:, active]) + # Use regular indexing for R + println(size(R_)) + R[1:cur_r_size, 1:cur_r_size] .= R_ # Copy values into R + itn = 1 + end - if state === nothing - state = zeros(Int, n) - end + # if state === nothing + # state = zeros(Int, n) + # end if R === nothing - R = Matrix{Float64}(undef, 0, 0) + R = Matrix{Float64}(undef, size(A,2), size(A,2)) + S = Matrix{Float64}(undef, size(A,1), size(A,2)) end - @info @sprintf("%4s %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm") + if actMax === nothing + actMax = size(A, 2) + end + if loglevel>0 + @info @sprintf("%4s %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm") + end # Main loop. while true @@ -125,24 +151,26 @@ function asp_omp( if itn == 0 x = Float64[] r = b - z = A' * r + z = A_T * r nprodAt += 1 zmax = norm(z, Inf) else - x,y = csne(R, S, vec(b)) + x,y = csne( (@view R[1:cur_r_size, 1:cur_r_size]), + (@view S[:,1:cur_r_size]), vec(b)) if norm(x, Inf) > 1e12 eFlag = :EXIT_SINGULAR_LS break end - - Sx = S * x + Sx = (@view S[:,1:cur_r_size]) * x r = b - Sx end rNorm = norm(r, 2) xNorm = norm(x, 1) - @info @sprintf("%4i %8i %12.5e %12.5e %12.5e", itn, p, zmax, rNorm, xNorm) + if loglevel>0 + @info @sprintf("%4i %8i %12.5e %12.5e %12.5e", itn, p, zmax, rNorm, xNorm) + end # Check exit conditions. if eFlag != :EXIT_UNKNOWN @@ -153,8 +181,10 @@ function asp_omp( eFlag = :EXIT_OPTIMAL elseif itn >= itnMax eFlag = :EXIT_TOO_MANY_ITNS + elseif itn-1 == actMax- int_ac + eFlag = :EXIT_ACTMAX end - + if eFlag != :EXIT_UNKNOWN break end @@ -163,38 +193,42 @@ function asp_omp( itn += 1 # Find step to the nearest inactive constraint - z = A' * r - + z = A_T * r + # mul!(z, A', r) nprodAt += 1 zmax, p = findmax(abs.(z)) + push!(active, p) - if z[p] < 0 - state[p] = -1 - else - state[p] = 1 - end + # if z[p] < 0 + # state[p] = -1 + # else + # state[p] = 1 + # end zerovec[p] = 1 # Extract a = A[:, p] - a = A * zerovec + a = A * zerovec # Compute A[:, p] nprodA += 1 zerovec[p] = 0 - R = qraddcol(S, R, a) # Update R - S = hcat(S, a) # Expand S, active + 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 - push!(tracer.iteration, itn) - push!(tracer.lambda, zmax) - sparse_x_full = SparseVector(n, copy(active), copy(x)) - push!(tracer.solution, copy(sparse_x_full)) - - push!(active, p) + if traceFlag + push!(tracer.iteration, itn) + push!(tracer.lambda, zmax) + sparse_x_full = spzeros(n) + sparse_x_full[copy(active)] = copy(x) + push!(tracer.solution, copy(sparse_x_full)) + end end #while true push!(tracer.iteration, itn) push!(tracer.lambda, zmax) - sparse_x_full = SparseVector(n, copy(active), copy(x)) + sparse_x_full = spzeros(n) + sparse_x_full[copy(active)] = copy(x) push!(tracer.solution, copy(sparse_x_full)) tottime = time() - time0 diff --git a/test/runtests.jl b/test/runtests.jl index 963436b..8864d38 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,5 +10,7 @@ println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡ @testset "ActiveSetPursuit.jl" begin @testset "BPDN" begin include("test_bpdn.jl") end @testset "Recover decaying coefficients" begin include("test_recover_decaying.jl") end + @testset "OMP" begin include("test_omp.jl") end + @testset "triminf" begin include("test_triminf.jl") end end diff --git a/test/test_omp.jl b/test/test_omp.jl new file mode 100644 index 0000000..d581c1c --- /dev/null +++ b/test/test_omp.jl @@ -0,0 +1,49 @@ + +# ------------------------------------------------------------------ +# Test orthogonal matching pursuit +# ------------------------------------------------------------------ + +using Test, LinearAlgebra, Random, SparseArrays, ActiveSetPursuit, LinearOperators + +function test_omp() + m = 600 + n = 2560 + k = 20 + + # Generate sparse solution + p = randperm(n)[1:k] # Position of nonzeros in x + x = zeros(n) + x[p] .= randn(k) + + A = randn(m, n) + + # Compute the RHS vector + b = A * x + + # Solve the basis pursuit problem + tracer = asp_omp(A, b, 0.0; loglevel =0); + + xx, _ = tracer[end] + pFeas = norm(A * xx - b, Inf) / max(1, norm(xx, Inf)) + @test pFeas <= 1e-6 + + # ------------------------ + # Linear operator + # ------------------------ + + OP = LinearOperator(A) + b_op = OP * x + + # Solve the basis pursuit problem + tracer_op = asp_omp(OP, b_op, 0.0; loglevel =0); + + xx_op, _ = tracer_op[end] + pFeas_op = norm(OP * xx_op - b_op, Inf) / max(1, norm(xx_op, Inf)) + @test pFeas_op <= 1e-6 +end + +Random.seed!(1234) + +for ntest = 1:10 + test_omp() +end \ No newline at end of file diff --git a/test/test_triminf.jl b/test/test_triminf.jl new file mode 100644 index 0000000..686e967 --- /dev/null +++ b/test/test_triminf.jl @@ -0,0 +1,39 @@ +using Test +using LinearAlgebra +using ActiveSetPursuit: triminf! + +# Test 1: No constraints to remove +active1 = [1, 2] +state1 = [0, 0] +m, n = 3, 2 +S1 = rand(m, n) +R1 = Matrix{Float64}(I, n, n) +bl1 = [0.0, 0.0] +bu1 = [0.0, 0.0] +g = zeros(m) +b = zeros(m) + +active_out1, state_out1, S_out1, R_out1 = + triminf!(copy(active1), copy(state1), copy(S1), copy(R1), bl1, bu1) + +@test active_out1 == active1 +@test state_out1 == state1 +@test S_out1 == S1 + +# Test 2: All constraints out of bounds +active2 = [1, 2] +state2 = [-1, 1] +S2 = rand(m, n) +R2 = Matrix{Float64}(I, n, n) +bl2 = [-1e11, 0.0] +bu2 = [0.0, 1e11] + +active_out2, state_out2, S_out2, R_out2 = + triminf!(copy(active2), copy(state2), copy(S2), copy(R2), bl2, bu2) + +@test isempty(active_out2) +@test state_out2 == [0, 0] +@test size(S_out2) == (m, 0) +@test size(R_out2) == (0, 0) +@test size(R_out2,1) == size(R_out2,2) +@test istriu(R_out2)