Calculus reference

Conjugate(f)

The convex conjugate f*(y) = supₓ ⟨x, y⟩ - f(x) of a convex f. Its proximal operator follows from that of f by Moreau's identity, prox_{γf*}(x) = x - γ prox_{f/γ}(x/γ), so every proximable function gives one; for a separable f so does prox_entry. Its value is known for the norms (the indicators of the dual balls), the balls (their support functions), the box, the nonnegative orthant, SqrNormL2, LinearFunction, ZeroFunction, and the losses CheckLoss, HingeLoss and EpsInsensitiveLoss, and for translations and strictly positive rescalings of these functions. The conjugate value of a zero rescaling is not implemented, since it depends on the domain of the underlying function.

Scaled(f, α)
α * f

The function α f(x) for α ≥ 0, with the roles of f: its gradient is α ∇f, and its proximal operator prox_{γαf}, so that 0.5 * f halves a loss and (1/n) * f averages one.

SeparableSum(f1, f2, ...)

The function Σₖ fₖ(xₖ) of a tuple x = (x1, x2, ...) of blocks, with the roles that all the fₖ share: the proximal operator, gradient and value act block by block on tuples. Sparse plus low-rank models, for example, take SeparableSum(NormL1(λ), TracePSD(μ)) on the pair (S, L).

Sum(f1, f2, ...)

The function Σₖ fₖ(x), smooth when all the fₖ are, with the sum of their gradients: a loss plus a smooth penalty, for example. A sum has no proximal operator in general.

Tilt(f, c, d = 0)

The function f(x) + ⟨c, x⟩ + d, with the roles of f. Its proximal operator is that of f at x - γc, and its gradient ∇f + c. The linear terms of augmented Lagrangians and of duality take this form.

Translate(f, b)

The function f(x - b), with the roles of f: its gradient is ∇f(x - b) and its proximal operator b + prox_{γf}(x - b). For example, Translate(NormL1(λ), β₀) is the penalty λ ‖β - β₀‖₁ that shrinks towards β₀.

Constrained(f, C)

The function f(x) on the set C, and Inf outside it, for a convex separable f and a set C that is

• a box: IndBox, IndNonnegative or IndBallLinf;
• a hyperplane or a halfspace: IndHyperplane or IndHalfspace;
• a simplex, IndSimplex: the box x ≥ 0 with the hyperplane Σᵢ xᵢ = r, or with the halfspace Σᵢ xᵢ ≤ r;
• for a norm, an ℓ2 ball or sphere, IndBallL2 or IndSphereL2, described below;
• for a quadratic, an affine set IndAffine, described below.

On a box, the proximal operator clamps that of f entry by entry, so Constrained(f, box) is separable, with prox_entry for coordinate descent. On a hyperplane ⟨a, x⟩ = b, the proximal operator is prox_{γf}(x - γνa) for the multiplier ν at which it satisfies the constraint, found by a search over the scalar ν; f may itself be constrained to a box. Each step of the search is a pass over x. It starts from the multiplier of the previous call, which the function keeps, with a first step from the exact slope of the constraint where the scalar proximal operators of f are piecewise linear, as for every function below. The calls of a solver on slowly changing points then take three to five passes, against 7 to 16 from zero. From zero, on 4096 entries or more, where the scalar proximal operators of f stop at a bound, as on a box, the nonnegative orthant or a simplex, the search fixes the entries that have reached one once few others stay, and its later passes skip them: when the multiplier lies far out in the tail of the entries, as for a simplex of radius 1, such a call takes about six passes and is 1.6 to 2.9 times faster. Do not share one function object between threads.

On a box, a hyperplane, a halfspace or a simplex, γ may also be an array of steps, one per entry: the proximal operator in the metric of Diagonal(1 ./ γ), by the same search with prox_{γᵢfᵢ}(xᵢ - γᵢνaᵢ) in entry i. For example:

• Constrained(NormL1(λ), IndHyperplane(1, 0)) is the lasso penalty on coefficients that sum to zero, as in regression on compositional data;
• Constrained(NormL1(λ), IndBox(0, u)) is the nonnegative lasso penalty with upper bounds;
• Constrained(IndBox(0, 1), IndHyperplane(1, k)) is the indicator of the capped simplex {x : 0 ≤ x ≤ 1, Σᵢ xᵢ = k};
• Constrained(SqrNormL2(λ), IndSimplex(1)) is the ridge penalty on the probability simplex.

The conjugate of the indicator of a box with scalar bounds intersected with a hyperplane or a halfspace with a scalar a, such as the capped simplex, has a value: the support function of the set. For the capped simplex it is the sum of the k largest entries, and SingularValue and Spectral turn it into the Ky Fan k-norm and the sum of the k largest eigenvalues.

f must be convex: for a nonconvex penalty, clamping its proximal operator does not give the constrained one. The exceptions are MCP, SCAD and NormL0 on a box, whose scalar problems are solved exactly by comparing the minimizers of their convex pieces with the ends of the interval; the result is separable but not convex. When the set misses the domain of f, the proximal operator throws an ArgumentError.

Constrained(f, IndBallL2(r))
Constrained(f, IndSphereL2(r))

The function f(x) on the Euclidean ball ‖x‖₂ ≤ r or on the sphere ‖x‖₂ = r (the Frobenius ball and sphere for matrices), and Inf outside it: the problems of sparse principal component and canonical correlation analysis, whose alternating steps put a shrunken vector on the ball or the sphere.

On the ball, f is positively homogeneous: a norm (NormL1, NormL1Nonneg, NormL2, NormLinf, GroupNormL2, GroupNormLinf, NuclearNorm, TracePSD, or a Scaled one) or the indicator of a cone (IndNonnegative, IndPSD, and the nonconvex IndBallL0 and IndBallRank). The proximal operator is the proximal point of f scaled onto the ball, and the value scales with it: f(tz) = t f(z) for t > 0.

On the sphere, f is NormL1, GroupNormL2 or IndBallL0, and the proximal operator is the global minimizer of the nonconvex problem: the proximal point of f scaled onto the sphere. When that point is zero, so that every entry of x is within its threshold, the minimizer puts all the mass on one entry (one group) where |xᵢ| exceeds the threshold the least, along the sign of xᵢ (the direction of x on the group).

u = prox(Constrained(NormL1(λ), IndBallL2(1)), X * v)        # the step of the penalized matrix decomposition
u = prox(Constrained(IndBallL0(k), IndSphereL2(1)), S * u)    # the truncated power method

Constrained(f, IndAffine(A, b))

The quadratic f, a LeastSquaresLoss with a vector response or a QuadraticFunction, on the affine set {x : Ax = b}, and Inf outside it. Its proximal operator is the equality-constrained ridge step

x = M(c - Aᵀλ),   (AMAᵀ)λ = AMc - b,

with M = (H + I/γ)⁻¹ for the Hessian H of f and c = v/γ - ∇f(0): for the least-squares loss ‖Y - Xβ‖²/2, H = X'X and c = X'Y + v/γ.

M is applied through the decomposition that f keeps for its own proximal operator, the thin SVD of X or the eigendecomposition of the quadratic's matrix, computed on first use. The m × m matrix AMAᵀ depends on γ, and the function keeps its Cholesky factor for the last γ. A call with the same γ costs the unconstrained step plus two products with an m × n matrix, and a new γ adds the factorization of AMAᵀ; neither allocates. A rank-deficient A takes the equivalent constraint of full row rank from the thin SVD that IndAffine caches. Do not share one function object between threads.

With ADMM, f = Constrained(LeastSquaresLoss(Y, X), IndAffine(A, b)) and g = NormL1(λ) give the lasso under equality constraints, such as log-contrast regression with several subcompositions, Aβ = 0 with a row of A for each.

Spectral(f)

The function F(X) = f(λ(X)) of a symmetric matrix X through its eigenvalues λ(X), for a proximable function f of a vector that does not depend on the order of the entries. Its proximal operator acts on the eigenvalues: for X = V diag(λ) V', prox_{γF}(X) = V diag(prox_{γf}(λ)) V' (Lewis, 1996), and F is convex when f is. For example:

• Spectral(Constrained(IndBox(0, 1), IndHyperplane(1, k))) is the indicator of the Fantope {X : 0 ⪯ X ⪯ I, tr(X) = k}, the convex hull of the projections of rank k, as in sparse principal component analysis (Vu, Cho, Lei and Rohe, 2013);
• Spectral(Conjugate(Constrained(IndBox(0, 1), IndHyperplane(1, k)))) is the sum of the k largest eigenvalues, the support function of the Fantope;
• Spectral(IndBox(l, u)) is the indicator of the matrices with eigenvalues in [l, u];
• Spectral(IndSimplex(1)) is the indicator of the density matrices, {X ⪰ 0 : tr(X) = 1};
• Spectral(NormL1(λ)) is λ times the trace norm, the nuclear norm of a symmetric matrix.

Spectral(IndNonnegative()) is IndPSD(), and Spectral(NormL1Nonneg(λ)) is TracePSD(λ). The proximal operator reads the upper triangle of its input, which should be symmetric, and the value is Inf at a matrix that is not. A function with weights, bounds or normals that differ between entries depends on their order and is rejected.

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.

MoreauEnvelope(f, γ)

The Moreau envelope e(x) = min_z f(z) + ‖z - x‖² / (2γ) of a convex, proximable f, for γ > 0: a smooth function below f that approaches it as γ → 0. With p = prox_{γf}(x), its value is f(p) + ‖x - p‖² / (2γ) and its gradient (x - p) / γ, Lipschitz with constant 1/γ, so a value and a gradient cost one proximal operator of f (and no second evaluation of f, which comes from value_and_prox!).

It is also proximable, by prox_{te}(x) = x + t/(γ + t) ⋅ (prox_{(γ + t)f}(x) - x), and separable when f is, so coordinate descent applies to it. The envelope of δ|⋅| is the Huber function, that of the indicator function of a set C the squared distance to C divided by 2γ. The value, the gradient and the proximal operator use a buffer that e owns; do not share e between threads, and pass a gradient output that is not x.

Precompose(f, A, μ; check = true)
Precompose(f, a)

The function h(x) = f(Ax) of a vector x, for a linear map A with AA' = μI and μ > 0. Its proximal operator follows from that of f:

prox_{γh}(x) = x + A'(prox_{μγf}(Ax) - Ax) / μ,

and since A ⋅ prox_{γh}(x) = prox_{μγf}(Ax), the value at the result comes from that of f. A is anything with mul!(y, A, x) and the five-argument mul!(y, A', z, α, β): a matrix, a sparse matrix or an operator object. For a matrix, the constructor checks AA' = μI unless check = false. Rows of the identity (a selection of entries), permutations, orthogonal transforms and tight frames qualify. Smooth when f is, with gradient A'∇f(Ax) and Lipschitz constant μ L_f; proximable and convex when f is. When f has second-order oracles, so does h, with Hessian product A'∇²f(Ax)Av.

With a scalar a ≠ 0, Precompose(f, a) is f(ax) on arrays of any shape (μ = a²), separable when f is, with coordinate curvatures scaled by a²; with Translate it gives TFOCS's f(ax + b). The value, the gradient and the proximal operator use buffers that h owns; do not share h between threads.

Regularize(f, ρ, c = nothing)

The function f(x) + (ρ/2)‖x - c‖² for ρ ≥ 0, with c = 0 when omitted, with the roles of f. Its proximal operator is that of f with a shorter step, at a shrunk point,

prox_{γ(f + ρ‖⋅ - c‖²/2)}(v) = prox_{γ/(1 + γρ) f}((v + γρc) / (1 + γρ)),

and the value at the result comes from that of f. Smooth when f is, with gradient ∇f + ρ(x - c) and Lipschitz constant L_f + ρ; separable when f is, so coordinate descent applies. Hessian products and quadratic models are available when f provides them, and coordinate curvatures increase by ρ. On an AtomIterate a regularized group norm updates a group. Regularize(NormL1(λ), ρ) is ElasticNet(λ, ρ), Regularize(GroupNormL2(λ, groups), ρ) the group elastic net, Regularize(NuclearNorm(λ), ρ) the matrix elastic net, and with a center c it is the strongly convex subproblem of proximal-point methods.