From 6b84ab7b42df74b9a8d9cb84faba19c942937305 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Sat, 7 Sep 2024 01:59:19 -0700 Subject: [PATCH 1/7] added modes to ASP --- test/test_linearsolvers.jl | 59 ++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 28 deletions(-) diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index 4535ba9..296b04d 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -1,6 +1,7 @@ using ACEfit using LinearAlgebra +using Random using PythonCall @info("Test Solver on overdetermined system") @@ -111,32 +112,34 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) -@info(" ... ASP_homotopy selected by error") -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["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) -@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["path"] -@show results["nnzs"] -@show norm(A * C - y) -@show norm(C) -@show norm(C - c_ref) +@info(" ... ASP") +shuffled_indices = shuffle(1:length(y)) +train_indices = shuffled_indices[1:round(Int, 0.85 * length(y))] +val_indices = shuffled_indices[round(Int, 0.85 * length(y)) + 1:end] +At = A[train_indices,:] +Av = A[val_indices,:] +yt = y[train_indices] +yv = y[val_indices] + +for (sel, mod) in [((:final,nothing),:basic ),( (:byerror,1.3),:basic ),((:bysize,73),:basic ) + ,((:val,nothing),:smart ),((:byerror,1.3),:smart ),((:bysize,73),:smart )] + solver = ACEfit.ASP(P=I, select= sel, mode=mod ,params = (loglevel=0, traceFlag=true)) + if mod == :basic + results = ACEfit.solve(solver, A, y) + C = results["C"] + full_path = results["path"] + @show results["nnzs"] + @show norm(A * C - y) + @show norm(C) + @show norm(C - c_ref) + elseif mod == :smart + results = ACEfit.solve(solver, At, yt, Av, yv) + C = results["C"] + full_path = results["path"] + @show results["nnzs"] + @show norm(Av * C - yv) + @show norm(C) + @show norm(C - c_ref) + end +end From 1647ac600d73074e500935ce9d2fd76c4b8c5773 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Sat, 7 Sep 2024 02:00:18 -0700 Subject: [PATCH 2/7] added modes to ASP --- src/solvers.jl | 99 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 68 insertions(+), 31 deletions(-) diff --git a/src/solvers.jl b/src/solvers.jl index fc35934..1c37523 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -198,89 +198,86 @@ 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. - ``` +`struct ASP` : Active Set Pursuit solver + solves the optimization problem using either the basic or smart homotopy approach. * Input - * `A` : `m`-by-`n` explicit matrix or linear operator. + * `A` : `m`-by-`n` matrix. * `b` : `m`-vector. + * `Aval` : `p`-by-`n` validation matrix (only for smart mode). + * `bval` : `p`- validation vector (only for smart mode). * Solver parameters * `min_lambda` : Minimum value for `λ`. Defaults to zero if not provided. * `loglevel` : Logging level. * `itnMax` : Maximum number of iterations. * `actMax` : Maximum number of active constraints. + * `mode` : Either `:basic` or `:smart`. Constructor ```julia - ACEfit.ASP(; P = I, select, params) + ACEfit.ASPCombined(; P = I, select, mode, params) ``` where - `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 + mode::Symbol params::NamedTuple end -function ASP(; P = I, select, params...) +function ASP(; P = I, select, mode=:basic, params...) params_tuple = NamedTuple(params) - return ASP(P, select, params_tuple) + return ASP(P, select, mode, params_tuple) end -function solve(solver::ASP, A, y) +function solve(solver::ASP, A, y, Aval=nothing, yval=nothing) # Apply preconditioning AP = A / solver.P tracer = asp_homotopy(AP, y; solver.params[1]...) - - new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, length(tracer)) + q = length(tracer) + new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, q) - for i in 1:length(tracer) + for i in 1:q new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2]) end - # Select the final solution based on the criterion - xs, in = select_solution(new_tracer, solver, A, y) - + if solver.mode == :basic + xs, in = select_solution(new_tracer, solver, A, y) + elseif solver.mode == :smart + xs, in = select_smart(new_tracer, solver, Aval, yval) + else + @error("Unknown mode: $solver.mode") + end + println("done.") - return Dict("C" => xs, "path" => new_tracer, "nnzs" => length((tracer[in][1]).nzind) ) + return Dict("C" => xs, "path" => new_tracer, "nnzs" => length((new_tracer[in][:solution]).nzind) ) end function select_solution(tracer, solver, A, y) criterion, threshold = solver.select if criterion == :final - return tracer[end][1], length(tracer) + return tracer[end][:solution], length(tracer) elseif criterion == :byerror - errors = [norm(A * t[1] - y) for t in tracer] + errors = [norm(A * t[:solution] - 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 + return tracer[i][:solution], i end end elseif criterion == :bysize for i in 1:length(tracer) - if length((tracer[i][1]).nzind) == threshold - return tracer[i][1], i + if length((tracer[i][:solution]).nzind) == threshold + return tracer[i][:solution], i end end else @@ -288,3 +285,43 @@ function select_solution(tracer, solver, A, y) end end + +function select_smart(tracer, solver, Aval,yval) + + best_metric = Inf + best_iteration = 0 + validation_metric = 0 + q = length(tracer) + errors = [norm(Aval * t[:solution] - yval) for t in tracer] + nnzss = [(t[:solution]).nzind for t in tracer] + best_iteration = argmin(errors) + validation_metric = errors[best_iteration] + validation_end = norm(Aval * tracer[end][:solution] - yval) + + if validation_end < validation_metric #make sure to check the last one too in case q<<100 + best_iteration = q + end + + criterion, threshold = solver.select + + if criterion == :val + return tracer[best_iteration][:solution], best_iteration + + elseif criterion == :byerror + for (i, error) in enumerate(errors) + if error <= threshold * validation_metric + return tracer[i][:solution], i + end + end + + elseif criterion == :bysize + first_index = findfirst(sublist -> threshold in sublist, nnzss) + relevant_errors = errors[1:first_index - 1] + min_error = minimum(relevant_errors) + min_error_index = findfirst(==(min_error), relevant_errors) + return tracer[min_error_index][:solution], min_error_index + + else + @error("Unknown selection criterion: $criterion") + end +end From e77d05dcc0f2422f2801331cd72a1e2e4ac87640 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Sat, 7 Sep 2024 02:09:42 -0700 Subject: [PATCH 3/7] added modes to ASP --- src/solvers.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/solvers.jl b/src/solvers.jl index 1c37523..543811d 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -200,7 +200,14 @@ end @doc raw""" `struct ASP` : Active Set Pursuit solver solves the optimization problem using either the basic or smart 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` matrix. * `b` : `m`-vector. From 085245a9fe0a2aa803fdf7a94eec3839101b7359 Mon Sep 17 00:00:00 2001 From: Tina Torabi Date: Sat, 7 Sep 2024 02:12:00 -0700 Subject: [PATCH 4/7] added modes to ASP --- src/solvers.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers.jl b/src/solvers.jl index 543811d..39bb32f 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -223,7 +223,7 @@ end Constructor ```julia - ACEfit.ASPCombined(; P = I, select, mode, params) + ACEfit.ASP(; P = I, select, mode, params) ``` where - `P` : right-preconditioner / tychonov operator From 07da5e608ae7c9382c6ac7e58374b6c67574a17a Mon Sep 17 00:00:00 2001 From: tinatorabi Date: Sat, 7 Sep 2024 02:19:50 -0700 Subject: [PATCH 5/7] added random --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index f0feef2..2f7f5d2 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73" Optim = "429524aa-4258-5aef-a3af-852621145aeb" ParallelDataTransfer = "2dcacdae-9679-587a-88bb-8b444fb7085b" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" From 21a9f8b4666fe0927e42bd86783cccd00dac950a Mon Sep 17 00:00:00 2001 From: Christoph Ortner Date: Sat, 7 Sep 2024 11:28:13 -0700 Subject: [PATCH 6/7] move random to tests --- Project.toml | 1 - test/Project.toml | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 2f7f5d2..f0feef2 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,6 @@ LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73" Optim = "429524aa-4258-5aef-a3af-852621145aeb" ParallelDataTransfer = "2dcacdae-9679-587a-88bb-8b444fb7085b" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/test/Project.toml b/test/Project.toml index 2284da9..53ae33c 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,4 +4,5 @@ MLJ = "add582a8-e3ab-11e8-2d5e-e98b27df1bc7" MLJLinearModels = "6ee0df7b-362f-4a72-a706-9e79364fb692" MLJScikitLearnInterface = "5ae90465-5518-4432-b9d2-8a1def2f0cab" PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From eb557591a6068036d2c939344ed2ce18449fe8b8 Mon Sep 17 00:00:00 2001 From: Christoph Ortner Date: Sat, 7 Sep 2024 20:56:06 -0700 Subject: [PATCH 7/7] reorganize ASP a little --- src/ACEfit.jl | 1 + src/asp.jl | 157 +++++++++++++++++++++++++++++++++++++ src/solvers.jl | 136 -------------------------------- test/test_linearsolvers.jl | 111 ++++++++++++++++++++------ 4 files changed, 247 insertions(+), 158 deletions(-) create mode 100644 src/asp.jl diff --git a/src/ACEfit.jl b/src/ACEfit.jl index 86c6841..acd3e14 100644 --- a/src/ACEfit.jl +++ b/src/ACEfit.jl @@ -4,6 +4,7 @@ include("bayesianlinear.jl") include("data.jl") include("assemble.jl") include("solvers.jl") +include("asp.jl") include("fit.jl") end diff --git a/src/asp.jl b/src/asp.jl new file mode 100644 index 0000000..3de00e9 --- /dev/null +++ b/src/asp.jl @@ -0,0 +1,157 @@ + +@doc raw""" + +`ASP` : Active Set Pursuit solver + +Solves the lasso optimization problem. +```math +\max_{y} \left( b^T y - \frac{1}{2} λ y^T y \right) +``` +subject to +```math + \|A^T y\|_{\infty} \leq 1. +``` + +### Constructor Keyword arguments +```julia +ACEfit.ASP(; P = I, select = (:byerror, 1.0), + params...) +``` + +* `select` : Selection criterion for the final solution (required) + * `:final` : final solution (largest computed basis) + * `(:byerror, q)` : solution with error within `q` times the minimum error + along the path; if training error is used and `q == 1.0`, then this is + equivalent to to `:final`. + * `(:bysize, n)` : best solution with at most `n` non-zero features; if + training error is used, then it will be the solution with exactly `n` + non-zero features. +* `P = I` : prior / regularizer (optional) + +The remaining kwarguments to `ASP` are parameters for the ASP homotopy solver. + +* `actMax` : Maximum number of active constraints. +* `min_lambda` : Minimum value for `λ`. (defaults to 0) +* `loglevel` : Logging level. +* `itnMax` : Maximum number of iterations. + +### Extended syntax for `solve` + +```julia +solve(solver::ASP, A, y, Aval=A, yval=y) +``` +* `A` : `m`-by-`n` design matrix. (required) +* `b` : `m`-vector. (required) +* `Aval = nothing` : `p`-by-`n` validation matrix (only for `validate` mode). +* `bval = nothing` : `p`- validation vector (only for `validate` mode). + +If independent `Aval` and `yval` are provided (instead of detaults `A, y`), +then the solver will use this separate validation set instead of the training +set to select the best solution along the model path. +""" +struct ASP + P + select + params +end + +function ASP(; P = I, select, mode=:train, params...) + return ASP(P, select, params) +end + +function solve(solver::ASP, A, y, Aval=A, yval=y) + # Apply preconditioning + AP = A / solver.P + + tracer = asp_homotopy(AP, y; solver.params...) + q = length(tracer) + new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, q) + + for i in 1:q + new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2]) + end + + xs, in = select_solution(new_tracer, solver, Aval, yval) + + # println("done.") + return Dict( "C" => xs, + "path" => new_tracer, + "nnzs" => length((new_tracer[in][:solution]).nzind) ) +end + + +function select_solution(tracer, solver, A, y) + if solver.select == :final + criterion = :final + else + criterion, p = solver.select + end + + if criterion == :final + return tracer[end][:solution], length(tracer) + end + + if criterion == :byerror + maxind = length(tracer) + threshold = p + elseif criterion == :bysize + maxind = findfirst(t -> length((t[:solution]).nzind) > p, + tracer) - 1 + threshold = 1.0 + else + error("Unknown selection criterion: $criterion") + end + + errors = [ norm(A * t[:solution] - y) for t in tracer[1:maxind] ] + min_error = minimum(errors) + for (i, error) in enumerate(errors) + if error <= threshold * min_error + return tracer[i][:solution], i + end + end + + error("selection failed for unknown reasons; please file an issue with a MWE to reproduce this error.") +end + + +#= +function select_smart(tracer, solver, Aval, yval) + + best_metric = Inf + best_iteration = 0 + validation_metric = 0 + q = length(tracer) + errors = [norm(Aval * t[:solution] - yval) for t in tracer] + nnzss = [(t[:solution]).nzind for t in tracer] + best_iteration = argmin(errors) + validation_metric = errors[best_iteration] + validation_end = norm(Aval * tracer[end][:solution] - yval) + + if validation_end < validation_metric #make sure to check the last one too in case q<<100 + best_iteration = q + end + + criterion, threshold = solver.select + + if criterion == :val + return tracer[best_iteration][:solution], best_iteration + + elseif criterion == :byerror + for (i, error) in enumerate(errors) + if error <= threshold * validation_metric + return tracer[i][:solution], i + end + end + + elseif criterion == :bysize + first_index = findfirst(sublist -> threshold in sublist, nnzss) + relevant_errors = errors[1:first_index - 1] + min_error = minimum(relevant_errors) + min_error_index = findfirst(==(min_error), relevant_errors) + return tracer[min_error_index][:solution], min_error_index + + else + @error("Unknown selection criterion: $criterion") + end +end +=# \ No newline at end of file diff --git a/src/solvers.jl b/src/solvers.jl index 39bb32f..9867271 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -196,139 +196,3 @@ function solve(solver::TruncatedSVD, A, y) return Dict{String, Any}("C" => solver.P \ θP) end - -@doc raw""" -`struct ASP` : Active Set Pursuit solver - solves the optimization problem using either the basic or smart 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` matrix. - * `b` : `m`-vector. - * `Aval` : `p`-by-`n` validation matrix (only for smart mode). - * `bval` : `p`- validation vector (only for smart mode). - - * Solver parameters - * `min_lambda` : Minimum value for `λ`. Defaults to zero if not provided. - * `loglevel` : Logging level. - * `itnMax` : Maximum number of iterations. - * `actMax` : Maximum number of active constraints. - * `mode` : Either `:basic` or `:smart`. - - Constructor - ```julia - ACEfit.ASP(; P = I, select, mode, params) - ``` - where - - `P` : right-preconditioner / tychonov operator - - `select`: Selection mode for the final solution. - - `params`: The solver parameters, passed as named arguments. -""" -struct ASP - P::Any - select::Tuple - mode::Symbol - params::NamedTuple -end - -function ASP(; P = I, select, mode=:basic, params...) - params_tuple = NamedTuple(params) - return ASP(P, select, mode, params_tuple) -end - -function solve(solver::ASP, A, y, Aval=nothing, yval=nothing) - # Apply preconditioning - AP = A / solver.P - - tracer = asp_homotopy(AP, y; solver.params[1]...) - q = length(tracer) - new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, q) - - for i in 1:q - new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2]) - end - - if solver.mode == :basic - xs, in = select_solution(new_tracer, solver, A, y) - elseif solver.mode == :smart - xs, in = select_smart(new_tracer, solver, Aval, yval) - else - @error("Unknown mode: $solver.mode") - end - - println("done.") - return Dict("C" => xs, "path" => new_tracer, "nnzs" => length((new_tracer[in][:solution]).nzind) ) -end - -function select_solution(tracer, solver, A, y) - criterion, threshold = solver.select - - if criterion == :final - return tracer[end][:solution], length(tracer) - - elseif criterion == :byerror - errors = [norm(A * t[:solution] - y) for t in tracer] - min_error = minimum(errors) - - for (i, error) in enumerate(errors) - if error <= threshold * min_error - return tracer[i][:solution], i - end - end - elseif criterion == :bysize - for i in 1:length(tracer) - if length((tracer[i][:solution]).nzind) == threshold - return tracer[i][:solution], i - end - end - else - @error("Unknown selection criterion: $criterion") - end -end - - -function select_smart(tracer, solver, Aval,yval) - - best_metric = Inf - best_iteration = 0 - validation_metric = 0 - q = length(tracer) - errors = [norm(Aval * t[:solution] - yval) for t in tracer] - nnzss = [(t[:solution]).nzind for t in tracer] - best_iteration = argmin(errors) - validation_metric = errors[best_iteration] - validation_end = norm(Aval * tracer[end][:solution] - yval) - - if validation_end < validation_metric #make sure to check the last one too in case q<<100 - best_iteration = q - end - - criterion, threshold = solver.select - - if criterion == :val - return tracer[best_iteration][:solution], best_iteration - - elseif criterion == :byerror - for (i, error) in enumerate(errors) - if error <= threshold * validation_metric - return tracer[i][:solution], i - end - end - - elseif criterion == :bysize - first_index = findfirst(sublist -> threshold in sublist, nnzss) - relevant_errors = errors[1:first_index - 1] - min_error = minimum(relevant_errors) - min_error_index = findfirst(==(min_error), relevant_errors) - return tracer[min_error_index][:solution], min_error_index - - else - @error("Unknown selection criterion: $criterion") - end -end diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index 296b04d..234be82 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -1,10 +1,14 @@ using ACEfit -using LinearAlgebra +using LinearAlgebra, Random, Test using Random using PythonCall +## + @info("Test Solver on overdetermined system") + +Random.seed!(1234) Nobs = 10_000 Nfeat = 100 A1 = randn(Nobs, Nfeat) / sqrt(Nobs) @@ -12,9 +16,12 @@ U, S1, V = svd(A1) S = 1e-4 .+ ((S1 .- S1[end]) / (S1[1] - S1[end])).^2 A = U * Diagonal(S) * V' c_ref = randn(Nfeat) -y = A * c_ref + 1e-3 * randn(Nobs) / sqrt(Nobs) +epsn = 1e-5 +y = A * c_ref + epsn * randn(Nobs) / sqrt(Nobs) P = Diagonal(1.0 .+ rand(Nfeat)) +## + @info(" ... QR") solver = ACEfit.QR() results = ACEfit.solve(solver, A, y) @@ -22,6 +29,10 @@ C = results["C"] @show norm(A * C - y) @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 100 * epsn + +## @info(" ... regularised QR, λ = 1e-5") solver = ACEfit.QR(lambda = 1e-5, P = P) @@ -30,6 +41,11 @@ C = results["C"] @show norm(A * C - y) @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 1000 * epsn + + +## @info(" ... regularised QR, λ = 1e-2") solver = ACEfit.QR(lambda = 1e-2, P = P) @@ -39,6 +55,11 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 1 +@test norm(C - c_ref) < 10 + +## + @info(" ... RRQR, rtol = 1e-15") solver = ACEfit.RRQR(rtol = 1e-15, P = P) results = ACEfit.solve(solver, A, y) @@ -47,6 +68,10 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 100 * epsn + +## @info(" ... RRQR, rtol = 1e-5") solver = ACEfit.RRQR(rtol = 1e-5, P = P) @@ -55,6 +80,10 @@ C = results["C"] @show norm(A * C - y) @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 100 * epsn + +## @info(" ... RRQR, rtol = 1e-3") solver = ACEfit.RRQR(rtol = 1e-3, P = P) @@ -64,6 +93,11 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 1 +@test norm(C - c_ref) < 1 + +## + @info(" ... LSQR") solver = ACEfit.LSQR(damp = 0, atol = 1e-6) results = ACEfit.solve(solver, A, y) @@ -72,6 +106,11 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 1 + +## + @info(" ... SKLEARN_BRR") solver = ACEfit.SKLEARN_BRR() results = ACEfit.solve(solver, A, y) @@ -80,6 +119,8 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +## + @info(" ... SKLEARN_ARD") solver = ACEfit.SKLEARN_ARD() results = ACEfit.solve(solver, A, y) @@ -88,6 +129,8 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +## + @info(" ... BLR") solver = ACEfit.BLR() results = ACEfit.solve(solver, A, y) @@ -96,6 +139,11 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 1 + +## + @info(" ... TruncatedSVD(; rtol = 1e-5)") solver = ACEfit.TruncatedSVD(; rtol = 1e-5) results = ACEfit.solve(solver, A, y) @@ -104,6 +152,11 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 100 * epsn + +## + @info(" ... TruncatedSVD(; rtol = 1e-4)") solver = ACEfit.TruncatedSVD(; rtol=1e-4) results = ACEfit.solve(solver, A, y) @@ -112,6 +165,10 @@ C = results["C"] @show norm(C) @show norm(C - c_ref) +@test norm(A * C - y) < 10 * epsn +@test norm(C - c_ref) < 1 + +## @info(" ... ASP") shuffled_indices = shuffle(1:length(y)) @@ -122,24 +179,34 @@ Av = A[val_indices,:] yt = y[train_indices] yv = y[val_indices] -for (sel, mod) in [((:final,nothing),:basic ),( (:byerror,1.3),:basic ),((:bysize,73),:basic ) - ,((:val,nothing),:smart ),((:byerror,1.3),:smart ),((:bysize,73),:smart )] - solver = ACEfit.ASP(P=I, select= sel, mode=mod ,params = (loglevel=0, traceFlag=true)) - if mod == :basic - results = ACEfit.solve(solver, A, y) - C = results["C"] - full_path = results["path"] - @show results["nnzs"] - @show norm(A * C - y) - @show norm(C) - @show norm(C - c_ref) - elseif mod == :smart - results = ACEfit.solve(solver, At, yt, Av, yv) - C = results["C"] - full_path = results["path"] - @show results["nnzs"] - @show norm(Av * C - yv) - @show norm(C) - @show norm(C - c_ref) - end +for (select, tolr, tolc) in [ (:final, 10*epsn, 1), + ( (:byerror,1.3), 10*epsn, 1), + ( (:bysize,73), 1, 10) ] + @show select + local solver, results, C + solver = ACEfit.ASP(P=I, select = select, loglevel=0, traceFlag=true) + # without validation + results = ACEfit.solve(solver, A, y) + C = results["C"] + full_path = results["path"] + @show results["nnzs"] + @show norm(A * C - y) + @show norm(C) + @show norm(C - c_ref) + + @test norm(A * C - y) < tolr + @test norm(C - c_ref) < tolc + + + # with validation + results = ACEfit.solve(solver, At, yt, Av, yv) + C = results["C"] + full_path = results["path"] + @show results["nnzs"] + @show norm(Av * C - yv) + @show norm(C) + @show norm(C - c_ref) + + @test norm(Av * C - yv) < tolr + @test norm(C - c_ref) < tolc end