From 1efa62da05bc087a65f4e6e2db750785766743ac Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Mon, 26 Aug 2024 16:37:15 -0700 Subject: [PATCH 1/8] added homotopy --- test/test_linearsolvers.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index 3a410a6..32ed249 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -111,3 +111,10 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@info(" ... ASP_homotopy") +solver = ACEfit.ASP(P = P) +results = ACEfit.solve(solver, A, y) +C = results["C"] +@show norm(A * C - y) +@show norm(C) +@show norm(C - c_ref) From 7feb94df7f5088b46a466f75cb5d2e0eae8eb41c Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Mon, 26 Aug 2024 16:37:50 -0700 Subject: [PATCH 2/8] added ASP --- src/solvers.jl | 44 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/solvers.jl b/src/solvers.jl index 077d17b..74bf49d 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -195,3 +195,47 @@ function solve(solver::TruncatedSVD, A, y) return Dict{String, Any}("C" => solver.P \ θP) end + +@doc raw""" +`struct ASP` : Active Set Pursuit sparse solver + solves the following optimization problem using the homotopy approach: + + ```math + \max_{y} \left( b^T y - \frac{1}{2} λ y^T y \right) + ``` + subject to + + ```math + \|A^T y\|_{\infty} \leq 1. + ``` + + * Input + * `A` : `m`-by-`n` explicit matrix or linear operator. + * `b` : `m`-vector. + * `min_lambda` : Minimum value for `λ`. Is set to zero if not input is given. + * `loglevel` : Logging level. + * `itnMax` : Maximum number of iterations. + * `feaTol` : Feasibility tolerance. + * `actMax` : Maximum number of active constraints. + + Constructor + ```julia + ACEfit.ASP(; P = I) + ``` + where + * `P` : right-preconditioner / tychonov operator +""" +struct ASP + P::Any +end + +ASP(; P = I) = ASP(P) + +function solve(solver::ASP, A, y; kwargs...) + AP = A / solver.P + tracer = asp_homotopy(AP, y; loglevel=0, kwargs...) + xs = tracer[end][1] + x_f = solver.P \ Array(xs) + println("done.") + return Dict{String, Any}("C" => x_f) +end From 6747520417af37dddd830ea5a4e26a68916fcf5b Mon Sep 17 00:00:00 2001 From: tinatorabi Date: Mon, 26 Aug 2024 16:51:33 -0700 Subject: [PATCH 3/8] added ASP --- Project.toml | 11 ++++++----- src/solvers.jl | 1 + 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 0cb99a6..29d3c82 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["William C Witt , Christoph Ortner Date: Mon, 26 Aug 2024 17:49:57 -0700 Subject: [PATCH 4/8] full path --- src/solvers.jl | 4 ++-- test/test_linearsolvers.jl | 2 ++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index 29eabe2..76c9157 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -234,9 +234,9 @@ ASP(; P = I) = ASP(P) function solve(solver::ASP, A, y; kwargs...) AP = A / solver.P - tracer = asp_homotopy(AP, y; loglevel=0, kwargs...) + tracer = asp_homotopy(AP, y; loglevel=0, traceFlag=true, kwargs...) xs = tracer[end][1] x_f = solver.P \ Array(xs) println("done.") - return Dict{String, Any}("C" => x_f) + return Dict{String, Any}("C" => x_f, "tracer" =>tracer) end diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index 32ed249..b81616e 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -115,6 +115,8 @@ C = results["C"] solver = ACEfit.ASP(P = P) results = ACEfit.solve(solver, A, y) C = results["C"] +full_path = results["tracer"] @show norm(A * C - y) @show norm(C) @show norm(C - c_ref) + From e18f839063eda41a96943d7251d65241c5d0771f Mon Sep 17 00:00:00 2001 From: tinatorabi Date: Tue, 27 Aug 2024 00:43:15 -0700 Subject: [PATCH 5/8] fixed ASP --- src/solvers.jl | 62 ++++++++++++++++++++++++++++++++------ test/test_linearsolvers.jl | 14 +++++++-- 2 files changed, 64 insertions(+), 12 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index 76c9157..fad579f 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -213,30 +213,72 @@ end * Input * `A` : `m`-by-`n` explicit matrix or linear operator. * `b` : `m`-vector. - * `min_lambda` : Minimum value for `λ`. Is set to zero if not input is given. + + * Solver parameters + * `min_lambda` : Minimum value for `λ`. Defaults to zero if not provided. * `loglevel` : Logging level. * `itnMax` : Maximum number of iterations. - * `feaTol` : Feasibility tolerance. * `actMax` : Maximum number of active constraints. Constructor ```julia - ACEfit.ASP(; P = I) + ACEfit.ASP(; P = I, select, params) ``` where - * `P` : right-preconditioner / tychonov operator + - `P` : right-preconditioner / tychonov operator + - `select`: Selection mode for the final solution. + - `(:byerror, th)`: Selects the smallest active set fit within a factor `th` of the smallest fit error. + - `(:final, nothing)`: Returns the final iterate. + - `params`: The solver parameters, passed as named arguments. """ struct ASP P::Any + select::Tuple + params::NamedTuple end -ASP(; P = I) = ASP(P) +function ASP(; P = I, select, params...) + params_tuple = NamedTuple(params) + return ASP(P, select, params_tuple) +end -function solve(solver::ASP, A, y; kwargs...) +function solve(solver::ASP, A, y) + # Apply preconditioning AP = A / solver.P - tracer = asp_homotopy(AP, y; loglevel=0, traceFlag=true, kwargs...) - xs = tracer[end][1] - x_f = solver.P \ Array(xs) + + tracer = asp_homotopy(AP, y; solver.params[1]...) + + new_tracer = Vector{typeof(tracer[1])}(undef, length(tracer)) + for i in 1:length(tracer) + new_tracer[i] = (solver.P \ Array(tracer[i][1]), tracer[i][2]) + end + + # Select the final solution based on the criterion + xs, in = select_solution(new_tracer, solver, A, y) + x_f = Array(xs) + println("done.") - return Dict{String, Any}("C" => x_f, "tracer" =>tracer) + return Dict("C" => x_f, "tracer" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) ) end + +function select_solution(tracer, solver, A, y) + criterion, threshold = solver.select + + if criterion == :final + return tracer[end][1], length(tracer) + + elseif criterion == :byerror + errors = [norm(A * t[1] - y) for t in tracer] + min_error = minimum(errors) + + # Find the solution with the smallest error within the threshold + for (i, error) in enumerate(errors) + if error <= threshold * min_error + return tracer[i][1], i + end + end + else + @error("Unknown selection criterion: $criterion") + end +end + diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index b81616e..614c6f4 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -111,12 +111,22 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) -@info(" ... ASP_homotopy") -solver = ACEfit.ASP(P = P) +@info(" ... ASP_homotopy selected by error") +solver = ACEfit.ASP(P = P, select = (:byerror,1.2), params = (loglevel=0, traceFlag=true)) results = ACEfit.solve(solver, A, y) C = results["C"] full_path = results["tracer"] +@show results["nnzs"] @show norm(A * C - y) @show norm(C) @show norm(C - c_ref) +@info(" ... ASP_homotopy final solution") +solver = ACEfit.ASP(P = P, select = (:final,nothing), params = (loglevel=0, traceFlag=true)) +results = ACEfit.solve(solver, A, y) +C = results["C"] +full_path = results["tracer"] +@show results["nnzs"] +@show norm(A * C - y) +@show norm(C) +@show norm(C - c_ref) From 93a1f325e053f50d08957d5cbc02c333d9297872 Mon Sep 17 00:00:00 2001 From: tinatorabi Date: Wed, 28 Aug 2024 21:13:13 -0700 Subject: [PATCH 6/8] added bysize --- src/solvers.jl | 13 ++++++++++--- test/test_linearsolvers.jl | 16 +++++++++++++--- 2 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index fad579f..9a38de6 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -248,9 +248,10 @@ function solve(solver::ASP, A, y) tracer = asp_homotopy(AP, y; solver.params[1]...) - new_tracer = Vector{typeof(tracer[1])}(undef, length(tracer)) + new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, length(tracer)) + for i in 1:length(tracer) - new_tracer[i] = (solver.P \ Array(tracer[i][1]), tracer[i][2]) + new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2]) end # Select the final solution based on the criterion @@ -258,7 +259,7 @@ function solve(solver::ASP, A, y) x_f = Array(xs) println("done.") - return Dict("C" => x_f, "tracer" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) ) + return Dict("C" => x_f, "path" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) ) end function select_solution(tracer, solver, A, y) @@ -277,6 +278,12 @@ function select_solution(tracer, solver, A, y) return tracer[i][1], i end end + elseif criterion == :bysize + for i in 1:length(tracer) + if length((tracer[i][1]).nzind) == threshold + return tracer[i][1], i + end + end else @error("Unknown selection criterion: $criterion") end diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index 614c6f4..4535ba9 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -112,10 +112,20 @@ C = results["C"] @show norm(C - c_ref) @info(" ... ASP_homotopy selected by error") -solver = ACEfit.ASP(P = P, select = (:byerror,1.2), params = (loglevel=0, traceFlag=true)) +solver = ACEfit.ASP(P = P, select = (:byerror,1.5), params = (loglevel=0, traceFlag=true)) results = ACEfit.solve(solver, A, y) C = results["C"] -full_path = results["tracer"] +full_path = results["path"] +@show results["nnzs"] +@show norm(A * C - y) +@show norm(C) +@show norm(C - c_ref) + +@info(" ... ASP_homotopy selected by size") +solver = ACEfit.ASP(P = P, select = (:bysize,50), params = (loglevel=0, traceFlag=true)) +results = ACEfit.solve(solver, A, y) +C = results["C"] +full_path = results["path"] @show results["nnzs"] @show norm(A * C - y) @show norm(C) @@ -125,7 +135,7 @@ full_path = results["tracer"] solver = ACEfit.ASP(P = P, select = (:final,nothing), params = (loglevel=0, traceFlag=true)) results = ACEfit.solve(solver, A, y) C = results["C"] -full_path = results["tracer"] +full_path = results["path"] @show results["nnzs"] @show norm(A * C - y) @show norm(C) From b9bcbb47c7643340229b7d62c67177b37e7f44e5 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Thu, 29 Aug 2024 09:45:55 -0700 Subject: [PATCH 7/8] Add files via upload From b706e7943527c8a72650ecfc6c027473f85fd6f2 Mon Sep 17 00:00:00 2001 From: tinatorabi Date: Thu, 29 Aug 2024 11:26:50 -0700 Subject: [PATCH 8/8] update --- src/solvers.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index 9a38de6..fc35934 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -256,10 +256,9 @@ function solve(solver::ASP, A, y) # Select the final solution based on the criterion xs, in = select_solution(new_tracer, solver, A, y) - x_f = Array(xs) println("done.") - return Dict("C" => x_f, "path" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) ) + return Dict("C" => xs, "path" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) ) end function select_solution(tracer, solver, A, y)