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.10902103801870822MLSolver 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.112s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Objective RMSE Dual Gap r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.682e+01 2.900e-01 9.256e+00 Inf Inf 1.000e+00 0.000
1 1.682e+01 2.900e-01 9.256e+00 0.000e+00 9.692e+00 1.000e+00 0.206
20 2.963e+00 3.613e-02 7.393e-02 1.897e-02 6.309e-02 1.000e+00 0.211
35 2.963e+00 3.602e-02 6.247e-04 6.305e-04 1.270e-03 1.000e+00 0.216
SOLVED in 0.216s, 35 iterations
Total time: 0.328s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.05093825Note 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.043s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Objective RMSE Dual Gap r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.682e+01 2.900e-01 9.256e+00 Inf Inf 1.000e+00 0.000
1 1.682e+01 2.900e-01 9.256e+00 0.000e+00 9.692e+00 1.000e+00 0.090
20 2.963e+00 3.613e-02 7.393e-02 1.897e-02 6.309e-02 1.000e+00 0.094
35 2.963e+00 3.602e-02 6.247e-04 6.305e-04 1.270e-03 1.000e+00 0.099
SOLVED in 0.099s, 35 iterations
Total time: 0.142s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.05093825QPSolver 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 9.934e+00 1.000e+00 0.001
20 -1.337e+01 3.609e-02 2.994e-01 1.000e+00 0.026
40 -1.380e+01 2.097e-02 1.168e-01 1.000e+00 0.044
60 -1.385e+01 6.188e-03 3.146e-02 1.000e+00 0.066
80 -1.385e+01 1.398e-03 6.595e-03 1.000e+00 0.091
89 -1.385e+01 6.617e-04 2.935e-03 1.000e+00 0.102
SOLVED in 0.102s, 89 iterations
Total time: 0.103s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.05123142Generic 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.083s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.682e+01 Inf Inf 1.000e+00 0.000
1 1.682e+01 0.000e+00 9.692e+00 1.000e+00 0.076
20 2.970e+00 1.897e-02 1.327e-02 1.000e+00 0.082
27 2.963e+00 2.518e-03 1.789e-03 1.000e+00 0.084
SOLVED in 0.084s, 27 iterations
Total time: 0.167s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.05084529ProximalOperators.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.010s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Iteration Obj Val r_primal r_dual ρ Time
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
0 1.682e+01 Inf Inf 1.000e+00 0.000
1 1.682e+01 0.000e+00 9.692e+00 1.000e+00 0.027
20 2.970e+00 1.897e-02 1.327e-02 1.000e+00 0.033
27 2.963e+00 2.518e-03 1.789e-03 1.000e+00 0.035
SOLVED in 0.035s, 27 iterations
Total time: 0.045s
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Final RMSE: 0.05084529This page was generated using Literate.jl.