Coordinate descent, losses and options

coordinate_descent!(x, f, g, options=CDOptions(); rng=Random.default_rng(), kwargs...)
solve!(x, f, g, options::CDOptions; rng=Random.default_rng(), kwargs...)

Minimize f(x) + g(x) by coordinate descent and store the solution in x. The loss f is one of CDLeastSquaresLoss, CDWeightedLSLoss, CDSqrtLassoLoss and CDQuadraticLoss, a loss of ProximalBase that ProximalSolvers's coordinate descent takes (LeastSquaresLoss, a least-squares GLMLoss, QuadraticFunction, KronQuadratic, possibly scaled as in (1/n) * LeastSquaresLoss(y, X)), a model of ProximalSolvers (ProximalSolvers.CDModel), or a user-defined CoordinateDifferentiableFunction. The penalty g separates across coordinates (NormL1, ZeroFunction, ElasticNet, MCP, SCAD, … from ProximalBase, or another penalty, see below) and x is a SparseIterate or a SymmetricSparseIterate, or g is a group lasso penalty GroupNormL2 and x an AtomIterate (see below).

The descent is the coordinate descent of ProximalSolvers (ProximalSolvers.CDSolver, called through ProximalSolvers.solve!) with the options of CDOptions. It alternates between full passes over all coordinates and passes over the nonzero coordinates only (the active set). It stops when a full pass is followed by a pass in which no coordinate changes by more than options.opt_tol, or after options.max_iter passes. In the second case x may be far from the minimizer, and a warning is issued.

If options.warm_start is true, the descent starts from the supplied x. Otherwise x is set to zero and, for NormL1 and GroupNormL2, the problem is solved along a geometric sequence of options.num_steps penalty levels that decreases from the smallest λ0 with zero solution to g.λ0, each one warm started from the previous solution. Coordinates with zero penalty weight are left out when finding that λ0. For a SparseIterate or an AtomIterate, each penalty level is screened with the sequential strong rule, as in LassoPath.

The keyword arguments kwargs are fields of CDOptions that replace those of options, as in coordinate_descent!(x, f, g; opt_tol=1e-8). With order = :random, the coordinates are visited in an order drawn from rng; pass a seeded generator, such as Xoshiro(1), for reproducible passes without reseeding the global one.

coordinate_descent! and solve! (a method of ProximalSolvers.solve!) return a ProximalSolvers.CDResults: the number of passes (including those on the path of penalties when warm_start is false), whether the stopping rule was met (converged), the largest scaled change of a coordinate in the last pass (maxchange) and the tolerance it was compared with, the objective minimum, and the fitted iterate minimizer. Passing a CDInfo as the fifth argument of the deprecated coordinateDescent! is supported for compatibility.

A dense Vector x is also accepted; it is converted to a SparseIterate, or to an AtomIterate for GroupNormL2, and overwritten with the solution.

Group lasso

For g = GroupNormL2(λ0, groups, w), g(x) = λ0 Σₖ wₖ ‖x[Gₖ]‖₂, the coordinates are the groups: x is an AtomIterate with the groups of g, such as AtomIterate(p, g.groups), and the active set is its set of active groups. Each step minimizes over one group G the quadratic with curvature L that majorizes f around x, plus the penalty: x[G] ← (1 - λ0 wₖ / (L ‖u‖₂))₊ u with u = x[G] - ∇_G f(x) / L (groupwise majorization descent; Yang and Zou, A fast unified algorithm for solving group-lasso penalized learning problems, 2015). For CDLeastSquaresLoss, with or without intercept, L is the largest eigenvalue of X[:, G]'X[:, G] / n (centered), and for CDQuadraticLoss that of A[G, G]; user-defined losses do not take groups. The stopping rule takes the largest change of an entry in the pass, as for single coordinates. The path starts from maxₖ ‖∇_Gₖ f(0)‖₂ / wₖ over the penalized groups, and the strong rule keeps the groups with ‖∇_G f‖₂ ≥ (2λ - λprev) wₖ. Groups with zero weight are not penalized, and the groups must cover all coordinates, so a coordinate that is not penalized needs a group of its own with zero weight.

Other penalties

A separable penalty g(x) = Σₖ gₖ(xₖ) of another type, a subtype of ProximalBase.AbstractFunction with ProximalBase.isseparable true for it, needs the proximal operator of one entry, ProximalBase.prox_entry(g, v, k, γ), which returns argmin_z gₖ(z) + (z - v)² / (2γ). The built-in losses take this proximal step on the quadratic in x[k], so they accept such a penalty, and so does a user-defined loss whose descend_coordinate! calls ProximalBase.prox_coordinate! (see CoordinateDifferentiableFunction). Such a penalty has no path of penalty levels: with warm_start = false, the descent starts from zero and solves at g directly.

Anderson acceleration

With anderson = K in CDOptions, which needs order = :cyclic, the passes over the active set are accelerated (Bertrand and Massias, Anderson acceleration of coordinate descent, 2021): every K + 1 passes over an unchanged active set, the iterates are extrapolated and the extrapolated point is kept if it lowers the objective. It needs a SparseIterate or a SymmetricSparseIterate, and a loss whose value the model knows (the built-in losses, or a callable user-defined loss). It pays off when coordinate descent needs many passes, as at small penalties and with correlated columns. It keeps the weakness of the cyclic order when the columns share a strong common factor, where order = :random is much faster.

Momentum

With momentum (momentum in CDOptions), each pass over the active set starts from the extrapolated point xₖ + βₖ (xₖ - xₖ₋₁), with Nesterov's βₖ = (k - 1) / (k + 2); the momentum restarts when the step of a pass points against it (O'Donoghue and Candès, 2015), after a full pass, and when the active set shrinks. The cache of a built-in loss moves with the point, at the cost of a few vector operations per pass; a user-defined loss is initialized at the extrapolated point. It needs a SparseIterate, a SymmetricSparseIterate or, with ProximalSolvers 0.3.1, an AtomIterate, whose active groups are extrapolated entry by entry while the curvatures of the groups are kept. Unlike Anderson acceleration, it works with the random order, so it also helps when the columns share a common factor; with the cyclic order it can stall.

By default (momentum = nothing) it is used with order = :random, a SparseIterate and a built-in loss other than KronQuadratic, and not otherwise; on an AtomIterate it stays off, except that multitask_lasso and MultitaskLassoPath ask for it on their stacked design. momentum = true asks for it with a user-defined loss or an AtomIterate as well; with a user-defined loss it can pay off when initializing the loss costs less than a pass. momentum = false turns it off.

solve!(x, f, g, options::CDOptions; rng=Random.default_rng(), kwargs...)

The method of ProximalSolvers.solve! for the options of this package, identical to coordinate_descent!. It returns a ProximalSolvers.CDResults, with the number of passes, whether the stopping rule was met (converged), the largest scaled change of a coordinate in the last pass (maxchange) and the tolerance it was compared with.

CoordinateDifferentiableFunction

Supertype of user-defined smooth losses f that coordinate_descent! can minimize. A subtype caches whatever it needs (for example residuals) to make coordinate updates cheap, and implements

• initialize!(f, x): set the cache for the iterate x
• gradient(f, x, k): coordinate k of the gradient of f at x
• num_coordinates(f): number of coordinates
• descend_coordinate!(f, g, x, k): update coordinate k of x and the cache

and optionally coordinate_scale(f, k) and response_scale(f), which make the stopping rule independent of the units of the data (see CDOptions), full_gradient!(grad, f, x), a faster way to compute the whole gradient, and f(x), the value of f at x from the cache, which Anderson acceleration needs (anderson in CDOptions). These functions are not exported; extend them as HighDimensionalModels.initialize! and so on. The four required functions have no methods of their own, so a missing method raises a MethodError that lists the methods that exist.

The penalty g enters only through descend_coordinate!. A loss whose step ends with the proximal operator of g on coordinate k, ProximalBase.prox_coordinate!(x, g, k, γ), works with every separable penalty of ProximalBase and with user-defined ones (see coordinate_descent!).

Coordinate descent runs on ProximalSolvers, which takes such a loss through a model (ProximalSolvers.CDModel); a loss written for ProximalSolvers directly, as a CDModel, is also accepted.

CDLeastSquaresLoss(y, X; intercept=false)

The least squares loss f(β) = ‖y - Xβ‖² / (2n), as the model of ProximalSolvers for (1/n) * LeastSquaresLoss(y, X) (see ProximalSolvers.cdmodel). It caches the residuals r = y - Xβ, so a coordinate update costs O(n), and reads X and y in place.

With intercept, the loss is minimized over an unpenalized intercept as well, f(β) = min_β₀ ‖y - β₀ - Xβ‖² / (2n), whose minimizer is β₀ = ȳ - x̄'β for the column means x̄ and the mean ȳ of y (ProximalSolvers.intercept(f, β)). The columns and y are centered on the fly, so X is not copied, and the residuals are y - β₀ - Xβ (ProximalSolvers.residuals(f)).

A column that the loss does not depend on, a zero column or, with intercept, a constant one, gets the minimizer of the penalty nearest zero. y and X must not contain NaN or Inf.

On an AtomIterate, for the group lasso, a group G is updated at once, with the largest eigenvalue of the centered X[:, G]'X[:, G] / n as the curvature.

CDQuadraticLoss(A, b)

The quadratic f(x) = x'Ax / 2 + b'x with symmetric A, as the model of ProximalSolvers for QuadraticFunction(A, b). It caches Ax, so a coordinate update costs O(p). With A = X'X/n and b = -X'y/n it gives the same minimizer as CDLeastSquaresLoss.

A coordinate with A[k, k] = 0, which f does not depend on when A is positive semidefinite, goes to the minimizer of the penalty nearest zero. A and b must not contain NaN or Inf. A dense Matrix A is used as it is, as in version 0.6; any other matrix is copied into one.

On an AtomIterate, for the group lasso, a group G is updated at once, with the largest eigenvalue of A[G, G] as the curvature.

CDSqrtLassoLoss(y, X; intercept=false)

The square-root lasso loss f(β) = ‖y - Xβ‖, as the model of ProximalSolvers for ProximalSolvers.SqrtLeastSquaresLoss(y, X). It caches the residuals r = y - Xβ and the squared column norms. With intercept, f(β) = min_β₀ ‖y - β₀ - Xβ‖, minimized by β₀ = ȳ - x̄'β; as in CDLeastSquaresLoss, the columns and y are centered on the fly.

This loss is not smooth where r = 0, so coordinate descent minimizes f + g over a coordinate exactly instead of using a quadratic approximation. It takes the penalty NormL1 only.

CDSqrtQuadraticLoss(A, b, c; check=true)

The square-root lasso loss from sufficient statistics, f(x) = √(x'Ax + 2b'x + c), with symmetric positive semidefinite A and a nonnegative quadratic, as the model of ProximalSolvers for ProximalSolvers.SqrtQuadraticFunction(A, b, c; check). With A = X'X/n, b = -X'y/n and c = y'y/n it is ‖y - Xx‖ / √n, the loss of CDSqrtLassoLoss divided by √n, without the data: it caches Ax and the value of the quadratic, so a coordinate update costs O(p) instead of O(n), and checking a coordinate that does not move costs O(1). It pays off when n > p, and when many fits share A. It also solves the scaled lasso of one variable on the others from their covariance matrix S: A = S, b = -S[:, j] and c = S[j, j], with an infinite penalty weight on coordinate j, which keeps it at zero (see scaled_lasso!).

As for CDSqrtLassoLoss, coordinate descent minimizes f + g over a coordinate exactly, for NormL1 only. A, b and c must not contain NaN or Inf, and A must be symmetric; checking A costs O(p²), as much as a few passes, so check = false skips it when many losses share a checked A, as the columns of a covariance matrix do. A dense Matrix A is used as it is.

CDWeightedLSLoss(y, X, w)

The weighted least squares loss f(β) = ∑ᵢ wᵢ (yᵢ - Xᵢ'β)² / (2n), as the model of ProximalSolvers for (1/n) * GLMLoss(ProximalBase.Squared(), y, X; weights = w). It caches the residuals r = y - Xβ and the weighted squared column norms.

y, X and w are read in place, so changing them changes the loss for the next solve. They are checked for NaN or Inf and the weights for negative entries only when the loss is built, so build it from initialized buffers. A coefficient whose column does not enter the loss (∑ᵢ wᵢ Xᵢₖ² = 0) goes to the minimizer of the penalty nearest zero.

coordinate_scale(f, k)

The scale of coordinate k for the relative stopping rule (see opt_tol in CDOptions): the change of x[k] is multiplied by it. Defaults to 1.

descend_coordinate!(f, g, x, k)

Update coordinate k of x to decrease f(x) + g(x), and return the change of the coordinate. The cache of f must be updated for the new x.

full_gradient!(grad, f, x)

Write the gradient of f at x into grad, from the cache of f, which must be current for x, and return grad. Used to find the largest useful penalty level when warm_start is false. Defaults to gradient(f, x, k) for every coordinate.

gradient(f, x, k)

Coordinate k of the gradient of f at x, computed from the cache of f, which must be current for x.

initialize!(f, x)

Set the internal cache of f (for example the residuals) for the iterate x. Called by coordinate_descent! before the first pass.

num_coordinates(f)

Number of coordinates (or blocks of coordinates) over which coordinate descent iterates.

response_scale(f)

The scale that opt_tol is multiplied by in the relative stopping rule (see CDOptions). Defaults to 1.

CDInfo()

Convergence information that the deprecated coordinateDescent! fills in when it is passed as the last argument. Deprecated: coordinate_descent!(x, f, g, options) (see coordinate_descent!) returns the same information, and more, in a ProximalSolvers.CDResults: passes for iterations, converged, maxchange for max_change, tolerance, and the objective minimum.

• iterations: number of passes over the coordinates, including those on the path of penalties when warm_start is false
• converged: whether the stopping rule was met before max_iter passes
• max_change: largest change of a coordinate in the last pass, scaled as described for opt_tol in CDOptions
• tolerance: the threshold max_change is compared with

CDOptions(; max_iter=100_000, opt_tol=1e-7, order=:random, warm_start=true,
            num_steps=50, relative_tol=true, anderson=0, momentum=nothing, gaptol=nothing)
CDOptions(o::CDOptions; changes...)

Options for coordinate_descent! and the lasso solvers, which run the coordinate descent of ProximalSolvers (ProximalSolvers.CDSolver) with them. The second form copies o with the fields named in changes replaced, as in CDOptions(o; warm_start=false). The solvers also take these fields as keyword arguments, lasso(X, y, λ; max_iter=10_000, opt_tol=1e-8), which change the options passed to them.

• max_iter: maximum number of passes over the coordinates; reaching it without convergence gives a warning
• opt_tol: convergence tolerance on the largest change of a coordinate in a pass. With relative_tol, the change of coordinate k is measured as |Δxₖ| ⋅ coordinate_scale(f, k) and compared with opt_tol ⋅ response_scale(f), so that the tolerance does not depend on the units of the data; for the least squares loss these scales are ‖X[:, k]‖ / √n and ‖y‖ / √n. Otherwise |Δxₖ| is compared with opt_tol
• order: the order in which a pass visits the coordinates
  • :random: a new random order in every pass, drawn from the random number generator rng that the solvers take as a keyword argument, by default the global one (Random.default_rng())
  • :cyclic: the order 1, 2, …
  The cyclic order can be faster on some designs, but much slower when the columns are strongly correlated with each other. randomize = true and randomize = false are the same as order = :random and order = :cyclic
• warm_start: start from the supplied iterate; otherwise start from zero and follow a path of decreasing penalties
• num_steps: number of penalty levels on that path when warm_start is false
• relative_tol: see opt_tol; losses that do not define the scales, such as user-defined ones, use the absolute rule either way
• anderson: the memory K of Anderson acceleration, 0 for none (see coordinate_descent!); it needs order = :cyclic
• momentum: Nesterov momentum on the passes over the active set (see coordinate_descent!). The default, nothing, uses it where it is known to help: with order = :random, a SparseIterate and a built-in loss, and on the stacked design of multitask_lasso and MultitaskLassoPath, which ask for it on an AtomIterate. true asks for it with any loss and on an AtomIterate (ProximalSolvers 0.3.1), and needs order = :random; false turns it off
• gaptol: the relative duality gap at which each level of a path stops, and the last level of the path of warm_start = false (see ProximalSolvers' solve_path), for the losses and penalties that have a gap (least squares and the GLMs with the ℓ1 norm, its weighted and nonnegative forms, the elastic net, the group lasso and the sparse group lasso). nothing, the default, takes ProximalSolvers' default PATH_GAPTOL = 1e-4; 1e-3 is the faster setting; 0 stops on the changes of the coordinates with opt_tol, as before. A positive gaptol also stops a single solve (warm_start = true), which otherwise uses opt_tol. The coefficients are within about gaptol of the solution, relative to the largest, and within about 50 gaptol at worst

IterLassoOptions(; max_iter=20, opt_tol=1e-2, init_procedure=:Screening,
                   sinit=5, σinit=1.0, options_cd=CDOptions())
IterLassoOptions(o::IterLassoOptions; changes...)

Options for feasible_lasso!, which alternates a lasso fit with an update of the penalty loadings, and for scaled_lasso!, which uses only options_cd (since version 0.6.1 it solves the square-root lasso in one run of coordinate descent instead of alternating a lasso fit with an update of the noise level). The second form copies o with the fields named in changes replaced. The two solvers also take these fields as keyword arguments, as in feasible_lasso!(x, X, y, λ0; max_iter=50).

• max_iter: maximum number of outer iterations
• opt_tol: tolerance on the relative change of the noise level or loadings
• init_procedure: how the procedure is initialized
  • :Screening: least squares on the sinit columns of X most correlated with y
  • :InitStd: noise level σinit
  • :WarmStart: the iterate passed in
• sinit, σinit: used by :Screening and :InitStd
• options_cd: CDOptions for coordinate descent; in feasible_lasso! each lasso fit is warm started from the previous one, so its warm_start is not used there