The source files for all examples can be found in /examples.
Lasso, Three Ways
This example shows how to use the MLSolver, QPSolver, and the GenericSolver interfaces to solve a lasso regression problem.
The lasso regression problem is
\[\begin{array}{ll} \text{minimize} & (1/2)\|Ax - b\|_2^2 + \gamma \|x\|_1 \end{array}\]
using GeNIOS
using Random, LinearAlgebra, SparseArraysGenerating the problem data
Random.seed!(1)
m, n = 200, 400
A = randn(m, n)
A .-= sum(A, dims=1) ./ m
normalize!.(eachcol(A))
xstar = sprandn(n, 0.1)
b = A*xstar + 1e-3*randn(m)
λ = 0.05*norm(A'*b, Inf)0.14048198582260793MLSolver interface
The MLSolver interface requires the per-sample loss and (optionally) its conjugate, if we want to use the dual gap as a stopping criterion (discussed in our paper). The conjugate function of $f$, defined as
\[f^*(y) = \sup_x \{yx - f(x)\},\]
Specifying the conjugate function is optional, and the solver will fall back to using the primal and dual residuals if it is not specified.
f(x) = 0.5*x^2
fconj(x) = 0.5*x^2
λ1 = λ
λ2 = 0.0
solver = GeNIOS.MLSolver(f, λ1, λ2, A, b; fconj=fconj)
res = solve!(solver; options=GeNIOS.SolverOptions(relax=true, use_dual_gap=true, dual_gap_tol=1e-3, verbose=true))
rmse = sqrt(1/m*norm(A*solver.zk - b, 2)^2)
println("Final RMSE: $(round(rmse, digits=8))")Starting setup...
Setup in 0.140s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Objective RMSE Dual Gap r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.995e+01 3.158e-01 9.256e+00 Inf Inf 1.000e+00 0.000
1 1.995e+01 3.158e-01 9.256e+00 0.000e+00 1.058e+01 1.000e+00 0.187
20 4.260e+00 5.491e-02 5.403e-02 2.000e-02 6.603e-02 1.000e+00 0.192
36 4.260e+00 5.408e-02 8.625e-04 4.108e-04 1.712e-03 1.000e+00 0.197
SOLVED in 0.197s, 36 iterations
Total time: 0.337s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.076475Note that the solver uses forward-mode automatic differentiation to compute the first and second derivatives of $f$. We can supply these arguments for additional speedup, if desired.
df(x) = x
d2f(x) = 1.0
solver = GeNIOS.MLSolver(f, df, d2f, λ1, λ2, A, b; fconj=fconj)
res = solve!(solver; options=GeNIOS.SolverOptions(relax=true, use_dual_gap=true, dual_gap_tol=1e-3, verbose=true))
rmse = sqrt(1/m*norm(A*solver.zk - b, 2)^2)
println("Final RMSE: $(round(rmse, digits=8))")Starting setup...
Setup in 0.038s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Objective RMSE Dual Gap r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.995e+01 3.158e-01 9.256e+00 Inf Inf 1.000e+00 0.000
1 1.995e+01 3.158e-01 9.256e+00 0.000e+00 1.058e+01 1.000e+00 0.076
20 4.260e+00 5.491e-02 5.403e-02 2.000e-02 6.603e-02 1.000e+00 0.081
36 4.260e+00 5.408e-02 8.625e-04 4.108e-04 1.712e-03 1.000e+00 0.087
SOLVED in 0.087s, 36 iterations
Total time: 0.125s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.076475QPSolver interface
The QPSolver interface requires us to specify the problem in the form
\[\begin{array}{ll} \text{minimize} & (1/2)x^TPx + q^Tx \\ \text{subject to} & l \leq Mx \leq u \end{array}\]
We introduce the new variable $t$ and inntroduce the constraint $-t \leq x \leq t$ to enforce the $\ell_1$ norm. The problem then becomes
\[\begin{array}{ll} \text{minimize} & (1/2)x^TA^TAx + b^TAx + \gamma \mathbf{1}^Tt \\ \text{subject to} & \begin{bmatrix}0 \\ -\infty\end{bmatrix} \leq \begin{bmatrix} I & I \\ I & -I \end{bmatrix} \begin{bmatrix}x \\ t\end{bmatrix} \leq \begin{bmatrix}\infty \\ 0\end{bmatrix} \end{array}\]
P = blockdiag(sparse(A'*A), spzeros(n, n))
q = vcat(-A'*b, λ*ones(n))
M = [
sparse(I, n, n) sparse(I, n, n);
sparse(I, n, n) -sparse(I, n, n)
]
l = [zeros(n); -Inf*ones(n)]
u = [Inf*ones(n); zeros(n)]
solver = GeNIOS.QPSolver(P, q, M, l, u)
res = solve!(
solver; options=GeNIOS.SolverOptions(
relax=true,
max_iters=1000,
eps_abs=1e-4,
eps_rel=1e-4,
verbose=true)
);
rmse = sqrt(1/m*norm(A*solver.xk[1:n] - b, 2)^2)
println("Final RMSE: $(round(rmse, digits=8))")Starting setup...
Setup in 0.000s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 0.000e+00 Inf Inf 1.000e+00 0.000
1 0.000e+00 0.000e+00 1.095e+01 1.000e+00 0.048
20 -1.530e+01 4.039e-02 2.905e-01 1.000e+00 0.073
40 -1.565e+01 1.575e-02 9.517e-02 1.000e+00 0.091
60 -1.568e+01 5.419e-03 3.517e-02 1.000e+00 0.112
80 -1.569e+01 2.035e-03 1.340e-02 1.000e+00 0.136
100 -1.569e+01 7.494e-04 4.885e-03 1.000e+00 0.161
109 -1.569e+01 4.726e-04 3.031e-03 1.000e+00 0.173
SOLVED in 0.173s, 109 iterations
Total time: 0.173s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.07684455Generic interface
Finally, we use the generic interface, which provides a large amount of control but is more complicated to use than the specialized interfaces demonstrated above.
First, we define the custom HessianOperator, which is used to solve the linear system in the $x$-subproblem. Since the Hessian of the objective is simply $A^TA$, this operator is simple for the Lasso problem. However, note that, since this is a custom type, we can speed up the multiplication to be more efficient than if we were to use $A^TA$ directly.
The update! function is called before the solution of the $x$-subproblem to update the Hessian, if necessary.
struct HessianLasso{T, S <: AbstractMatrix{T}} <: HessianOperator
A::S
vm::Vector{T}
end
function LinearAlgebra.mul!(y, H::HessianLasso, x)
mul!(H.vm, H.A, x)
mul!(y, H.A', H.vm)
return nothing
end
function update!(::HessianLasso, ::Solver)
return nothing
endupdate! (generic function with 1 method)Now, we define $f$, its gradient, $g$, and its proximal operator.
params = (; A=A, b=b, tmp=zeros(m), λ=λ)
function f(x, p)
A, b, tmp = p.A, p.b, p.tmp
mul!(tmp, A, x)
@. tmp -= b
return 0.5 * sum(w->w^2, tmp)
end
function grad_f!(g, x, p)
A, b, tmp = p.A, p.b, p.tmp
mul!(tmp, A, x)
@. tmp -= b
mul!(g, A', tmp)
return nothing
end
Hf = HessianLasso(A, zeros(m))
g(z, p) = p.λ*sum(x->abs(x), z)
function prox_g!(v, z, ρ, p)
λ = p.λ
@inline soft_threshold(x::T, κ::T) where {T <: Real} = sign(x) * max(zero(T), abs(x) - κ)
v .= soft_threshold.(z, λ/ρ)
endprox_g! (generic function with 1 method)Finally, we can solve the problem.
solver = GeNIOS.GenericSolver(
f, grad_f!, Hf, # f(x)
g, prox_g!, # g(z)
I, zeros(n); # M, c: Mx + z = c
params=params
)
res = solve!(solver; options=GeNIOS.SolverOptions(relax=true, verbose=true))
rmse = sqrt(1/m*norm(A*solver.zk - b, 2)^2)
println("Final RMSE: $(round(rmse, digits=8))")Starting setup...
Setup in 0.110s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.995e+01 Inf Inf 1.000e+00 0.000
1 1.995e+01 0.000e+00 1.058e+01 1.000e+00 0.108
20 4.268e+00 1.994e-02 2.279e-02 1.000e+00 0.113
30 4.260e+00 2.118e-03 8.626e-04 1.000e+00 0.116
SOLVED in 0.116s, 30 iterations
Total time: 0.226s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.07648496GPU support
Check out the page GPU Support for an example showing how to use GPUs with GeNIOS.
ProximalOperators.jl
We could alternatively use ProximalOperators.jl to define the proximal operator for g:
using ProximalOperators
prox_func = NormL1(λ)
gp(x, p) = prox_func(x)
prox_gp!(v, z, ρ, p) = prox!(v, prox_func, z, ρ)
# We see that this give the same result
solver = GeNIOS.GenericSolver(
f, grad_f!, Hf, # f(x)
gp, prox_gp!, # g(z)
I, zeros(n); # M, c: Mx + z = c
params=params
)
res = solve!(solver; options=GeNIOS.SolverOptions(relax=true, verbose=true))
rmse = sqrt(1/m*norm(A*solver.zk - b, 2)^2)
println("Final RMSE: $(round(rmse, digits=8))")Starting setup...
Setup in 0.008s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.995e+01 Inf Inf 1.000e+00 0.000
1 1.995e+01 0.000e+00 1.058e+01 1.000e+00 0.028
20 4.268e+00 1.994e-02 2.279e-02 1.000e+00 0.033
30 4.260e+00 2.118e-03 8.626e-04 1.000e+00 0.037
SOLVED in 0.037s, 30 iterations
Total time: 0.045s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.07648496This page was generated using Literate.jl.