Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

ASP #78

Merged
merged 8 commits into from
Aug 29, 2024
Merged

ASP #78

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
1efa62d
added homotopy
tinatorabi Aug 26, 2024
7feb94d
added ASP
tinatorabi Aug 26, 2024
6747520
added ASP
tinatorabi Aug 26, 2024
a2a6589
full path
tinatorabi Aug 27, 2024
e18f839
fixed ASP
tinatorabi Aug 27, 2024
93a1f32
added bysize
tinatorabi Aug 29, 2024
b9bcbb4
Add files via upload
tinatorabi Aug 29, 2024
b706e79
update
tinatorabi Aug 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 6 additions & 5 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ authors = ["William C Witt <[email protected]>, Christoph Ortner <christophor
version = "0.2.0"

[deps]
ActiveSetPursuit = "d86c1dc8-ba26-4c98-b330-3a8efc174d20"
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -21,25 +22,25 @@ MLJScikitLearnInterface = "5ae90465-5518-4432-b9d2-8a1def2f0cab"
PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"

[extensions]
ACEfit_PythonCall_ext = "PythonCall"
ACEfit_MLJLinearModels_ext = [ "MLJ", "MLJLinearModels" ]
ACEfit_MLJLinearModels_ext = ["MLJ", "MLJLinearModels"]
ACEfit_MLJScikitLearnInterface_ext = ["MLJScikitLearnInterface", "PythonCall", "MLJ"]
ACEfit_PythonCall_ext = "PythonCall"

[compat]
julia = "1.9"
IterativeSolvers = "0.9.2"
LowRankApprox = "0.5.3"
MLJ = "0.19"
MLJLinearModels = "0.9"
MLJScikitLearnInterface = "0.7"
LowRankApprox = "0.5.3"
Optim = "1.7"
ParallelDataTransfer = "0.5.0"
ProgressMeter = "1.7"
PythonCall = "0.9"
StaticArrays = "1.5"
julia = "1.9"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test", ]
test = ["Test"]
45 changes: 45 additions & 0 deletions src/solvers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ using LowRankApprox: pqrfact
using IterativeSolvers
using .BayesianLinear
using LinearAlgebra: SVD, svd
using ActiveSetPursuit

@doc raw"""
`struct QR` : linear least squares solver, using standard QR factorisation;
Expand Down Expand Up @@ -195,3 +196,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.
cortner marked this conversation as resolved.
Show resolved Hide resolved
* `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, traceFlag=true, kwargs...)
xs = tracer[end][1]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's also revisit this: I don't think that selecting the final iterate is normally desired. Maybe there could be a few options given by an additional parameters

ASP(; P = I, select = ( :byerror, 1.2 ), params... )
ASP(; P = I, select = ( :bysize, 1000 ), params... )
ASP(; P = I, select = ( :final, ), params... )

(:byerror, 1.2) would select the smallest active set fit that is within a factor 1.2 of the smallest fit error (not test error) achieved along the path. (:bysize, 1000) would select the smallest error fit that that less than 1000 active basis functions; while (:final, ) would do exactly what you did.

The only question is what the default should be. Maybe there should be no default and the user has to specify the mode?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess :bysize here doesn't really make sense unless we just integrate it into SmartASP and bring in test/validation sets, because using homotopy, the residual on the training set monotonically decreases so the smallest error fit that has <= 1000 active basis functions happens at actMax=1000.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

other than this, I think I've made all of the changes requested! Please let me know if I should change anything else :)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

monotonically decreases so the smallest error fit

ok, fair point. I need to think about it for a bit ...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should still add this option. I can envision that I want the solution with N parameters but have access to the model path up to 3*N parameters (or similar) for analysis ...

x_f = solver.P \ Array(xs)
println("done.")
return Dict{String, Any}("C" => x_f, "tracer" =>tracer)
cortner marked this conversation as resolved.
Show resolved Hide resolved
end
9 changes: 9 additions & 0 deletions test/test_linearsolvers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,3 +111,12 @@ 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"]
full_path = results["tracer"]
@show norm(A * C - y)
@show norm(C)
@show norm(C - c_ref)

Loading