Nonmonotone variant of linesearch methods

[aldma]

This PR adds the possibility to run algorithms with a linesearch procedure in a nonmonotone variant. In particular,

• this feature addresses algorithms such as PANOC, PANOCplus, ZeroFPR and DRLS;
• a kwarg monotonicity is added to control the level of monotonicity required, with 0 < monotonicity <= 1;
• the default monotonicity = 1 coincides with the current implementation (monotone descent);
• a scalar value merit is included in the state, monitoring progress instead of objective or envelope values.

[lostella]

@aldma thanks, I should be able to review this over the weekend.

I’m wondering about the test failures: the test dependency on AbstractOperators we could remove; regarding RecursiveArrayTools, not sure, I’ll need to look into it.

[lostella]

Rebasing should get rid of the failing tests!

[lostella]

Why do the assignment FBE_x = state.merit here? It seems to be giving an inappropriate name to this quantity. Should you not directly use state.merit where needed in the lines that follow?

[lostella]

Also wondering whether merit is a good name here: this quantity appears to be an exponential moving average of the values of the merit function we're minimizing, is that correct?

[aldma]

You are right, I just wanted to minimize the changes relative to the monotone version. The assignment FBE_x = state.merit has been removed now with dbdd5c1, where possible.
In ZeroFPR and PANOCplus I left an assignment to FBE_x before the loop because the tests where failing otherwise; I did not understand why the variable FBE_x was not available at the final update of state.merit.

[aldma]

Yes, and the value of merit is the one considered as baseline for accepting or rejecting an update in the linesearch. Also the objective function and the envelope are used as merit functions, but they can decrease nonmonotonically, whereas what is called merit decreases monotonically also in the nonmonotone variant.
I am happy to change nomenclature if you have a more appropriate one.

[aldma]

I hope I didn't mess things up with the rebasing. I am still learning gitology.

[aldma]

Hi @lostella , anything missing before this can be merged?

[lostella]

@aldma it seems like there is slightly different behavior in ZeroFPR when monotonicity=1 and adaptiveness is on. See number of iters from the attached script, in the two branches. Not necessarily better or worse, but it seems like there is just different behavior by default.

On master:

problem                            algorithm  iters  status objective
lasso, underdetermined 30x60       PANOC        130  ok     5.059791265
lasso, underdetermined 30x60       ZeroFPR       61  ok     5.059791265
lasso, overdetermined 60x30        PANOC         39  ok     19.84802782
lasso, overdetermined 60x30        ZeroFPR       21  ok     19.84802782
sparse logistic 40x100             PANOC        631  ok     2.311038119
sparse logistic 40x100             ZeroFPR      164  ok     2.311038119
sparse logistic, ill-cond. 80x200  PANOC       1997  ok     2.737235383
sparse logistic, ill-cond. 80x200  ZeroFPR      521  ok     2.737235383
nonconvex box QP n=30              PANOC         32  ok     -28.31841284
nonconvex box QP n=30              ZeroFPR       12  ok     -28.31841284

On PR 96:

problem                            algorithm  iters  status objective
lasso, underdetermined 30x60       PANOC        130  ok     5.059791265
lasso, underdetermined 30x60       ZeroFPR       63  ok     5.059791265
lasso, overdetermined 60x30        PANOC         39  ok     19.84802782
lasso, overdetermined 60x30        ZeroFPR       23  ok     19.84802782
sparse logistic 40x100             PANOC        631  ok     2.311038119
sparse logistic 40x100             ZeroFPR      161  ok     2.311038119
sparse logistic, ill-cond. 80x200  PANOC       1997  ok     2.737235383
sparse logistic, ill-cond. 80x200  ZeroFPR      543  ok     2.737235383
nonconvex box QP n=30              PANOC         32  ok     -28.31841284
nonconvex box QP n=30              ZeroFPR       12  ok     -28.31841284

To reproduce (note: Claude-authored):

using ProximalAlgorithms
using ProximalCore
using LinearAlgebra
using Random
using Printf

const MAXIT = 5000
const TOL = 1e-8

# smooth terms, with analytic gradients (extension point `value_and_gradient`)

struct LeastSquaresLoss{M,V}   # 1/2 ||A x - b||^2
    A::M
    b::V
end
(f::LeastSquaresLoss)(x) = (r = f.A * x - f.b; dot(r, r) / 2)
function ProximalAlgorithms.value_and_gradient(f::LeastSquaresLoss, x)
    r = f.A * x - f.b
    return dot(r, r) / 2, f.A' * r
end

struct LogisticLoss{M,V}       # sum_i log(1 + exp(-(A x - b)_i))
    A::M
    b::V
end
softplus(z) = max(-z, zero(z)) + log1p(exp(-abs(z)))   # log(1 + exp(-z)), overflow-safe
(f::LogisticLoss)(x) = sum(softplus, f.A * x - f.b)
function ProximalAlgorithms.value_and_gradient(f::LogisticLoss, x)
    z = f.A * x - f.b
    return sum(softplus, z), -(f.A' * (1 ./ (1 .+ exp.(z))))
end

struct QuadraticLoss{M,V}      # 1/2 x'Qx + q'x
    Q::M
    q::V
end
(f::QuadraticLoss)(x) = dot(x, f.Q * x) / 2 + dot(f.q, x)
function ProximalAlgorithms.value_and_gradient(f::QuadraticLoss, x)
    Qx = f.Q * x
    return dot(x, Qx) / 2 + dot(f.q, x), Qx + f.q
end

# nonsmooth terms, with their proximal mappings

struct NormL1{T}               # lambda ||x||_1
    lambda::T
end
(g::NormL1)(x) = g.lambda * norm(x, 1)
function ProximalCore.prox!(y, g::NormL1, x, gamma)
    y .= sign.(x) .* max.(abs.(x) .- gamma * g.lambda, 0)
    return g(y)
end

struct IndBox{T}               # indicator of [lo, hi]^n
    lo::T
    hi::T
end
(g::IndBox)(x) = all(xi -> g.lo <= xi <= g.hi, x) ? zero(g.lo) : oftype(g.lo, Inf)
function ProximalCore.prox!(y, g::IndBox, x, gamma)
    y .= clamp.(x, g.lo, g.hi)
    return zero(g.lo)
end

function lasso(seed, m, n)
    rng = Random.Xoshiro(seed)
    A, b = randn(rng, m, n), randn(rng, m)
    return (f = LeastSquaresLoss(A, b), g = NormL1(0.1 * norm(A' * b, Inf)), n = n)
end

function logistic(seed, m, n; lam, decay = 0.0)
    rng = Random.Xoshiro(seed)
    A = randn(rng, m, n) * Diagonal(10.0 .^ range(0, -decay, length = n))
    b = randn(rng, m)
    return (f = LogisticLoss(A, b), g = NormL1(lam), n = n)
end

function boxqp(seed, n)
    rng = Random.Xoshiro(seed)
    S = randn(rng, n, n) / sqrt(n)
    Q = (S + S') / 2                      # symmetric indefinite: nonconvex
    return (f = QuadraticLoss(Q, randn(rng, n)), g = IndBox(-1.0, 1.0), n = n)
end

const PROBLEMS = [
    ("lasso, underdetermined 30x60", lasso(1, 30, 60)),
    ("lasso, overdetermined 60x30", lasso(2, 60, 30)),
    ("sparse logistic 40x100", logistic(3, 40, 100; lam = 0.1)),
    ("sparse logistic, ill-cond. 80x200", logistic(4, 80, 200; lam = 0.02, decay = 1.5)),
    ("nonconvex box QP n=30", boxqp(5, 30)),
]

const ALGORITHMS = [
    ("PANOC", ProximalAlgorithms.PANOC),
    ("ZeroFPR", ProximalAlgorithms.ZeroFPR),
]

@printf("%-34s %-9s %6s  %-6s %s\n", "problem", "algorithm", "iters", "status", "objective")

for (pname, p) in PROBLEMS, (aname, alg) in ALGORITHMS
    solver = alg(; tol = TOL, maxit = MAXIT, adaptive = true)
    x, it = solver(; x0 = zeros(p.n), f = p.f, g = p.g)
    @printf("%-34s %-9s %6d  %-6s %.10g\n", pname, aname, it,
            it >= MAXIT ? "MAXIT" : "ok", p.f(x) + p.g(x))
end

[aldma]

Thanks @lostella ! Bug fixed and the behavior seems the same on master and #96 with the script above. Do you know why the benchmark action is failing now?

[lostella]

@aldma looks like this may be due to the Julia 1.13 release: developing the package in the benchmark workflow causes Project.toml to be edited, leaving the repo in a dirty state. Then PkgBenchmark refuses to compare branches from a dirty state.

I'll propose a fix, then we can rebase this.