Function reference

KronQuadratic(A, B, C = nothing)

The quadratic f(Θ) = ⟨Θ, AΘB⟩ / 2 + ⟨C, Θ⟩ of a p × q matrix Θ, for symmetric positive semidefinite A (p × p) and B (q × q), with C = 0 when omitted: vec(Θ)'(B ⊗ A) vec(Θ) / 2 + vec(C)'vec(Θ), without forming the pq × pq matrix B ⊗ A. Smooth, with gradient AΘB + C; size(f) is (p, q). A and B are stored as dense matrices; one that is symmetric only up to rounding, such as an inverse computed by inv, is symmetrized.

Θ may be a p × q matrix, the vector vec(Θ), a SparseIterate of either shape or, when p = q, a SymmetricSparseIterate. On a SymmetricSparseIterate the gradient is the one over symmetric matrices, (AΘB + BΘA + C + C')/2, so that a solver that steps on one triangle moves correctly. The value and the gradient use buffers that f owns and allocate nothing; do not share f between threads.

Examples:

• the direct difference (D-trace) loss tr(ΣxΔΣyΔ)/2 - tr(Δ(Σy - Σx)) of a symmetric Δ (CovSel) is KronQuadratic(Σx, Σy, Σx - Σy);
• the loss tr(ΣxΘΣyΘ')/2 - Θ[a, b] (DiffPrecTest), whose minimizer is the column (b - 1)p + a of (Σy ⊗ Σx)⁻¹, is KronQuadratic(Σx, Σy, -E), with E the indicator of the entry (a, b); DiffPrecTest's reduced loss on svec(Θ) is this function on a symmetric Θ, in coordinates that scale the off-diagonal entries by √2.

LinearFunction(c, d = 0)

The linear function f(x) = ⟨c, x⟩ + d, with gradient c. Its proximal operator is the shift x - γc.

QuadraticFunction(A, b = nothing, c = 0)

The quadratic f(x) = x'Ax / 2 + b'x + c for a symmetric positive semidefinite A, with b = 0 when omitted. Smooth, with gradient Ax + b, and proximable: the first proximal operator computes the eigendecomposition A = QΛQ' once, after which prox_{γf}(x) = Q (I + γΛ)⁻¹ Q'(x - γb) costs two products with Q for any γ.

SqrNormL2(λ = 1)
SqrNormL2(λ, y)

The squared Euclidean norm f(x) = (λ/2) ‖x - y‖², with y = 0 when omitted: the ridge penalty, and for λ = 1 the loss ‖x - y‖² / 2. Smooth, with gradient λ(x - y), separable, and proximable, with proximal operator (x + γλy) / (1 + γλ). For a matrix x, ‖⋅‖ is the Frobenius norm.

ZeroFunction()

The zero function f(x) = 0: smooth, with gradient zero, and proximable, with the identity as proximal operator. The penalty of an unpenalized problem.

LeastSquaresLoss(Y, X)

The least-squares loss f(β) = ‖Y - Xβ‖² / 2, smooth with gradient X'(Xβ - Y). Y may be a vector or an n × k matrix with one column per response; then β and the gradient are p × k matrices and ‖⋅‖ is the Frobenius norm, as in multi-response regression with NuclearNorm or the row groups of GroupNormL2. β may be a dense array, a SparseIterate or an AtomIterate, whose products cost O(n ⋅ nnz).

It is also proximable: the first proximal operator computes the thin SVD X = U Diagonal(s) V' and Diagonal(s) U'Y once, after which prox_{γf}(v) = (X'X + I/γ)⁻¹ (X'Y + v/γ) costs two products with V for any γ, also when p > n. ADMM changes γ freely at no cost.

GLMLoss(ℓ, y, X; weights = nothing, offset = nothing)

The loss f(β) = Σᵢ wᵢ ℓ(ηᵢ, yᵢ) of the linear predictor η = Xβ + offset, for a ProximalBase.ScalarLoss ℓ, with all wᵢ = 1 and no offset when omitted. Smooth, with gradient X' (wᵢ ℓ'(ηᵢ, yᵢ))ᵢ. LogisticLoss, PoissonLoss and HuberLoss construct the common ones.

f keeps X (field X) and a buffer for η (field η), so value_and_gradient! costs one product with X, one pass over the observations and one product with X', and allocates nothing. β may be a dense vector, a SparseIterate or an AtomIterate, whose products with X cost O(n ⋅ nnz). restrict keeps the columns S of X.

HuberLoss(y, X, δ; weights = nothing, offset = nothing)

The Huber regression loss Σᵢ wᵢ h_δ(ηᵢ - yᵢ), η = Xβ + offset, with h_δ(r) = r²/2 for |r| ≤ δ and δ(|r| - δ/2) beyond: a GLMLoss with ProximalBase.Huber.

LogisticLoss(y, X; weights = nothing, offset = nothing)

The logistic regression loss Σᵢ wᵢ (log(1 + exp(ηᵢ)) - yᵢηᵢ), η = Xβ + offset, for responses yᵢ ∈ {0, 1}: a GLMLoss with ProximalBase.Logistic.

PoissonLoss(y, X; weights = nothing, offset = nothing)

The Poisson regression loss with the log link Σᵢ wᵢ (exp(ηᵢ) - yᵢηᵢ), η = Xβ + offset: a GLMLoss with ProximalBase.Poisson. Its gradient is not globally Lipschitz; use a line search.

LogDetLoss(S, w = 1)

The Gaussian negative log-likelihood f(X) = w (tr(SX) - log det X) of a precision matrix X, for a symmetric S such as a sample covariance, and Inf when X is not positive definite. The loss of the graphical lasso.

It is smooth on the positive definite matrices, with gradient w(S - X⁻¹), which value_and_gradient! computes from one Cholesky factorization; the value never throws, so that a line search can reject a step that leaves the cone. It is also proximable: with ρ = 1/(γw) and the eigendecomposition ρX - S = U Diagonal(t) U', prox_{γf}(X) = U Diagonal(μ) U' with μᵢ = (tᵢ + √(tᵢ² + 4ρ)) / (2ρ), and value_and_prox! returns the value from the μᵢ. Both use workspaces that f owns, sized by S, and allocate nothing; do not share f between threads. The element type of X must be that of S, Float32 or Float64.

ElasticNet(λ1, λ2)
ElasticNet(λ1, λ2, w)

The elastic-net penalty f(x) = λ1 Σᵢ wᵢ |xᵢ| + (λ2/2) ‖x‖², with all wᵢ = 1 when w is omitted. Its proximal operator is shrink(x, γλ1w) / (1 + γλ2).

IndBox(l, u)

The indicator of the box {x : l ≤ x ≤ u}: 0 inside and Inf outside. The bounds are scalars or arrays of the size of x, and may be infinite. Its proximal operator is the projection clamp.(x, l, u).

IndNonnegative()

The indicator of the nonnegative orthant {x : x ≥ 0}. Its proximal operator is max.(x, 0).

MCP(λ, a = 3)

The minimax concave penalty of Zhang (2010), nonconvex: fᵢ(t) = λ|t| - t²/(2a) for |t| ≤ aλ and aλ²/2 beyond, for a > 0. It matches the ℓ1 penalty at zero and stops penalizing large entries. For a step γ < a its proximal operator is firm thresholding; for larger steps it is the global minimizer of the nonconvex problem.

NormL0(λ = 1)

The ℓ0 penalty f(x) = λ ⋅ #{i : xᵢ ≠ 0}, nonconvex. Its proximal operator is hard thresholding: xᵢ is kept when |xᵢ| > √(2γλ) and set to zero otherwise.

NormL1(λ0 = 1)
NormL1(λ0, w)

The weighted ℓ1 norm f(x) = λ0 Σᵢ wᵢ |xᵢ|, with all wᵢ = 1 when w is omitted. Its proximal operator is entry-wise soft thresholding (shrink). Separable, so prox_coordinate! updates any coordinate, and on a sparse iterate the value and the proximal operator touch only the stored entries. On a SymmetricSparseIterate the weights are a p × p matrix (or anything indexed by the linear indices of one), as for a dense matrix.

NormL1Nonneg(λ0 = 1)
NormL1Nonneg(λ0, w)

The weighted ℓ1 norm restricted to the nonnegative orthant: f(x) = λ0 Σᵢ wᵢ xᵢ for x ≥ 0, and Inf otherwise, the penalty of the nonnegative lasso. Its proximal operator is max(x - γλ0w, 0).

SCAD(λ, a = 3.7)

The smoothly clipped absolute deviation penalty of Fan and Li (2001), nonconvex: fᵢ(t) = λ|t| for |t| ≤ λ, (2aλ|t| - t² - λ²) / (2(a - 1)) for λ < |t| ≤ aλ, and λ²(a + 1)/2 beyond, for a > 2. For a step γ < a - 1 its proximal operator is the SCAD thresholding rule; for larger steps it is the global minimizer of the nonconvex problem.

IndBallLinf(r)

The indicator of the ℓ∞ ball {x : maxᵢ |xᵢ| ≤ r}, the box IndBox(-r, r).

CheckLoss(τ)
CheckLoss(τ, w)

The check loss of quantile regression, f(x) = Σᵢ wᵢ ρ_τ(xᵢ) with ρ_τ(u) = max(τu, (τ - 1)u) for τ ∈ (0, 1), and all wᵢ = 1 when w is omitted. Its proximal operator moves an entry above γwᵢτ down by γwᵢτ, one below -γwᵢ(1 - τ) up by γwᵢ(1 - τ), and sets those in between to zero.

It is a function of residuals. Since ρ_τ(-u) = ρ_{1-τ}(u), the loss Σᵢ ρ_τ(yᵢ - zᵢ) of predictions z of responses y is Translate(CheckLoss(1 - τ), y), which a splitting method such as ADMM applies to z = Xβ. Its conjugate is the indicator of the box [wᵢ(τ - 1), wᵢτ], the dual feasible set of quantile regression. For coordinate descent and Newton steps, ProximalBase.SmoothCheck and ProximalBase.ConvolutionCheck are smooth versions of the loss of Xβ for GLMLoss.

EpsInsensitiveLoss(ε)
EpsInsensitiveLoss(ε, w)

The ε-insensitive loss of support vector regression, f(x) = Σᵢ wᵢ max(0, |xᵢ| - ε) for ε ≥ 0, with all wᵢ = 1 when w is omitted: residuals within ε cost nothing. Its proximal operator keeps an entry with |v| ≤ ε, puts one with ε < |v| ≤ ε + γwᵢ at ±ε, and moves a larger one towards zero by γwᵢ. As for CheckLoss, the loss of predictions z of y is Translate(EpsInsensitiveLoss(ε), y). Its conjugate is ε Σᵢ |sᵢ| on the box |sᵢ| ≤ wᵢ.

HingeLoss()
HingeLoss(y)
HingeLoss(y, w)

The hinge loss of support vector machines, f(x) = Σᵢ wᵢ max(0, 1 - yᵢxᵢ), for nonzero labels yᵢ, usually ±1, and weights wᵢ, all 1 when omitted. Without labels, x holds the margins. For the label 1, its proximal operator moves an entry below 1 - γwᵢ up by γwᵢ, puts one in between at 1, and keeps one above 1; a label yᵢ acts on the margin yᵢxᵢ, with the step γwᵢyᵢ². Its conjugate is Σᵢ sᵢ/yᵢ on the box -wᵢ ≤ sᵢ/yᵢ ≤ 0, and Inf off it.

Berhu(λ, M)
Berhu(λ, M, w)

The reverse Huber penalty of Owen (2007): fᵢ(t) = λwᵢ|t| for |t| ≤ M and λwᵢ(t² + M²)/(2M) beyond, for M > 0, with all wᵢ = 1 when w is omitted. It is the ℓ1 norm near zero, for sparsity, and a ridge penalty for large entries, which groups correlated predictors; convex, with the slope λwᵢ at ±M on both sides. With c = γλwᵢ, its proximal operator is

0                      for |v| ≤ c,
v - c sign(v)          for c < |v| ≤ M + c,
vM / (M + c)           for |v| > M + c.

At zero its subdifferential is that of the ℓ1 norm, so kkt_violation is |gᵢ| - λwᵢ.

KLDivergence(q; λ = 1)

The Kullback–Leibler divergence f(x) = λ Σᵢ (xᵢ log(xᵢ/qᵢ) - xᵢ + qᵢ) from a positive reference q (a scalar, or an array of the size of x) on x ≥ 0, and Inf for a negative entry, for λ > 0. Separable and convex, and zero at x = q. It is the Poisson deviance of q as a function of the mean x. With μ = γλ, its proximal operator is

μ ω(v/μ + log(qᵢ/μ)),

with the Wright omega function ω, as for NegEntropy.

NegEntropy(λ = 1)

The negative entropy f(x) = λ Σᵢ xᵢ log xᵢ on x ≥ 0 (with 0 log 0 = 0), and Inf for a negative entry, for λ > 0. Separable and convex. With μ = γλ, its proximal operator is

μ ω(v/μ - 1 - log μ),

where ω(u) = W₀(eᵘ) is the Wright omega function, computed without forming e^(v/μ), so the proximal operator stays accurate for a small step. The result is positive, so the function does not preserve zeros.

NegLog(λ = 1)

The log barrier f(x) = -λ Σᵢ log xᵢ on x > 0, and Inf elsewhere, for λ > 0. Separable and convex. With μ = γλ, its proximal operator is the positive root (v + √(v² + 4μ))/2 of z² - vz - μ = 0, computed without cancellation for negative v, the root that the proximal operator of LogDetLoss takes for each eigenvalue.

GroupNormL2(λ0, groups)
GroupNormL2(λ0, groups, w)

The group-lasso penalty f(x) = λ0 Σₖ wₖ ‖x[Gₖ]‖₂ over disjoint groups Gₖ = groups[k] of (linear) indices of x, of any sizes, such as 1:5 or [2, 7, 9]. The weights default to wₖ = √|Gₖ|, so that groups of one index give NormL1(λ0). Indices in no group are not penalized. The proximal operator shrinks each group towards zero in Euclidean norm.

The rows of a p × k matrix are the groups [i:p:p*k for i in 1:p], which gives the multi-task group lasso with a matrix response in LeastSquaresLoss.

On an AtomIterate, whose coordinates are groups, prox_coordinate! updates a group and kkt_violation measures one. AtomIterate(x, f) builds an iterate that shares the groups of f, which must then partition 1:length(x) (an index in no group goes in a group of its own with weight 0); with shared groups, the value and the proximal operator touch only the active groups. On ordinary arrays, active-set solvers can use restrict and kkt_violation.

Example

f = GroupNormL2(0.5, [1:3, 4:5, [6, 8]])    # x[7] is not penalized
prox(f, randn(8))

GroupNormLinf(λ0, groups)
GroupNormLinf(λ0, groups, w)

The ℓ1/ℓ∞ penalty f(x) = λ0 Σₖ wₖ maxᵢ∈Gₖ |xᵢ| over disjoint groups, as GroupNormL2 but with the ℓ∞ norm of each group, for simultaneous selection in multi-task problems. The weights default to wₖ = 1. The proximal operator clips each group, by Moreau's identity with the projection onto an ℓ1 ball, and its optimality test (kkt_violation) uses the dual norm, ‖g[Gₖ]‖₁ - λ0 wₖ.

NormL1Fused(λ1, λ2)

The penalty f(x1, x2) = Σᵢ λ1 (|x1ᵢ| + |x2ᵢ|) + λ2 |x1ᵢ - x2ᵢ| of a pair of arrays of the same size, given as a tuple (x1, x2), as in the joint estimation of two sparse models that should be close. Its proximal operator acts entry by entry on the pairs (x1ᵢ, x2ᵢ); prox(f, (u, v), γ) gives it for two scalars.

On a vector v of length K it is λ1 Σₖ |vₖ| + λ2 Σₖ<ₗ |vₖ - vₗ|, which fuses all pairs of entries; for K = 2 this is the function of one pair. As the fiber function of FiberSum, FiberSum(NormL1Fused(λ1, λ2)) on an n × n × K array is the penalty of the fused graphical lasso of K classes (Danaher, Wang and Witten, 2014). With the entries sorted, v₍₁₎ ≤ … ≤ v₍K₎, the fusion term is Σᵢ (2i - K - 1) v₍ᵢ₎, linear on the vectors sorted the same way, and its proximal operator keeps the order of v. So the proximal operator sorts v, subtracts γλ2 (2i - K - 1) from the i-th smallest entry, pools adjacent violators into a nondecreasing sequence, restores the order, and soft-thresholds by γλ1: exact, in O(K log K), in buffers that the function owns.

NormL2(λ = 1)

The Euclidean norm f(x) = λ ‖x‖₂ (the Frobenius norm for a matrix). Its proximal operator shrinks x towards zero: max(1 - γλ/‖x‖, 0) x. For the group lasso, a sum of such norms over groups of entries, use GroupNormL2.

NormLinf(λ = 1)

The ℓ∞ norm f(x) = λ maxᵢ |xᵢ|. Its proximal operator clips the entries at a level τ, sign(xᵢ) min(|xᵢ|, τ): by Moreau's identity it is x minus the projection onto the ℓ1 ball of radius γλ, which takes O(n) expected time (Condat's algorithm).

FiberSum(f; diagonal = f, symmetric = false)

The function Σᵢⱼ f(x[i, j, :]) of an m × n × K array x: a vector function f of each fiber across the K slices, as in the penalties of K related precision matrices that make them share a sparsity pattern or be close. diagonal is the function of the fibers x[i, i, :] of square slices: nothing leaves them unpenalized, and another function gives them their own rule. The proximal operator applies the proximal operator of f to each fiber, on views, and allocates nothing beyond what that of f does. On a matrix, the fibers are its rows.

With symmetric = true the slices are taken to be symmetric, as precision matrices are: only the fibers with i ≥ j are read, and the proximal operator copies each result to the fiber (j, i), which halves the work and keeps the result exactly symmetric. The value counts each off-diagonal fiber twice.

With x[:, :, k] the precision matrix of class or time point k:

and diagonal = nothing leaves the diagonal free, diagonal = NormL1Fused(0, λ2) keeps only the fusion on it.

SliceSum(f₁, …, f_K)
SliceSum(fs::AbstractVector)

The function Σₖ fₖ(xₖ) of an array x whose slices along the last dimension are x₁, …, x_K: the matrices x[:, :, k] of an m × n × K array, or the columns of a matrix. The value, the gradient and the proximal operator act slice by slice, on views, and allocate nothing beyond what the fₖ do. SliceSum has the traits that all the fₖ share. The fₖ are given as arguments, or as a vector when there are many of one type.

With fₖ = LogDetLoss(Sₖ, nₖ/n) it is the smooth part of the joint graphical lasso of K classes, whose penalty on the entries across the classes is a FiberSum.

IndBallRank(r; partial = false)

The indicator of the matrices of rank at most r, nonconvex. Its proximal operator is the truncated singular value decomposition, which keeps the r largest singular triplets.

With partial = true the proximal operator computes only the leading singular triplets or eigenpairs, by randomized subspace iteration started from those of the previous call, in O(p²k) instead of O(p³) for k of them. It accepts the result when the residuals of the pairs kept are below eps^(3/4) times the largest value, and otherwise, or when the pairs kept would be more than about a quarter of the dimension, falls back to the full decomposition. This pays when the rank of the result is small and the proximal operator is called on a sequence of nearby matrices, as in a solver; see Partial decompositions.

The function caches a LAPACK workspace for the size and element type (Float32 or Float64) of its input, so that repeated calls allocate nothing; do not share one function object between threads.

IndPSD(; partial = false)
IndPSD(; rank = r, partial = false)

The indicator of the positive semidefinite matrices, and with rank of those of rank at most r, which is nonconvex. Its proximal operator is the projection of a symmetric matrix (its upper triangle is read) onto the set: the eigendecomposition with the negative eigenvalues set to zero, and for a rank r only the r largest eigenpairs, which LAPACK computes without the others. The result is rebuilt from the eigenvectors that remain, in O(n²k) for k of them.

With partial = true the proximal operator computes only the leading singular triplets or eigenpairs, by randomized subspace iteration started from those of the previous call, in O(p²k) instead of O(p³) for k of them. It accepts the result when the residuals of the pairs kept are below eps^(3/4) times the largest value, and otherwise, or when the pairs kept would be more than about a quarter of the dimension, falls back to the full decomposition. This pays when the rank of the result is small and the proximal operator is called on a sequence of nearby matrices, as in a solver; see Partial decompositions. Without a rank, the partial projection pays when few eigenvalues are positive.

The function caches a LAPACK workspace for the size and element type (Float32 or Float64) of its input, so that repeated calls allocate nothing; do not share one function object between threads.

NuclearNorm(λ = 1; partial = false)

The nuclear norm f(X) = λ ‖X‖_*, the sum of the singular values of X. Its proximal operator soft-thresholds the singular values, and value_and_prox! returns the value from the thresholded singular values without a second decomposition.

With partial = true the proximal operator computes only the leading singular triplets or eigenpairs, by randomized subspace iteration started from those of the previous call, in O(p²k) instead of O(p³) for k of them. It accepts the result when the residuals of the pairs kept are below eps^(3/4) times the largest value, and otherwise, or when the pairs kept would be more than about a quarter of the dimension, falls back to the full decomposition. This pays when the rank of the result is small and the proximal operator is called on a sequence of nearby matrices, as in a solver; see Partial decompositions. The value f(X) alone needs every singular value and always takes the full decomposition.

The function caches a LAPACK workspace for the size and element type (Float32 or Float64) of its input, so that repeated calls allocate nothing; do not share one function object between threads.

TracePSD(λ = 1; partial = false)

The function f(X) = λ tr(X) on the positive semidefinite matrices, Inf elsewhere: the nuclear norm restricted to them, as in the low-rank part of latent-variable graphical models. Its proximal operator shifts the eigenvalues of a symmetric matrix down by γλ and sets the negative ones to zero.

With partial = true the proximal operator computes only the leading singular triplets or eigenpairs, by randomized subspace iteration started from those of the previous call, in O(p²k) instead of O(p³) for k of them. It accepts the result when the residuals of the pairs kept are below eps^(3/4) times the largest value, and otherwise, or when the pairs kept would be more than about a quarter of the dimension, falls back to the full decomposition. This pays when the rank of the result is small and the proximal operator is called on a sequence of nearby matrices, as in a solver; see Partial decompositions.

The function caches a LAPACK workspace for the size and element type (Float32 or Float64) of its input, so that repeated calls allocate nothing; do not share one function object between threads.

IndBallL0(s; diagonal = true)

The indicator of the arrays with at most s nonzero entries, nonconvex. Its proximal operator keeps the s entries of largest magnitude (exactly s when there are ties) and sets the others to zero: a selection in O(n) expected time (Floyd and Rivest's algorithm), without allocating, and in O(nnz) on a SparseIterate.

For a symmetric matrix, a Symmetric or a SymmetricSparseIterate, s counts the entries of one triangle, so that an entry and its transpose are kept or dropped together, and the result is symmetric; the output may be a dense matrix or a SymmetricSparseIterate. The result is the projection in the Frobenius norm: an off-diagonal entry stands for a pair, whose removal costs twice its square, so it competes with the diagonal by √2|xᵢⱼ|. With diagonal = false, the diagonal of a matrix does not count and is kept, as for the off-diagonal sparsity of a precision matrix:

prox!(S, IndBallL0(s; diagonal = false), Symmetric(Θ, :L))   # S::SymmetricSparseIterate

IndBallL1(r = 1)

The indicator of the ℓ1 ball {x : Σᵢ |xᵢ| ≤ r}. Its proximal operator is the projection, soft thresholding at the level that puts the result on the boundary, found in O(n) expected time by Condat's algorithm.

IndBallL2(r = 1)

The indicator of the Euclidean ball {x : ‖x‖₂ ≤ r} (the Frobenius ball for matrices). Its proximal operator scales x onto the ball when it is outside.

IndHalfspace(a, b)

The indicator of the halfspace {x : ⟨a, x⟩ ≤ b}, a ≠ 0. Its proximal operator is the projection, x - (max(⟨a, x⟩ - b, 0) / ‖a‖²) a. A scalar a stands for an array of as of the size of x, as in IndHyperplane.

IndHyperplane(a, b)

The indicator of the hyperplane {x : ⟨a, x⟩ = b}, a ≠ 0. Its proximal operator is the projection x - ((⟨a, x⟩ - b) / ‖a‖²) a. A scalar a stands for an array of as of the size of x: IndHyperplane(1, 0) is the set of arrays whose entries sum to zero, of any size.

IndSimplex(r = 1; equality = true)

The indicator of the simplex {x ≥ 0 : Σᵢ xᵢ = r}, or with equality = false of {x ≥ 0 : Σᵢ xᵢ ≤ r}. Its proximal operator is the projection, max.(x .- τ, 0) for the level τ found in O(n) expected time by Condat's algorithm.

IndSphereL2(r = 1)

The indicator of the Euclidean sphere {x : ‖x‖₂ = r} (the Frobenius sphere for matrices), for r > 0, nonconvex. Its proximal operator scales x onto the sphere, r x / ‖x‖₂; every point of the sphere is nearest to x = 0, which goes to r times the first unit vector. With a norm, Constrained(f, IndSphereL2(r)) is the step of sparse principal component analysis.

FusedLasso(λ1, λ2)

The penalty f(x) = λ1 Σᵢ |xᵢ| + λ2 Σᵢ |xᵢ₊₁ - xᵢ| of the fused lasso signal approximator, on a vector. Its proximal operator soft-thresholds the total-variation denoising of x by γλ1 (Friedman, Hastie, Höfling and Tibshirani, 2007).

IndIsotonic(; decreasing = false)

The indicator of the nondecreasing vectors {x : x₁ ≤ x₂ ≤ … ≤ xₙ}, or with decreasing = true of the nonincreasing ones. Its proximal operator, the projection, is isotonic regression, computed by pooling adjacent violators in O(n) time, in buffers that the function owns.

SLOPE(λ)

The sorted ℓ1 norm f(x) = Σᵢ λᵢ |x|₍ᵢ₎, where |x|₍₁₎ ≥ |x|₍₂₎ ≥ … are the magnitudes of the entries in decreasing order and λ₁ ≥ λ₂ ≥ … ≥ 0: the penalty of SLOPE (Bogdan, van den Berg, Sabatti, Su and Candès, 2015), which controls the false discovery rate of the selected variables with the Benjamini–Hochberg sequence of λ. With all λᵢ equal it is the ℓ1 norm, and with λ = (λ₁, 0, …, 0) the ℓ∞ norm. Its proximal operator sorts the magnitudes, subtracts γλ, pools adjacent violators into a nonincreasing sequence, clips it at zero, and restores the order and the signs: O(n log n), in buffers that the function owns.

TotalVariation1D(λ = 1)

The total variation f(x) = λ Σᵢ |xᵢ₊₁ - xᵢ| of a vector. Its proximal operator, total-variation denoising, is computed exactly and in place by Condat's direct algorithm (2013), in O(n) time on most inputs and O(n²) in the worst case.

CappedL1(λ, a)
CappedL1(λ, a, w)

The capped ℓ1 penalty of Zhang (2010), nonconvex: fᵢ(t) = λ wᵢ min(|t|, a) for a > 0, with all wᵢ = 1 when w is omitted. It is the ℓ1 penalty up to a and constant beyond, so that entries larger than a are not shrunk. Its proximal operator is the global minimizer, the better of soft thresholding clipped at a and the identity clipped below at a. At zero its subdifferential is that of the ℓ1 norm, so kkt_violation is |gᵢ| - λwᵢ.

GroupMCP(λ, groups; a = 3)
GroupMCP(λ, groups, w; a = 3)

The group minimax concave penalty of Huang, Breheny and Ma (2012), nonconvex: f(x) = Σₖ pₖ(‖x[Gₖ]‖₂), where pₖ is the MCP with level λwₖ, over disjoint groups Gₖ = groups[k] of (linear) indices. Groups are selected as a whole, as by GroupNormL2, but a group of large norm is not shrunk. The weights default to wₖ = √|Gₖ|, and indices in no group are not penalized.

Its proximal operator is radial: the scalar proximal operator of pₖ at ‖v[Gₖ]‖ rescales v[Gₖ], which gives the global minimizer for every step. At zero its subdifferential is that of the group lasso, so kkt_violation is ‖g[Gₖ]‖₂ - λwₖ. On an AtomIterate with its groups, prox_coordinate! updates a group, the block step of coordinate descent.

GroupSCAD(λ, groups; a = 3.7)
GroupSCAD(λ, groups, w; a = 3.7)

The group SCAD penalty, nonconvex: f(x) = Σₖ pₖ(‖x[Gₖ]‖₂), where pₖ is the SCAD penalty with level λwₖ, over disjoint groups, as GroupMCP is for the MCP; a > 2. Its proximal operator is radial and exact for every step.

IndBallGroupL0(s, groups; k = nothing)

The indicator of the vectors with at most s nonzero groups, over disjoint groups Gₖ = groups[k] of (linear) indices, nonconvex; indices in no group are free. Its proximal operator keeps the s groups of largest Euclidean norm (exactly s when there are ties, the earlier groups first) and zeroes the others, by a selection in O(n) expected time, without allocating. Projected gradient descent with it is iterative hard thresholding of groups, as in the group IHT of MendelIHT (Chu et al., 2020). On an AtomIterate with its groups, only the active groups are visited, and the dropped ones are deactivated.

With k, each group also has at most k nonzero entries: the doubly sparse set of group IHT with sparsity within the groups. The projection ranks the groups by the norm of their k largest entries, keeps the s best, and within each the k largest entries (exactly k when there are ties, the earlier entries first).

LogPenalty(λ, ε)
LogPenalty(λ, ε, w)

The log penalty fᵢ(t) = λ wᵢ log(1 + |t|/ε) for ε > 0, nonconvex, with all wᵢ = 1 when w is omitted: the penalty whose linearization gives the reweighted ℓ1 minimization of Candès, Wakin and Boyd (2008). It grows only logarithmically, so it shrinks large entries little, and it approaches the ℓ0 penalty as ε → 0. Its proximal operator is the global minimizer: zero or the larger root of a quadratic, whichever is better. Its slope at zero is λwᵢ/ε, so kkt_violation is |gᵢ| - λwᵢ/ε.

NormLq(λ, q)
NormLq(λ, q, w)

The ℓq penalty fᵢ(t) = λ wᵢ |t|^q for 0 < q < 1, nonconvex, with all wᵢ = 1 when w is omitted. It lies between the ℓ0 penalty (NormL0, the limit q → 0) and the ℓ1 norm (q = 1), and shrinks large entries less than the ℓ1 norm. Its proximal operator is the global minimizer. With μ = γλwᵢ it is zero for |v| up to the threshold

τ = (2 - q)/(2 - 2q) ⋅ (2μ(1 - q))^(1/(2 - q)),

where zero and the nonzero candidate tie, and above it the root z > 0 of z + μq z^(q-1) = |v| that exceeds (2μ(1 - q))^(1/(2 - q)), with the sign of v. The root has a closed form for q = 1/2 (half thresholding: Xu, Chang, Xu and Zhang, 2012) and q = 2/3 (Cao, Sun and Xu, 2013), and comes from Newton's method otherwise.

Its slope at zero is infinite, so zero is a local minimizer of every problem h + f, and kkt_violation is -Inf: an active-set method that grows the support from zero by that test never adds a coordinate. Start such a method from a nonzero point, or use the proximal operator of the whole vector, as proximal gradient methods do.

SingularValue(f)

The function F(X) = f(σ(X)) of the singular values σ(X) of a matrix, for a vector function f with a proximal operator that is absolutely symmetric: invariant under permutations and sign changes of its argument. Its proximal operator acts on the singular values,

prox_{γF}(X) = U Diagonal(prox_{γf}(σ)) Vᵀ,   X = U Diagonal(σ) Vᵀ (thin SVD),

(Lewis, 1995), and value_and_prox! returns f at the new singular values, without a second decomposition. For a nonconvex separable f whose scalar proximal operator is nondecreasing on [0, ∞), such as MCP, SCAD or NormL0, the same formula gives a global minimizer (generalized singular value thresholding; Lu, Zhu, Xu, Yan and Lin, 2015).

The Ky Fan norm is the support function of {w : 0 ≤ wᵢ ≤ λ, Σᵢ wᵢ ≤ λk} at the singular values. It needs the halfspace: the conjugate of the capped simplex, with the hyperplane Σᵢ wᵢ = k, is the sum of the k largest entries, which is not absolutely symmetric.

isconvex and the indicator property come from f. In the value F(X), singular values within a relative √ε of the largest count as zero, as they do for IndBallRank: the value at a proximal point of low rank then agrees with the one value_and_prox! returns, also for a discontinuous f such as NormL0. The value of Conjugate(SingularValue(f)) is f*(σ(Y)), known when that of Conjugate(f) is, as for the nuclear-norm ball and the spectral norm, which are conjugate to each other.

The function caches a LAPACK workspace for the size and element type (Float32 or Float64) of its input, so that repeated calls allocate nothing; do not share one function object between threads.

IndBanded(k, p)
IndBanded(k, (m, n))
IndBanded((kl, ku), (m, n))

The indicator of the banded m × n matrices, p × p in the first form: those whose entries more than kl below the diagonal or more than ku above it are zero, with kl = ku = k for a single bandwidth. k = 0 gives the diagonal matrices. Its proximal operator zeroes the entries outside the band, which for a sample covariance is banding (Bickel and Levina, 2008).

The size is part of the function, so that the entry with linear index i has a known row and column: like IndPattern, it is separable and preserves zeros, without storing a pattern, and it works on dense and Symmetric matrices of that size and on the sparse iterates. On a SymmetricSparseIterate the band must be symmetric, kl = ku.

IndPattern(P)

The indicator of the arrays supported on a pattern: those whose entries outside the pattern are zero. P is an array of Bools, whose true entries form the pattern, or any other array, whose nonzero entries do, such as a SparseMatrixCSC, a SymmetricSparseIterate or an estimate whose support a refit keeps. The pattern is stored as a BitArray; a BitArray argument is used as it is, not copied.

Its proximal operator zeroes the entries outside the pattern. It is separable and preserves zeros, so prox_coordinate! updates any coordinate, and on a sparse iterate the value and the proximal operator touch only the stored entries. With LogDetLoss in ADMM it gives the Gaussian maximum likelihood estimate on a known graph. On a SymmetricSparseIterate the pattern must be symmetric.

IndToeplitz(; symmetric = false)

The indicator of the Toeplitz matrices, which are constant along each diagonal: Xᵢⱼ = c_{j-i}. Its proximal operator replaces each diagonal by its mean. With symmetric = true the set is that of the symmetric Toeplitz matrices, which are square, and the projection averages the diagonals d and -d together. These are the covariance matrices of a stationary series. The function is not separable. The projection takes two passes over each diagonal, works in place, and allocates nothing. Its output may be a SymmetricSparseIterate when the result is symmetric: with symmetric = true, or from a symmetric input.

IndAffine(A, b)
IndAffine(A, b; μ)

The indicator of the affine set {x : Ax = b} for an m × n matrix A (dense, sparse, or anything with mul!) and a vector b of length m. Its proximal operator is the projection

x - Aᵀ(AAᵀ)⁻¹(Ax - b),

from a factorization cached at construction:

• full row rank: a Cholesky factorization of the m × m matrix AAᵀ;
• otherwise: a thin SVD A = U Diagonal(s) Vᵀ of the rank r of A, with the projection x - V(Vᵀx - Diagonal(s)⁻¹Uᵀb). b must be in the range of A, or the set is empty and the constructor throws;
• with μ: A is a tight frame, AAᵀ = μI, and nothing is factorized.

The rank-deficient path is taken when a pivot of the Cholesky factorization is below eps^(1/4) times the largest: rounding leaves a dependent row a pivot of about √eps times the largest, and the factor of AAᵀ squares the condition number of A. A call costs two products with A (or with V) and allocates nothing. The value is 0 when ‖Ax - b‖ is within a relative √eps of the scale of the problem, and Inf otherwise.

lla_weight(f, v, i)

The slope pᵢ′(|v|) of the nonconvex penalty f(x) = Σᵢ pᵢ(|xᵢ|) at |v|, taken from the right: the weight of entry i in the local linear approximation (LLA) of f at a point whose entry i is v. i is a linear index, as in value_entry, and the entry weights of f are included. For GroupMCP and GroupSCAD, i is a group and v its Euclidean norm. lla_weights! gives the weights of a whole point.

At zero the slope is pᵢ′(0+), so the first stage of LLA from x = 0 is the lasso. For the capped ℓ1 penalty at |v| = a the weight is 0, the rule of the multi-stage convex relaxation of T. Zhang (2010): λwᵢ for |v| < a and 0 otherwise. The slopes of the MCP and SCAD are continuous away from zero.

NormLq has an infinite slope at zero, so every stage would keep the zero entries at zero; it throws an ArgumentError, as do functions that are not of this form.

lla_weights!(w, f, x)
lla_weights(f, x)

Store in w the weights of the local linear approximation (LLA) of the nonconvex penalty f at x, and return w. The next stage of LLA, the weighted lasso that majorizes f at x up to a constant, has the penalty

• NormL1(1, w) for a separable f (CappedL1, LogPenalty, MCP, SCAD): w has the size of x, and w[i] =lla_weight(f, x[i], i), the slope pᵢ′(|xᵢ|) from the right;
• GroupNormL2(1, f.groups, w) for GroupMCP and GroupSCAD: w has one entry per group, w[k] = pₖ′(‖x[Gₖ]‖₂).

The weights include the level λ and the entry or group weights of f; to keep the level apart, divide them by λ and use NormL1(λ, w). x may be a SparseIterate, a SymmetricSparseIterate (with a p × p array w) or an AtomIterate. lla_weights! allocates nothing; lla_weights(f, x) allocates w.

Each LLA stage is convex, so its minimizer, λmax, the strong rule and safe screening are those of a weighted lasso. Alternating stages and weights is a majorization-minimization algorithm: each stage solved exactly decreases the loss plus f.

Example

f = MCP(0.5, 3)
w = lla_weights(f, [0.0, 0.6, 2.0])     # [0.5, 0.3, 0.0]
NormL1(1, w)                            # the penalty of the next stage

lla_weights!(w, f, x)
lla_weights(f, x)

Store in w the weights of the local linear approximation (LLA) of the nonconvex penalty f at x, and return w. The next stage of LLA, the weighted lasso that majorizes f at x up to a constant, has the penalty

• NormL1(1, w) for a separable f (CappedL1, LogPenalty, MCP, SCAD): w has the size of x, and w[i] =lla_weight(f, x[i], i), the slope pᵢ′(|xᵢ|) from the right;
• GroupNormL2(1, f.groups, w) for GroupMCP and GroupSCAD: w has one entry per group, w[k] = pₖ′(‖x[Gₖ]‖₂).

The weights include the level λ and the entry or group weights of f; to keep the level apart, divide them by λ and use NormL1(λ, w). x may be a SparseIterate, a SymmetricSparseIterate (with a p × p array w) or an AtomIterate. lla_weights! allocates nothing; lla_weights(f, x) allocates w.

Each LLA stage is convex, so its minimizer, λmax, the strong rule and safe screening are those of a weighted lasso. Alternating stages and weights is a majorization-minimization algorithm: each stage solved exactly decreases the loss plus f.

Example

f = MCP(0.5, 3)
w = lla_weights(f, [0.0, 0.6, 2.0])     # [0.5, 0.3, 0.0]
NormL1(1, w)                            # the penalty of the next stage

SparseGroupNormL2(λ1, λ2, groups)
SparseGroupNormL2(λ1, λ2, groups, w)

The sparse-group-lasso penalty

f(x) = λ1 Σᵢ |xᵢ| + λ2 Σₖ wₖ ‖x[Gₖ]‖₂

over disjoint groups Gₖ = groups[k] of (linear) indices of x, with the weights wₖ = √|Gₖ| by default, as in GroupNormL2. The group term selects groups and the ℓ1 term zeroes entries inside a selected group (Simon, Friedman, Hastie and Tibshirani, 2013); λ1 = 0 gives the group lasso. Indices in no group are not penalized.

Its proximal operator is exact: soft thresholding at γλ1, then the shrinkage of each group at γλ2wₖ, in one pass over a group, which also gives the value. With one group covering x, SparseGroupNormL2(λ1, λ2, [1:n], [1.0]) is the sparse group lasso of a single vector. As the fiber function of a FiberSum, FiberSum(SparseGroupNormL2(λ1, λ2, [1:K], [1.0])) is the penalty λ1 Σᵢⱼₖ |xᵢⱼₖ| + λ2 Σᵢⱼ ‖x[i, j, :]‖₂ of the group graphical lasso of K precision matrices (Danaher, Wang and Witten, 2014).

As for GroupNormL2, on an AtomIterate that shares its groups (AtomIterate(x, f)), prox_coordinate! takes the block step on a group and kkt_violation measures a group at zero, which is optimal when ‖soft(g[Gₖ], λ1)‖₂ ≤ λ2wₖ. Inside a nonzero group, an entry at zero is optimal when |gᵢ| ≤ λ1.

ProximalBase.ScalarLoss

Abstract supertype of the losses ℓ(η, y) of one observation y as a function of its linear predictor η, the building blocks of GLMLoss. A loss is callable, ℓ(η, y), and implements ProximalBase.deriv, ProximalBase.deriv2 and ProximalBase.curvature, and ProximalBase.value_and_deriv when the value and the derivative share work. The losses are ProximalBase.Squared, ProximalBase.Logistic, ProximalBase.Poisson, ProximalBase.Huber, ProximalBase.SquaredHinge, ProximalBase.SmoothCheck and ProximalBase.ConvolutionCheck; they are not exported, as names such as Poisson belong to distributions.

ProximalBase.Squared()

The squared loss ℓ(η, y) = (η - y)² / 2.

ProximalBase.Logistic()

The logistic loss ℓ(η, y) = log(1 + exp(η)) - yη, the negative log-likelihood of logistic regression with y ∈ {0, 1} (or a proportion), computed without overflow.

ProximalBase.Poisson()

The Poisson loss with the log link, ℓ(η, y) = exp(η) - yη, the negative log-likelihood up to a term in y. Its curvature is unbounded.

ProximalBase.Huber(δ)

The Huber loss of the residual r = η - y: r²/2 for |r| ≤ δ and δ(|r| - δ/2) beyond, for robust regression.

ProximalBase.SquaredHinge()

The squared hinge loss ℓ(η, y) = max(0, 1 - yη)² for a label y ∈ {-1, 1}, the loss of the L2-loss support vector machine. Smooth, with curvature 2, so coordinate descent and gradient methods minimize a GLMLoss with it directly.

ProximalBase.SmoothCheck(τ, h)

The check loss of quantile regression smoothed at the scale h > 0: the Moreau envelope min_s ρ_τ(s) + (s - u)²/(2h) of ρ_τ(u) = max(τu, (τ - 1)u) at the residual u = y - η. It is τu - hτ²/2 for u ≥ hτ, (τ - 1)u - h(1 - τ)²/2 for u ≤ -h(1 - τ), and u²/(2h) in between: an asymmetric Huber loss, within h/2 of ρ_τ, with curvature 1/h. With it, coordinate descent and proximal Newton fit penalized quantile regression through GLMLoss, without a linear program.

ProximalBase.ConvolutionCheck(τ, h; kernel = :epanechnikov)

The check loss of quantile regression convolved with a kernel of bandwidth h > 0, the smoothing of the conquer estimator (He, Pan, Tan and Zhou, 2023): ∫ ρ_τ(u - t) K(t/h)/h dt at the residual u = y - η. The kernels are :uniform, K = 1/2 on [-1, 1], and :epanechnikov, K(t) = 3(1 - t²)/4 on [-1, 1]. The loss equals ρ_τ(u) for |u| ≥ h and is a polynomial in between, with second derivative K(u/h)/h, so its curvature is 1/(2h) or 3/(4h). Unlike ProximalBase.SmoothCheck, it is centered: its derivative is τ - 1/2 at u = 0.

ProximalBase.Smoothed(φ, γ)

The Moreau envelope of a convex scalar function φ of the residual r = η - y, as a ProximalBase.ScalarLoss for GLMLoss:

ℓ(η, y) = min_t φ(t) + (t - r)² / (2γ),   γ > 0.

φ is a SeparableFunction with scalar parameters, of which the term of entry 1 is used. With p = prox_{γφ}(r), the value is φ(p) + (r - p)²/(2γ), the derivative (r - p)/γ and the second derivative (1 - ∂p/∂r)/γ, which is 1/γ or 0 on the pieces of a piecewise linear proximal operator; the curvature is at most 1/γ. So a nonsmooth loss such as the check loss of quantile regression or the hinge, smoothed, runs in coordinate descent with the curvature bound 1/γ, and in proximal Newton. Smoothed(NormL1(δ), 1) is ProximalBase.Huber(δ).

ProximalBase.deriv(ℓ, η, y)

The derivative ∂ℓ(η, y)/∂η of a ProximalBase.ScalarLoss.

ProximalBase.deriv2(ℓ, η, y)

The second derivative ∂²ℓ(η, y)/∂η² of a ProximalBase.ScalarLoss: the weight of the observation in a Newton step.

ProximalBase.curvature(ℓ)

An upper bound on the second derivative of the ProximalBase.ScalarLoss ℓ over all η and y, Inf when it is unbounded: with it, the gradient of Σᵢ wᵢ ℓ((Xβ)ᵢ, yᵢ) is Lipschitz with constant curvature(ℓ) ⋅ maxᵢ wᵢ ⋅ ‖X‖₂², and a coordinate step of size 1 / (curvature(ℓ) ⋅ Σᵢ wᵢ Xᵢⱼ²) decreases the loss.

ProximalBase.value_and_deriv(ℓ, η, y)

(ℓ(η, y), deriv(ℓ, η, y)), computed together.

IndIntersection(C₁, C₂, …; tol = 1e-10, maxiter = 10_000, warmstart = true)

The indicator of the intersection of the sets C₁, C₂, …, given by their indicator functions, such as IndPSD, IndBox or IndHyperplane. Its proximal operator, the projection onto the intersection, has no closed form even when each set has one. It is computed by Dykstra's algorithm (Boyle and Dykstra, 1986), which projects onto the sets in turn and corrects each projection by the increment of the previous round, so that it converges to the projection onto the intersection of convex sets, not merely to some point of it as plain alternating projections do.

The proximal operator is iterative and so inexact: it stops when a round over the sets moves the point by less than tol relative to its norm (at least ten machine epsilons in the precision of the output), or after maxiter rounds, with a warning. The point it returns lies exactly in the last set and within the tolerance of the others. f.iterations and f.converged report the last call.

With warmstart, each call starts from the corrections of the previous call, which Dykstra's algorithm allows, since it is block coordinate ascent on the dual of the projection (Tibshirani, 2017). Inside a solver the successive inputs are close, so the rounds per call drop as the solver converges. The corrections and the iterate live in buffers that f owns, sized by the first input, so later calls allocate nothing; do not share f between threads.

Each set is visited through its own prox!, its projection. With a nonconvex set, such as IndBallL0, the iteration is only a heuristic: it may stop at a point of the intersection that is not the projection, or cycle without converging.

For example, the nearest correlation matrix (Higham, 2002) is the projection onto the positive semidefinite matrices with unit diagonal:

L = fill(-Inf, p, p); L[diagind(L)] .= 1
U = fill(Inf, p, p);  U[diagind(U)] .= 1
C = prox(IndIntersection(IndPSD(), IndBox(L, U)), R)