QR refactoring

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

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
 
 

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()
 
@@ -133,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
@@ -198,7 +213,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 +236,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)

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))

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))