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, SparseArrays
Generating 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.10902103801870822
MLSolver 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.05093825
Note 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.05093825
QPSolver 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.05123142
Generic 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
end
update! (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, λ/ρ)
end
prox_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.05084529
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.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.05084529
This page was generated using Literate.jl.