From a64d8d8aa4644747279664df86fa48383fc582e6 Mon Sep 17 00:00:00 2001 From: Christoph Ortner Date: Sun, 8 Sep 2024 21:49:53 -0700 Subject: [PATCH 1/3] experimental tsvd postprocessing --- test/test_asp.jl | 113 +++++++++++++++++++++++++++++++++++++ test/test_linearsolvers.jl | 42 -------------- 2 files changed, 113 insertions(+), 42 deletions(-) create mode 100644 test/test_asp.jl diff --git a/test/test_asp.jl b/test/test_asp.jl new file mode 100644 index 0000000..1b81703 --- /dev/null +++ b/test/test_asp.jl @@ -0,0 +1,113 @@ +using ACEfit +using LinearAlgebra, Random, Test +using Random + +## + +@info("Test Solver on overdetermined system") + +Random.seed!(1234) +Nobs = 10_000 +Nfeat = 300 +A1 = randn(Nobs, Nfeat) / sqrt(Nobs) +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) +epsn = 1e-5 +y = A * c_ref + epsn * randn(Nobs) / sqrt(Nobs) +P = Diagonal(1.0 .+ rand(Nfeat)) + +## + +@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 (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 + +## + +# Experimental Implementation of tsvd postprocessing + + +using SparseArrays + +function solve_tsvd(At, yt, Av, yv) + Ut, Σt, Vt = svd(At); zt = Ut' * yt + Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv + @assert issorted(Σt, rev=true) + + Rv_Vt = Rv * Vt + + θv = zeros(size(Av, 2)) + θv[1] = zt[1] / Σt[1] + rv = Rv_Vt[:, 1] * θv[1] - zv + + tsvd_errs = Float64[] + push!(tsvd_errs, norm(rv)) + + for k = 2:length(Σt) + θv[k] = zt[k] / Σt[k] + rv += Rv_Vt[:, k] * θv[k] + push!(tsvd_errs, norm(rv)) + end + + imin = argmin(tsvd_errs) + θv[imin+1:end] .= 0 + return Vt * θv, Σt[imin] +end + +function post_asp_tsvd(path, At, yt, Av, yv) + Qt, Rt = qr(At); zt = Matrix(Qt)' * yt + Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv + + post = [] + for (θ, λ) in path + if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end + inz = θ.nzind + θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) + θ2 = copy(θ); θ2[inz] .= θ1 + push!(post, (θ = θ2, λ = λ, σ = σ)) + end + return identity.(post) +end + +solver = ACEfit.ASP(P=I, select = :final, loglevel=0, traceFlag=true) +result = ACEfit.solve(solver, At, yt); +post = post_asp_tsvd(result["path"], At, yt, Av, yv); diff --git a/test/test_linearsolvers.jl b/test/test_linearsolvers.jl index 234be82..9c03b53 100644 --- a/test/test_linearsolvers.jl +++ b/test/test_linearsolvers.jl @@ -168,45 +168,3 @@ C = results["C"] @test norm(A * C - y) < 10 * epsn @test norm(C - c_ref) < 1 -## - -@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 (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 From 59723c5a044442f3b939045a47ec04a790cecc97 Mon Sep 17 00:00:00 2001 From: tinatorabi Date: Mon, 9 Sep 2024 06:01:19 -0700 Subject: [PATCH 2/3] integrated tsvd --- Project.toml | 1 + src/asp.jl | 99 ++++++++++++++++++++----------------- test/test_asp.jl | 126 +++++++++++++++++++++++++++++++++-------------- 3 files changed, 144 insertions(+), 82 deletions(-) diff --git a/Project.toml b/Project.toml index f0feef2..c49320f 100644 --- a/Project.toml +++ b/Project.toml @@ -13,6 +13,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb" ParallelDataTransfer = "2dcacdae-9679-587a-88bb-8b444fb7085b" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [weakdeps] diff --git a/src/asp.jl b/src/asp.jl index 3de00e9..56de1b4 100644 --- a/src/asp.jl +++ b/src/asp.jl @@ -48,15 +48,19 @@ solve(solver::ASP, A, y, Aval=A, yval=y) 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 + mode::Symbol + tsvd::Bool + nstore::Integer params end -function ASP(; P = I, select, mode=:train, params...) - return ASP(P, select, params) +function ASP(; P = I, select, mode=:train, tsvd=false, nstore=100, params...) + return ASP(P, select, mode, tsvd, nstore, params) end function solve(solver::ASP, A, y, Aval=A, yval=y) @@ -64,18 +68,24 @@ function solve(solver::ASP, A, y, Aval=A, yval=y) AP = A / solver.P tracer = asp_homotopy(AP, y; solver.params...) + q = length(tracer) + every = max(1, q ÷ solver.nstore) new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, q) + new_tracer = [(solution = solver.P \ tracer[i][1], λ = tracer[i][2]) for i in [1:every:q; q]] - for i in 1:q - new_tracer[i] = (solution = solver.P \ tracer[i][1], λ = tracer[i][2]) + if solver.tsvd # Post-processing if tsvd is true + post = post_asp_tsvd(new_tracer, A, y, Aval, yval) + new_post = [(solution = p.θ, λ = p.λ) for p in post] + else + new_post = new_tracer end - xs, in = select_solution(new_tracer, solver, Aval, yval) + xs, in = select_solution(new_post, solver, Aval, yval) # println("done.") - return Dict( "C" => xs, - "path" => new_tracer, + return Dict( "C" => xs, + "path" => new_post, "nnzs" => length((new_tracer[in][:solution]).nzind) ) end @@ -114,44 +124,45 @@ function select_solution(tracer, solver, A, y) 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 +using SparseArrays - 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 +function solve_tsvd(At, yt, Av, yv) + Ut, Σt, Vt = svd(At); zt = Ut' * yt + Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv + @assert issorted(Σt, rev=true) - 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 + Rv_Vt = Rv * Vt - else - @error("Unknown selection criterion: $criterion") - end + θv = zeros(size(Av, 2)) + θv[1] = zt[1] / Σt[1] + rv = Rv_Vt[:, 1] * θv[1] - zv + + tsvd_errs = Float64[] + push!(tsvd_errs, norm(rv)) + + for k = 2:length(Σt) + θv[k] = zt[k] / Σt[k] + rv += Rv_Vt[:, k] * θv[k] + push!(tsvd_errs, norm(rv)) + end + + imin = argmin(tsvd_errs) + θv[imin+1:end] .= 0 + return Vt * θv, Σt[imin] +end + +function post_asp_tsvd(path, At, yt, Av, yv) + Qt, Rt = qr(At); zt = Matrix(Qt)' * yt + Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv + + post = [] + for (θ, λ) in path + if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end + inz = θ.nzind + θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) + θ2 = copy(θ); θ2[inz] .= θ1 + push!(post, (θ = θ2, λ = λ, σ = σ)) + end + return identity.(post) end -=# \ No newline at end of file diff --git a/test/test_asp.jl b/test/test_asp.jl index 1b81703..7fc993c 100644 --- a/test/test_asp.jl +++ b/test/test_asp.jl @@ -61,53 +61,103 @@ for (select, tolr, tolc) in [ (:final, 10*epsn, 1), @test norm(C - c_ref) < tolc end +## + + +# I didn't wanna add more tsvd tests to yours so I just wrote this one +# I only wanted to naïvely demonstrate that tsvd actually does make a difference! :) + +for (select, tolr, tolc) in [ (:final, 20*epsn, 1.5), + ( (:byerror,1.3), 20*epsn, 1.5), + ( (:bysize,73), 1, 10) ] + @show select + local solver, results, C + solver_tsvd = ACEfit.ASP(P=I, select=select, mode=:train, tsvd=true, + nstore=100, loglevel=0, traceFlag=true) + + solver = ACEfit.ASP(P=I, select=select, mode=:train, tsvd=false, + nstore=100, loglevel=0, traceFlag=true) + # without validation + results_tsvd = ACEfit.solve(solver_tsvd, A, y) + results = ACEfit.solve(solver, A, y) + C_tsvd = results_tsvd["C"] + C = results["C"] + + @show results["nnzs"] + @show norm(A * C - y) + @show norm(A * C_tsvd - y) + if norm(A * C_tsvd - y)< norm(A * C - y) + @info "tsvd made improvements!" + else + @warn "tsvd did NOT make any improvements!" + end + + + # with validation + results_tsvd = ACEfit.solve(solver_tsvd, At, yt, Av, yv) + results = ACEfit.solve(solver, At, yt, Av, yv) + C_tsvd = results_tsvd["C"] + C = results["C"] + @show results["nnzs"] + @show norm(A * C - y) + @show norm(A * C_tsvd - y) + + if norm(A * C_tsvd - y)< norm(A * C - y) + @info "tsvd made improvements!" + else + @warn "tsvd did NOT make any improvements!" + end +end + + ## # Experimental Implementation of tsvd postprocessing -using SparseArrays +# using SparseArrays -function solve_tsvd(At, yt, Av, yv) - Ut, Σt, Vt = svd(At); zt = Ut' * yt - Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv - @assert issorted(Σt, rev=true) +# function solve_tsvd(At, yt, Av, yv) +# Ut, Σt, Vt = svd(At); zt = Ut' * yt +# Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv +# @assert issorted(Σt, rev=true) - Rv_Vt = Rv * Vt +# Rv_Vt = Rv * Vt - θv = zeros(size(Av, 2)) - θv[1] = zt[1] / Σt[1] - rv = Rv_Vt[:, 1] * θv[1] - zv +# θv = zeros(size(Av, 2)) +# θv[1] = zt[1] / Σt[1] +# rv = Rv_Vt[:, 1] * θv[1] - zv - tsvd_errs = Float64[] - push!(tsvd_errs, norm(rv)) +# tsvd_errs = Float64[] +# push!(tsvd_errs, norm(rv)) - for k = 2:length(Σt) - θv[k] = zt[k] / Σt[k] - rv += Rv_Vt[:, k] * θv[k] - push!(tsvd_errs, norm(rv)) - end +# for k = 2:length(Σt) +# θv[k] = zt[k] / Σt[k] +# rv += Rv_Vt[:, k] * θv[k] +# push!(tsvd_errs, norm(rv)) +# end - imin = argmin(tsvd_errs) - θv[imin+1:end] .= 0 - return Vt * θv, Σt[imin] -end +# imin = argmin(tsvd_errs) +# θv[imin+1:end] .= 0 +# return Vt * θv, Σt[imin] +# end + +# function post_asp_tsvd(path, At, yt, Av, yv) +# Qt, Rt = qr(At); zt = Matrix(Qt)' * yt +# Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv + +# post = [] +# for (θ, λ) in path +# if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end +# inz = θ.nzind +# θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) +# θ2 = copy(θ); θ2[inz] .= θ1 +# push!(post, (θ = θ2, λ = λ, σ = σ)) +# end +# return identity.(post) +# end + +# solver = ACEfit.ASP(P=I, select = :final, loglevel=0, traceFlag=true) +# result = ACEfit.solve(solver, At, yt); +# post = post_asp_tsvd(result["path"], At, yt, Av, yv); -function post_asp_tsvd(path, At, yt, Av, yv) - Qt, Rt = qr(At); zt = Matrix(Qt)' * yt - Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv - - post = [] - for (θ, λ) in path - if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end - inz = θ.nzind - θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) - θ2 = copy(θ); θ2[inz] .= θ1 - push!(post, (θ = θ2, λ = λ, σ = σ)) - end - return identity.(post) -end - -solver = ACEfit.ASP(P=I, select = :final, loglevel=0, traceFlag=true) -result = ACEfit.solve(solver, At, yt); -post = post_asp_tsvd(result["path"], At, yt, Av, yv); From 83dbee521fb8261b4c8cbcbfc292bce4abd19772 Mon Sep 17 00:00:00 2001 From: Christoph Ortner Date: Mon, 9 Sep 2024 12:42:49 -0700 Subject: [PATCH 3/3] cleanup, activate asp tests --- src/asp.jl | 37 +++++++++++++++++++++------------ test/runtests.jl | 2 ++ test/test_asp.jl | 54 +----------------------------------------------- 3 files changed, 27 insertions(+), 66 deletions(-) diff --git a/src/asp.jl b/src/asp.jl index 56de1b4..88a2aa7 100644 --- a/src/asp.jl +++ b/src/asp.jl @@ -69,24 +69,24 @@ function solve(solver::ASP, A, y, Aval=A, yval=y) tracer = asp_homotopy(AP, y; solver.params...) - q = length(tracer) - every = max(1, q ÷ solver.nstore) - new_tracer = Vector{NamedTuple{(:solution, :λ), Tuple{Any, Any}}}(undef, q) - new_tracer = [(solution = solver.P \ tracer[i][1], λ = tracer[i][2]) for i in [1:every:q; q]] + q = length(tracer) + every = max(1, q ÷ solver.nstore) + istore = unique([1:every:q; q]) + new_tracer = [ (solution = solver.P \ tracer[i][1], λ = tracer[i][2], σ = 0.0 ) + for i in istore ] if solver.tsvd # Post-processing if tsvd is true post = post_asp_tsvd(new_tracer, A, y, Aval, yval) - new_post = [(solution = p.θ, λ = p.λ) for p in post] + new_post = [ (solution = p.θ, λ = p.λ, σ = p.σ) for p in post ] else new_post = new_tracer end xs, in = select_solution(new_post, solver, Aval, yval) - # println("done.") return Dict( "C" => xs, "path" => new_post, - "nnzs" => length((new_tracer[in][:solution]).nzind) ) + "nnzs" => length( (new_tracer[in][:solution]).nzind) ) end @@ -156,13 +156,24 @@ function post_asp_tsvd(path, At, yt, Av, yv) Qt, Rt = qr(At); zt = Matrix(Qt)' * yt Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv - post = [] - for (θ, λ) in path - if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end + function _post(θλ) + (θ, λ) = θλ + if isempty(θ.nzind); return (θ = θ, λ = λ, σ = Inf); end inz = θ.nzind θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) θ2 = copy(θ); θ2[inz] .= θ1 - push!(post, (θ = θ2, λ = λ, σ = σ)) - end - return identity.(post) + return (θ = θ2, λ = λ, σ = σ) + end + + return _post.(path) + +# post = [] +# for (θ, λ) in path +# if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end +# inz = θ.nzind +# θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) +# θ2 = copy(θ); θ2[inz] .= θ1 +# push!(post, (θ = θ2, λ = λ, σ = σ)) +# end +# return identity.(post) end diff --git a/test/runtests.jl b/test/runtests.jl index 2039593..1c46419 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,5 +8,7 @@ using Test @testset "Linear Solvers" begin include("test_linearsolvers.jl") end + @testset "ASP" begin include("test_asp.jl") end + @testset "MLJ Solvers" begin include("test_mlj.jl") end end diff --git a/test/test_asp.jl b/test/test_asp.jl index 7fc993c..4cacb64 100644 --- a/test/test_asp.jl +++ b/test/test_asp.jl @@ -8,7 +8,7 @@ using Random Random.seed!(1234) Nobs = 10_000 -Nfeat = 300 +Nfeat = 100 A1 = randn(Nobs, Nfeat) / sqrt(Nobs) U, S1, V = svd(A1) S = 1e-4 .+ ((S1 .- S1[end]) / (S1[1] - S1[end])).^2 @@ -109,55 +109,3 @@ for (select, tolr, tolc) in [ (:final, 20*epsn, 1.5), end end - -## - -# Experimental Implementation of tsvd postprocessing - - -# using SparseArrays - -# function solve_tsvd(At, yt, Av, yv) -# Ut, Σt, Vt = svd(At); zt = Ut' * yt -# Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv -# @assert issorted(Σt, rev=true) - -# Rv_Vt = Rv * Vt - -# θv = zeros(size(Av, 2)) -# θv[1] = zt[1] / Σt[1] -# rv = Rv_Vt[:, 1] * θv[1] - zv - -# tsvd_errs = Float64[] -# push!(tsvd_errs, norm(rv)) - -# for k = 2:length(Σt) -# θv[k] = zt[k] / Σt[k] -# rv += Rv_Vt[:, k] * θv[k] -# push!(tsvd_errs, norm(rv)) -# end - -# imin = argmin(tsvd_errs) -# θv[imin+1:end] .= 0 -# return Vt * θv, Σt[imin] -# end - -# function post_asp_tsvd(path, At, yt, Av, yv) -# Qt, Rt = qr(At); zt = Matrix(Qt)' * yt -# Qv, Rv = qr(Av); zv = Matrix(Qv)' * yv - -# post = [] -# for (θ, λ) in path -# if isempty(θ.nzind); push!(post, (θ = θ, λ = λ, σ = Inf)); continue; end -# inz = θ.nzind -# θ1, σ = solve_tsvd(Rt[:, inz], zt, Rv[:, inz], zv) -# θ2 = copy(θ); θ2[inz] .= θ1 -# push!(post, (θ = θ2, λ = λ, σ = σ)) -# end -# return identity.(post) -# end - -# solver = ACEfit.ASP(P=I, select = :final, loglevel=0, traceFlag=true) -# result = ACEfit.solve(solver, At, yt); -# post = post_asp_tsvd(result["path"], At, yt, Av, yv); -