diff --git a/src/BPDual.jl b/src/BPDual.jl index ec7452b..85e6484 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) + print(io, "ASPTracer(maxactive = $nvars, nsteps = $nsteps)") +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 @@ -101,11 +110,11 @@ function bpdual( m, n = size(A) - work = rand(n) - work2 = rand(n) - work3 = rand(n) - work4 = rand(n) - work5 = rand(m) + work = Vector{Float64}(undef, n) + work2 = Vector{Float64}(undef, n) + work3 = Vector{Float64}(undef, n) + work4 = Vector{Float64}(undef, n) + work5 = Vector{Float64}(undef, m) tracer = ASPTracer( Int[], # iteration @@ -178,14 +187,14 @@ function bpdual( numtrim = 0 nprodA = 0 nprodAt = 0 - cur_r_size = 0 - # ------------------------------------------------------------------ # Cold/warm-start initialization. # ------------------------------------------------------------------ if coldstart x = zeros(0) + # cur_r_size = 0 + nact = 0 if homotopy y = b / λ z = z / λ @@ -201,6 +210,7 @@ function bpdual( y = restorefeas(y, active, state, S, R, bl, bu) z = A'*y nprodAt += 1 + # cur_r_size = nact end sL, sU = infeasibilities(bl, bu, vec(z)) @@ -220,26 +230,32 @@ function bpdual( sL, sU = infeasibilities(bl, bu, z) g = b - λ*y # Steepest-descent direction - if isempty((@view R[1:cur_r_size,1:cur_r_size])) + # if isempty((@view R[1:cur_r_size,1:cur_r_size])) + if isempty((@view R[1:nact,1:nact])) 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]))) + # 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]))) + rmin = minimum(diag((@view R[1:nact, 1:nact]))) + rmax = maximum(diag((@view R[1:nact, 1:nact]))) 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((@view S[:, 1:cur_r_size]), 2) - size(x, 1) + # npad = size((@view S[:, 1:cur_r_size]), 2) - size(x, 1) + npad = size((@view S[:, 1:nact]), 2) - size(x, 1) x = [x; zeros(npad)] else - dx, dy = newtonstep((@view S[:,1:cur_r_size]), (@view R[1:cur_r_size, 1:cur_r_size]), g, x, λ) + dx, dy = newtonstep((@view S[:,1:nact]), (@view R[1:nact, 1:nact]), 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[:,1:cur_r_size]*x - + # r = b - S[:,1:cur_r_size]*x + r = b - S[:,1:nact]*x + # Print to log. yNorm = norm(y, 2) rNorm = norm(r, 2) @@ -282,9 +298,9 @@ 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:nact]), (@view R[1:nact, 1:nact]), active, state, g, b, λ, feaTol, optTol, loglevel) + # 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. @@ -297,7 +313,8 @@ function bpdual( p = q = 0 if homotopy - 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) + # 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) + x, dy, dz, step, λ, p = htpynewlam(active, state, A, (@view R[1:nact, 1:nact]), (@view S[:,1:nact]), x, y, sL, sU, λ, lamFinal) nprodAt += 1 else if norm(dy, Inf) < eps() @@ -348,8 +365,20 @@ function bpdual( a = A * zerovec nprodA += 1 zerovec[p] = 0 - qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5) # Update R - cur_r_size +=1 + # qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5) # Update R + + qraddcol!(S, R, a, nact, work, work2, work3, work4, work5) # Update R + + # cur_r_size +=1 + nact +=1 + if itn % refactor_freq == 0 && nact > 0 + F = qr!(S[:, 1:nact]) + @views R[1:nact, 1:nact] = F.R + end + # 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 +402,16 @@ function bpdual( deleteat!(x, qa) # R = qrdelcol(R, qa) qrdelcol!(S, R, qa) - cur_r_size -=1 + # cur_r_size -=1 + nact -=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 itn % refactor_freq == 0 && nact > 0 + F = qr!(S[:, 1:nact]) + @views R[1:nact, 1:nact] = 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..3662eba 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() @@ -55,7 +56,7 @@ function asp_omp( tracer = ASPTracer( Int[], # iteration Float64[], # lambda - Vector{SparseVector{Float64}}() # now stores full sparse solutions + Vector{SparseVector{Float64}}() ) if loglevel > 0 @@ -66,8 +67,6 @@ function asp_omp( @info "-"^124 end - - # Initialize local variables. EXIT_INFO = Dict( :EXIT_OPTIMAL => "Optimal solution found -- full Newton step", :EXIT_TOO_MANY_ITNS => "Too many iterations", @@ -84,14 +83,14 @@ function asp_omp( x = zeros(Float64, 0) zerovec = zeros(Float64, n) p = 0 - cur_r_size = 0 + # cur_r_size = 0 + int_ac = 0 if norm(b, Inf) == 0 r = zeros(m) eFlag = :EXIT_RHS_ZERO end - # Solution is unconstrained for lambda large. zmax = norm(z, Inf) if eFlag == :EXIT_UNKNOWN && zmax < λin r = b @@ -105,12 +104,13 @@ function asp_omp( 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] + # cur_r_size = length(active) + # S[:, 1:cur_r_size] .= A[:, active] + S[:, 1:int_ac] .= 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 + # println(size(R_)) + # R[1:cur_r_size, 1:cur_r_size] .= R_ + R[1:int_ac, 1:int_ac] .= R_ itn = 1 end @@ -133,6 +133,17 @@ function asp_omp( # Main loop. while true # Compute dual obj gradient g, search direction dy, and residual r. + if isempty((@view R[1:int_ac,1:int_ac])) + # 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]))) + rmin = minimum(diag((@view R[1:int_ac, 1:int_ac]))) + rmax = maximum(diag((@view R[1:int_ac, 1:int_ac]))) + condS = rmax / rmin + end + if itn == 0 x = Float64[] r = b @@ -140,26 +151,40 @@ function asp_omp( nprodAt += 1 zmax = norm(z, Inf) else - x,y = csne( (@view R[1:cur_r_size, 1:cur_r_size]), - (@view S[:,1:cur_r_size]), vec(b)) + # x,y = csne( (@view R[1:cur_r_size, 1:cur_r_size]), + # (@view S[:,1:cur_r_size]), vec(b)) + x,y = csne( (@view R[1:int_ac, 1:int_ac]), + (@view S[:,1:int_ac]), vec(b)) if norm(x, Inf) > 1e12 eFlag = :EXIT_SINGULAR_LS break end - Sx = (@view S[:,1:cur_r_size]) * x + Sx = (@view S[:,1:int_ac]) * x + # Sx = (@view S[:,1:cur_r_size]) * x + r = b - Sx end rNorm = norm(r, 2) xNorm = norm(x, 1) + if traceFlag + push!(tracer.iteration, itn) + push!(tracer.lambda, zmax) + sparse_x_full = spzeros(n) + @views act_now = (itn == 0) ? Int[] : active[1:int_ac] + # @views act_now = (itn == 0) ? Int[] : active[1:cur_r_size] + @assert length(act_now) == length(x) + sparse_x_full[act_now] = x + push!(tracer.solution, copy(sparse_x_full)) + end + 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 - # Already set. Don't test the other exits. elseif zmax <= λin eFlag = :EXIT_LAMBDA elseif rNorm <= optTol @@ -173,41 +198,52 @@ function asp_omp( if eFlag != :EXIT_UNKNOWN break end - - # New iteration starts here. + + # New iter itn += 1 # Find step to the nearest inactive constraint z = A_T * r - # mul!(z, A', r) nprodAt += 1 - zmax, p = findmax(abs.(z)) - - # if z[p] < 0 - # state[p] = -1 - # else - # state[p] = 1 - # end + z_masked = copy(z) + z_masked[active] .= 0 + zmax, p = findmax(abs.(z_masked)) + if p in active + continue + end + + push!(active, p) - zerovec[p] = 1 # Extract a = A[:, p] - a = A * zerovec # Compute A[:, p] + zerovec[p] = 1 # get a = A[:, p] + a = A * zerovec nprodA += 1 zerovec[p] = 0 + qraddcol!(S, R, a, int_ac, work, work2, work3, work4, work5) # Update R - qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5) # Update R + # 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 + # cur_r_size +=1 + int_ac+=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) - sparse_x_full = spzeros(n) - sparse_x_full[copy(active)] = copy(x) - push!(tracer.solution, copy(sparse_x_full)) + if itn % refactor_freq == 0 && int_ac > 0 + F = qr!(S[:, 1:int_ac]) + @views R[1:int_ac, 1:int_ac] = F.R end - push!(active, p) + if condS > 1e+6 + F = qr!(S[:, 1:int_ac]) + @views R[1:int_ac, 1:int_ac] = F.R + end + + # if condS > 1e+6 + # F = qr!(S[:, 1:cur_r_size]) + # @views R[1:cur_r_size, 1:cur_r_size] = F.R + # end end #while true push!(tracer.iteration, itn) @@ -218,7 +254,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))