Merge pull request #4 from MPF-Optimization-Laboratory/tt

Added new test + fixed tracer

src/ActiveSetPursuit.jl

@@ -6,14 +6,8 @@ using DataFrames
 using LinearAlgebra, SparseArrays, LinearOperators, Printf
 using QRupdate, Random
 using LinearAlgebra
-export as_topy
 
-
-export bpdual, sparsity, newtonstep, objectives, infeasibilities
-export trimx, triminf, restorefeas, htpynewlam, find_step
-export bpdual, sparsity, newtonstep, objectives, infeasibilities
-export trimx, triminf, restorefeas, htpynewlam, find_step
-export as_topy
+export bpdual, asp_homotopy
 
 include("BPDual.jl")
 include("helpers.jl")

src/BPDual.jl

@@ -1,8 +1,24 @@
+struct ASPTracer{T}
+    iteration::Vector{Int}
+    lambda::Vector{T}
+    active::Vector{Vector{Int}}
+    activesoln::Vector{Vector{T}}
+    N::Int
+end
 
+Base.length(t::ASPTracer) = length(t.iteration)
 
+function Base.getindex(t::ASPTracer, i::Integer)
+    as = t.active[i]
+    x = zeros(t.N)
+    x[as] = t.activesoln[i]
+    return x, t.lambda[i]
+end
+
+Base.lastindex(t::ASPTracer) = lastindex(t.active)
 
 function bpdual(
-    A::AbstractMatrix,
+    A::Union{AbstractMatrix,AbstractLinearOperator},
     b::Vector,
     λin::Real,
     bl::Vector,
@@ -24,6 +40,45 @@ function bpdual(
     pivTol::Real = 1e-12,
     actMax::Real = Inf)
 
+    """
+    Solve the optimization problem:
+
+    DP:       minimize_y   -b'y  +  1/2 * λ * y'y
+              subject to   bl <= A'y <= bu
+
+    using given `A`, `b`, `bl`, `bu`, and `λ`. When `bl = -e` and `bu = e = ones(n, 1)`,
+    DP is the dual Basis Pursuit problem:
+
+    BPdual:   maximize_y   b'y  -  1/2 * λ * y'y
+            subject to   ||A'y||_inf  <=  1.
+
+    # Input
+    - `A` : `m`-by-`n` explicit matrix or operator.
+    - `b` : `m`-vector.
+    - `λin` : Nonnegative scalar.
+    - `bl`, `bu` : `n`-vectors (bl lower bound, bu upper bound).
+    - `active`, `state`, `y`, `S`, `R` : May be empty or output from `BPdual` with a previous value of `λ`.
+    - `loglevel` : Logging level.
+    - `coldstart` : Boolean indicating if a cold start should be used.
+    - `homotopy` : Boolean indicating if homotopy should be used.
+    - `λmin` : Minimum value for `λ`.
+    - `trim` : Number of constraints to trim.
+    - `itnMax` : Maximum number of iterations.
+    - `feaTol` : Feasibility tolerance.
+    - `optTol` : Optimality tolerance.
+    - `gapTol` : Gap tolerance.
+    - `pivTol` : Pivot tolerance.
+    - `actMax` : Maximum number of active constraints.
+
+    # Output
+    - `tracer` : A structure to store trace information at each iteration of the optimization process.
+        It contains:
+            - `active::Vector{Int}`: The indices of the active constraints.
+            - `activesoln::Vector{T}`: The solutions corresponding to the current active set.
+            - `lambda::Vector{T}`: Lambda values.
+            - `N::Int`: The total number of variables in the solution vector.
+    """
+
     # start 
     time0 = time()
 
@@ -32,14 +87,13 @@ function bpdual(
     # ------------------------------------------------------------------
     m, n = size(A)
 
-    tracer = DataFrame(
-    Itn = Int[],
-    Active = String[],
-    X = String[],
-    Lambda = Float64[],
-    Xnorm = Float64[],
-    Rnorm = Float64[],
-    CondS = Float64[])
+    tracer = ASPTracer(
+        Int[],                  # iteration
+        Float64[],              # lambda
+        Vector{Vector{Int}}(),  # active
+        Vector{Vector{Float64}}(), # activesoln
+        n                       # N
+    )
 
 
     if coldstart || isnothing(active) || isnothing(state) || isnothing(y) || isnothing(S) || isnothing(R)
@@ -239,7 +293,9 @@ function bpdual(
         if eFlag == :EXIT_OPTIMAL || eFlag == :EXIT_SMALL_DGAP
             # Optimal trimming. λ0 may be different from λ.
             # Recompute gradient just in case.
-            println("\nOptimal solution found. Trimming multipliers...")
+            if loglevel > 0 
+                println("\nOptimal solution found. Trimming multipliers...")
+            end
             g = b - λin*y
             trimx(x,S,R,active,state,g,b,λ,feaTol,optTol,loglevel)
             numtrim = nact - length(active)
@@ -352,19 +408,20 @@ function bpdual(
 
         end # if step < 1
 
-        if itn % 1000 == 0 #save every 1000 itns 
-            push!(tracer, (
-                Itn = itn,
-                Active = join(active, ","),
-                X = join(x, ","),
-                Lambda = λ,
-                Xnorm = xNorm,
-                Rnorm = rNorm,
-                CondS = condS))
-        end
+
+        push!(tracer.iteration, itn)
+        push!(tracer.lambda, λ)
+        push!(tracer.active, copy(active))
+        push!(tracer.activesoln, copy(x))
 
     end # while true
 
+    ## one last Update
+    push!(tracer.iteration, itn)
+    push!(tracer.lambda, λ)
+    push!(tracer.active, copy(active))
+    push!(tracer.activesoln, copy(x))
+
     tottime = time() - time0
     if loglevel > 0
         @printf("\n EXIT BPdual -- %s\n\n",EXIT_INFO[eFlag])
@@ -375,6 +432,6 @@ function bpdual(
         @printf(" %-20s: %8.1e %5s","Solution time (sec)",tottime,"")
         @printf("\n\n")
     end
-    return active,state,x,y,S,R,tracer
+    return tracer
 end # function bpdual
 

src/homotopy.jl

@@ -4,11 +4,9 @@ function asp_homotopy(A, b; min_lambda = 0.0,
     n = size(A, 2)
     bl = -ones(n)
     bu = +ones(n)
-    active, _, xx, _, _, _, tracer = bpdual(A, b, min_lambda, bl, bu;
+    tracer = bpdual(A, b, min_lambda, bl, bu;
                              homotopy = homotopy, 
                              kwargs...)
-    # BPdual's solution x is short. Make it full length.
-    x = zeros(n)
-    x[active] = xx
-    return x, tracer
+
+    return tracer
 end

test/runtests.jl

@@ -9,4 +9,6 @@ println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡
 
 @testset "ActiveSetPursuit.jl" begin
     @testset "BPDN" begin include("test_bpdn.jl") end
+    @testset "Recover decaying coefficients" begin include("test_recover_decaying.jl") end
 end
+

test/test_bpdn.jl

@@ -23,17 +23,12 @@ function test_bpdn()
     bu = +ones(n)
 
     # Solve the basis pursuit problem
-    active,state,xs,y,S,R,tracer = bpdual(A, b, 0., bl, bu, homotopy = false, loglevel =0)
-    xx = zeros(n)
-    xx[Int.(active)] = xs
-
+    tracer = bpdual(A, b, 0., bl, bu, homotopy = false, loglevel =0)
+
+    xx, λ = tracer[end]
     pFeas = norm(A * xx - b, Inf) / max(1, norm(xx, Inf))
-    dFeas = max(0, norm(A' * y, Inf) - 1)
-    dComp = abs(norm(xx, 1) - dot(b, y))
-
     @test pFeas <= 1e-6
-    @test dFeas <= 1e-6
-    @test dComp <= 1e-6
+
 end
 
 

test/test_recover_decaying.jl

@@ -0,0 +1,45 @@
+# ------------------------------------------------------------------
+#  Test Basis pursuit with decaying and permuted coefficients
+# ------------------------------------------------------------------ 
+
+using Test, LinearAlgebra, Random, SparseArrays, ActiveSetPursuit 
+
+function test_recover_decaying()
+    m = 600
+    n = 2560
+    threshold_percentage = 0.9   
+
+    # Generate solution with decaying coefficients and permutate
+    x = randn(n) ./ (1:n).^2
+    perm = randperm(n)
+    x = x[perm]
+
+    A = randn(m, n)
+
+    # Compute the RHS vector
+    b = A * x
+    bl = -ones(n)
+    bu = +ones(n)
+
+    # Solve the basis pursuit problem
+    tracer = asp_homotopy(A, b, min_lambda = 0.0, itnMax =300, loglevel =0)
+    xs, λ = tracer[end]
+
+    cumulative_norm = cumsum(abs.(x[sortperm(abs.(x), rev=true)]))
+    indices_to_recover = sortperm(abs.(x), rev=true)[1:findfirst(cumulative_norm .>= threshold_percentage * cumulative_norm[end])]
+
+    # Get the corresponding values in the recovered solution
+    recovered_indices = sortperm(abs.(xs), rev=true)[1:length(indices_to_recover)]
+
+    # Check if all the significant indices are correctly recovered
+    indices_correct = all(in(indices_to_recover).(recovered_indices))
+    values_correct = all((abs.(x[indices_to_recover]) .- abs.(xs[recovered_indices])) .< 1e-3)
+
+    @test indices_correct
+    @test values_correct
+end
+
+
+for ntest = 1:10
+    test_recover_decaying()
+end