Merge pull request #89 from ACEsuit/asp

ASP and OMP with D&C QR

Author: cortner

Project.toml

@@ -28,7 +28,7 @@ ACEfit_MLJScikitLearnInterface_ext = ["MLJScikitLearnInterface", "PythonCall", "
 ACEfit_PythonCall_ext = "PythonCall"
 
 [compat]
-ActiveSetPursuit = "0.0.2"
+ActiveSetPursuit = "0.0.4"
 IterativeSolvers = "0.9.2"
 LowRankApprox = "0.5.3"
 MLJ = "0.19"

src/asp.jl

@@ -1,3 +1,5 @@
+using Distributed
+using LinearAlgebra
 
 @doc raw"""
 
@@ -66,7 +68,18 @@ function solve(solver::ASP, A, y, Aval=A, yval=y)
     # Apply preconditioning
     AP = A / solver.P
     AvalP = Aval / solver.P
-    tracer = asp_homotopy(AP, y; solver.params..., traceFlag = true)
+
+    F = qr!(AP)
+    m, n = size(AP)
+    if m < n
+        error("ASP requires m >= n, but got m = $m, n = $n")
+    end
+    Qtb = F.Q' * y     
+    Qtb1 = Qtb[1:n]     
+
+    tracer = asp_homotopy(F.R, Qtb1; solver.params..., traceFlag = true)
+
+    # tracer = asp_homotopy(AP, y; solver.params..., traceFlag = true)
 
     q = length(tracer) 
     every = max(1, q / solver.nstore)
@@ -75,7 +88,7 @@ function solve(solver::ASP, A, y, Aval=A, yval=y)
                    for i in istore ]
 
     if solver.tsvd  # Post-processing if tsvd is true
-        post = post_asp_tsvd(new_tracer, AP, y, AvalP, yval)
+        post = post_asp_tsvd(new_tracer, F.R, Qtb1, AvalP, yval)
         new_post = [ (solution = solver.P \ p.θ, λ = p.λ, σ = p.σ) 
                      for p in post ]
     else
@@ -162,3 +175,105 @@ end
 #                 "path" => tracer, 
 #                 "nnzs" => length( (tracer[in][:solution]).nzind) )  
 # end
+
+
+
+@doc raw"""
+
+`OMP` : Orthogonal Matching 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), tsvd = false, nstore=100, 
+            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
+* `bval = nothing` : `p`- validation vector
+
+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 OMP
+    P
+    select
+    tsvd::Bool
+    nstore::Integer
+    params
+end
+
+function OMP(; P = I, select, tsvd=false, nstore=100, params...)
+    return OMP(P, select, tsvd, nstore, params)
+end
+
+function solve(solver::OMP, A, y, Aval=A, yval=y)
+    # Apply preconditioning
+    AP = A / solver.P
+    AvalP = Aval / solver.P
+
+    F = qr!(AP)
+    m, n = size(AP)
+    if m < n
+        error("OMP requires m >= n, but got m = $m, n = $n")
+    end
+    Qtb = F.Q' * y     
+    Qtb1 = Qtb[1:n]     
+
+    tracer = asp_omp(F.R, Qtb1, 0.0; traceFlag=true, loglevel=0, solver.params...)
+
+    q = length(tracer) 
+    every = max(1, q / solver.nstore)
+    istore = unique(round.(Int, [1:every:q; q]))
+    new_tracer = [ (solution = 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, F.R, Qtb1, AvalP, yval)
+        new_post = [ (solution = solver.P \ p.θ, λ = p.λ, σ = p.σ) 
+                     for p in post ]
+    else
+        new_post = [ (solution = solver.P \ p.solution, λ = p.λ, σ = 0.0) 
+                     for p in new_tracer ]
+    end
+
+    tracer_final = _add_errors(new_post, Aval, yval)
+    xs, in = asp_select(tracer_final, solver.select)
+
+    return Dict(   "C" => xs, 
+                "path" => tracer_final, )
+end