Interface reference

ProximalBase

Building blocks for proximal, coordinate-descent and active-set methods: functions with gradients and proximal operators, and sparse iterates.

A function f is callable, f(x) being its value, and its traits say what it provides: issmooth functions have gradient! and value_and_gradient!, isproximable ones prox! and value_and_prox!, and isseparable ones the scalar proximal operators prox_entry, with which prox_coordinate! updates one coordinate of an iterate. A function may be both smooth and proximable.

Proximal operators follow the convention

prox_{γf}(x) = argmin_z  f(z) + ‖z - x‖² / (2γ),

and mutating functions take their output first, as mul! does: prox!(y, f, x, γ), gradient!(g, f, x), prox_coordinate!(x, f, k, γ).

AbstractFunction

Abstract supertype of the functions of ProximalBase: losses, penalties and indicator functions. A function f is callable: f(x) is its value, Inf outside its domain. Its traits say which oracles it provides:

A function can have several roles: SqrNormL2, QuadraticFunction, LeastSquaresLoss and LogDetLoss are both smooth and proximable. isconvex tells a solver whether convex theory applies.

To add a function, subtype AbstractFunction, make it callable, declare its traits for its type, as in ProximalBase.issmooth(::Type{<:MyLoss}) = true, and implement the oracles of those traits.

gradient!(g, f, x)

Store the gradient of f at x in g and return g.

gradient(f, x)

Allocating version of gradient!.

gradient_mapping!(G, f, g, x, γ)
gradient_mapping(f, g, x, γ)

The gradient mapping G = (x - prox_{γg}(x - γ∇f(x))) / γ of the problem minimize f(x) + g(x), for a smooth f and a proximable g. When f and g are convex it is zero exactly at the minimizers, for every γ > 0; with g = 0 it is ∇f(x). Its norm is a stopping rule that needs neither the optimal value nor the structure of g, and a proximal gradient step from x is x - γG. G must not be x; gradient_mapping! allocates nothing beyond what gradient! and prox! do.

gradient_mapping!(G, f, g, x, γ)
gradient_mapping(f, g, x, γ)

The gradient mapping G = (x - prox_{γg}(x - γ∇f(x))) / γ of the problem minimize f(x) + g(x), for a smooth f and a proximable g. When f and g are convex it is zero exactly at the minimizers, for every γ > 0; with g = 0 it is ∇f(x). Its norm is a stopping rule that needs neither the optimal value nor the structure of g, and a proximal gradient step from x is x - γG. G must not be x; gradient_mapping! allocates nothing beyond what gradient! and prox! do.

isconvex(f)

Whether f is convex. For a nonconvex f, prox! returns one of the minimizers that define the proximal operator, and accelerated methods such as FISTA lose their guarantees.

isproximable(f)

Whether f has an efficient proximal operator, computed by prox! and value_and_prox!.

isseparable(f)

Whether f(x) = Σᵢ fᵢ(xᵢ) is a sum over the entries of x. A separable f has the scalar proximal operators prox_entry, with which prox_coordinate! updates a coordinate of any iterate, and kkt_violation needs only the gradient.

issmooth(f)

Whether f is differentiable on the interior of its domain, with gradient! and value_and_gradient!. The gradient need not be globally Lipschitz: LogDetLoss is smooth on the positive definite matrices, and its value is Inf elsewhere.

kkt_violation(f, g, i)
kkt_violation(f, g, x, k)

How much a coordinate at zero violates the optimality conditions of h + f, where g = ∇h(x) is the gradient of the smooth part: a positive value means that moving the coordinate away from zero decreases the objective, and a value ≤ 0 that zero is optimal. For NormL1 it is |gᵢ| - λ0 wᵢ.

For a separable f the condition involves only gᵢ, and the first form suffices. The second form takes a coordinate k of the point x: for a GroupNormL2, the condition of an index involves the rest of its group, and on an AtomIterate k is a group, whose violation is ‖g[Gₖ]‖₂ - λ0 wₖ (for a separable f, the largest over the group).

lipschitz(f)

A Lipschitz constant L of the gradient of the smooth function f, so that a gradient or proximal gradient step with step 1/L decreases the objective. It is exact for the closed forms (SqrNormL2(λ) gives λ); for a loss of a matrix it is the largest eigenvalue involved, such as ‖X‖₂² for LeastSquaresLoss, computed by the Lanczos process with products by X and X' only, to a relative accuracy of 1e-10 and from below. It is Inf for a smooth function whose gradient is not globally Lipschitz, such as LogDetLoss and PoissonLoss, which call for a line search, and an error for a function that is not smooth. Sums add the constants, Scaled scales them, and SeparableSum takes their maximum.

prox(f, x, γ = 1)

Allocating version of prox!.

prox!(y, f, x, γ = 1)

Store the proximal operator

prox_{γf}(x) = argmin_z  f(z) + ‖z - x‖² / (2γ)

in y and return y. y may be x. For a nonconvex f the result is one of the minimizers.

prox_coordinate!(x, f, k, γ)

Apply the proximal operator of γf to coordinate k of the iterate x in place, as in one step of coordinate descent, and return the new value of the coordinate. The coordinates of a SparseIterate are its entries, those of a SymmetricSparseIterate the entries of its lower triangle, and those of an AtomIterate its groups, for which prox_coordinate! returns a view of the group (see ncoordinates).

A separable f updates any coordinate through prox_entry; a group penalty such as GroupNormL2 updates a group of an AtomIterate that shares its groups. A coordinate that stays zero does not become stored or active; one that becomes nonzero does.

prox_entry(f, v, i, γ)

For a separable f(x) = Σᵢ fᵢ(xᵢ) (see isseparable), the proximal operator of γfᵢ at the scalar v: the value that a coordinate-descent step gives entry i when the gradient step lands at v. i is a linear index of the array, so that weights and bounds may depend on the entry; for a SymmetricSparseIterate, entry (r, c) has index (c - 1)p + r. It touches no iterate: prox_coordinate! updates a coordinate of one in place.

restrict(f, S)

The function of x[S] that f becomes on the vectors x supported on the coordinates S, as a function of the same kind: an active-set solver minimizes the restriction over the active coordinates S. At x[S] the restriction has the value f(x) and the gradient ∇f(x)[S], except that SqrNormL2 with a center y drops the constant Σⱼ∉S yⱼ².

value_and_gradient!(g, f, x)

Store the gradient of f at x in g and return f(x). Usually cheaper than calling f(x) and gradient! separately.

value_and_gradient(f, x)

Allocating version of value_and_gradient!: returns (f(x), ∇f(x)).

value_and_prox!(y, f, x, γ = 1)

Store prox_{γf}(x) in y, as prox! does, and return f(y), the value at the result. Many functions get the value from the computation of the proximal operator: a nuclear norm from the singular values that it thresholds, an indicator function as 0. This is then cheaper than prox! followed by f(y), and solvers that track the objective should call it.

hashessian(f)

Whether the smooth function f provides second-order oracles: hessian_vec!, the product of its Hessian with a direction, and quadratic_model, its second-order model at a point as another function. Newton–CG, trust-region and proximal Newton methods need them; a solver checks this trait as it checks issmooth.

hessian_vec!(Hv, f, x, v)

Store the product ∇²f(x) v of the Hessian of f at x with the direction v in Hv and return Hv, for a function with hashessian. Use an output that does not overlap x or v, as for mul!. Hv is a dense array of the shape of x; x and v may be sparse iterates. The product costs about as much as a gradient: two products with X for LeastSquaresLoss and GLMLoss, w X⁻¹VX⁻¹ from one Cholesky factorization for LogDetLoss (NaNs when x is not positive definite), and it allocates nothing. On a SymmetricSparseIterate the product of KronQuadratic is the one over symmetric matrices, as its gradient is. Several products at the same x are cheaper through quadratic_model, which does the work that depends on x once.

hessian_vec(f, x, v)

Allocating version of hessian_vec!.

quadratic_model(f, x)
quadratic_model!(q, f, x)

The second-order model of the smooth function f at x, as a smooth function q of ProximalBase whose gradient at x is ∇f(x) and whose Hessian is ∇²f(x):

q(z) = f(x) + ⟨∇f(x), z - x⟩ + ⟨z - x, ∇²f(x)(z - x)⟩ / 2   (up to a constant)

A proximal Newton method minimizes q + g with any solver and penalty g, then searches along the step on f + g. quadratic_model! updates a model returned by quadratic_model(f, x₀) to the model at x, in place and without allocating, and returns it.

For a GLMLoss, ℓ''ᵢ is floored at √eps, as in IRLS, so that the working responses stay finite where the curvature vanishes (a Huber loss beyond δ); the gradient of q at x is still ∇f(x) exactly. For LogDetLoss, q is written in the new point z rather than in the step z - x: minimizing q(z) + λ‖z‖₁ over symmetric z is the Newton step of the QUIC algorithm for the graphical lasso (Hsieh, Sustik, Dhillon and Ravikumar, 2014), and coordinate descent on a SymmetricSparseIterate solves it. x must then be positive definite.

quadratic_model(f, x)
quadratic_model!(q, f, x)

The second-order model of the smooth function f at x, as a smooth function q of ProximalBase whose gradient at x is ∇f(x) and whose Hessian is ∇²f(x):

q(z) = f(x) + ⟨∇f(x), z - x⟩ + ⟨z - x, ∇²f(x)(z - x)⟩ / 2   (up to a constant)

A proximal Newton method minimizes q + g with any solver and penalty g, then searches along the step on f + g. quadratic_model! updates a model returned by quadratic_model(f, x₀) to the model at x, in place and without allocating, and returns it.

For a GLMLoss, ℓ''ᵢ is floored at √eps, as in IRLS, so that the working responses stay finite where the curvature vanishes (a Huber loss beyond δ); the gradient of q at x is still ∇f(x) exactly. For LogDetLoss, q is written in the new point z rather than in the step z - x: minimizing q(z) + λ‖z‖₁ over symmetric z is the Newton step of the QUIC algorithm for the graphical lasso (Hsieh, Sustik, Dhillon and Ravikumar, 2014), and coordinate descent on a SymmetricSparseIterate solves it. x must then be positive definite.

coordinate_lipschitz!(L, f)

Fill L, an array of the size of the argument of f, with the per-coordinate Lipschitz constants coordinate_lipschitz(f, k), and return L. It allocates nothing.

coordinate_lipschitz(f, k)

A Lipschitz constant Lₖ of the partial derivative ∂f/∂xₖ of the smooth function f along coordinate k, a linear index of its argument: the curvature of f along the coordinate, for which a coordinate step 1/Lₖ decreases f. It is the diagonal entry of the Hessian for the quadratic losses, a bound on it for GLMLoss, and Inf where none exists:

Scaled scales the constants, Tilt and Translate keep them, and Sum adds them. coordinate_lipschitz! computes all of them.

These constants are what coordinate descent needs. A gradient method that moves all coordinates at once can use them as a diagonal metric D = Diagonal(L), with the steps γ = t ./ L per entry in prox!, but the diagonal of a Hessian does not majorize it: the scalar t must come from a line search or be 1/λmax(D^(-1/2) ∇²f D^(-1/2)), which standardizing the columns of X amounts to.

reset_cache!(f)

Invalidate data-dependent caches of f and return f. Call this after changing the entries of data stored by reference in a LeastSquaresLoss (X or Y) or QuadraticFunction (A). Array sizes must stay unchanged; construct a new function for different dimensions. The next proximal call rebuilds the decomposition and may allocate; subsequent calls reuse it.

Call this on the outermost function you use: it descends through compositions, and an affine-constrained loss also resets the cache of its constrained solve. If several wrappers share the same loss, reset each wrapper after a data change. Values and gradients read the data directly and need no reset.

An IndAffine must instead be reconstructed when its data change: its factorization type depends on the rank of A. For an affine-constrained loss, reset_cache! supports changes to the loss data with the set unchanged. Other functions without data-dependent caches return themselves. Custom functions that cache derived data should implement this method.

ProximalBase.Testing

The checks that ProximalBase's own test suite runs on every function, for authors of new ones: a SeparableFunction of your own, or the functions that packages build on ProximalBase.

The checks run @test inside the caller's @testset, so a failure names the failing condition and its values. They are defined by a package extension on the Test standard library: load it with using Test. Random points come from the global random number generator, or from rng.

using Test, ProximalBase
const PT = ProximalBase.Testing
@testset "my function" begin
    f = MyPenalty(0.5)
    x = randn(10)
    PT.test_traits(f; proximable = true, separable = true, convex = true)
    PT.test_prox(f, x, 0.7)
    PT.test_separable(f, x, 0.7)
    PT.test_kkt(f, x, 0.7)
    @test PT.allocations(prox!, similar(x), f, x, 0.7) == 0
end

test_prox(f, x, γ; rtol = 1e-8, atol = 1e-10, local_points = true, symmetric = false, rng = nothing)

Test the proximal operator of f at x with step γ:

• prox, prox! and value_and_prox! agree, prox! returns its output and may overwrite its input, and the value is f at the result;
• the result p is optimal: no nearby point, nor the proximal point of another input, has a smaller proximal objective f(z) + ‖z - x‖²/(2γ), and for a convex f the subgradient inequality f(z) ≥ f(p) + ⟨(x - p)/γ, z - p⟩ holds at those points.

With symmetric = true the perturbations are symmetric matrices, for the functions of symmetric matrices; with local_points = false only proximal points of other inputs are compared (for functions whose domain is a thin set). Returns p.

test_gradient(f, x; h = 1e-6, rtol = 1e-5, directions = 3, rng = nothing)

Test the gradient of the smooth function f at x against central finite differences along random directions, and gradient!, gradient, value_and_gradient! and value_and_gradient against each other and against f(x). Returns the gradient.

test_separable(f, x, γ; groups = nothing, symmetric = false, rng = nothing)

Test a separable f at the vector or matrix x with step γ:

• prox_entry agrees with prox! entry by entry, and for a SeparableFunction the terms value_entry sum to f(x);
• prox_coordinate! on every coordinate of a SparseIterate (a copy of x with some entries set to zero) gives the proximal point;
• when iszeropreserving(f), zero is fixed by every scalar proximal operator and has value 0, and the value and the proximal operator on a SparseIterate, which touch only its stored entries, agree with the dense ones;
• with groups, prox_coordinate! on each group of an AtomIterate gives the proximal point;
• with symmetric = true and a square x, prox_coordinate! on each coordinate of a SymmetricSparseIterate gives the proximal point of the symmetric part of x, which needs parameters that are equal at (r, c) and (c, r).

test_conjugate(f, x, γ; conjugate = Conjugate(f), rtol = 1e-8)

Test the convex f against its conjugate at x with step γ:

• Moreau's decomposition x = prox_{γf}(x) + γ prox_{f*/γ}(x/γ), which tests a conjugate given with its own proximal operator (for Conjugate(f) it holds by construction);
• when the value of the conjugate is known, the Fenchel–Young equality f(p) + f*(y) = ⟨p, y⟩ at p = prox_{γf}(x) and its subgradient y = (x - p)/γ, and the inequality f(z) + f*(y) ≥ ⟨z, y⟩ at other points.

test_kkt(f, x, γ; trials = 5, rng = nothing)

Test kkt_violation against the proximal operator: from zero, a proximal gradient step with a random gradient g keeps a coordinate at zero exactly when its violation is ≤ 0. The coordinates are the entries of x for a separable f, and the groups for an AtomIterate x and a group penalty. For a nonconvex f, take γ below its convexity limit, where the scalar proximal operator at zero is the local one.

test_traits(f; smooth = false, proximable = false, separable = false, convex = false)

Test the traits of f, on the instance and on its type.

allocations(f, args...)

The bytes allocated by f(args...), measured inside a function after a warm-up call, as a solver makes the call: allocations(prox!, y, f, x, γ) == 0 is the contract of the hot calls. Counts differ between Julia versions; test on all that you support.

brute_force_prox(fᵢ, v, γ; lo, hi, n = 200_001)

The minimizer of fᵢ(t) + (t - v)²/(2γ) for a scalar function fᵢ, over a fine grid on [lo, hi] refined by golden-section search around the best grid point, and the minimum: (t, value). A reference for scalar proximal operators, nonconvex ones included.

SeparableFunction

The supertype of the separable functions f(x) = Σᵢ fᵢ(xᵢ), and the way to add one of your own. A subtype implements, for the entry with linear index i,

• prox_entry(f, v, i, γ): the proximal operator of γfᵢ at the scalar v;
• value_entry(f, v, i): the value fᵢ(v).

From these it gets the value f(x), prox! and value_and_prox! on any array, and prox_coordinate! on the three iterates. Declare iszeropreserving when fᵢ(0) = 0 and the proximal operator maps 0 to 0; then the value and the proximal operator on a sparse iterate touch only its stored entries. isseparable and isproximable are true; declare isconvex when it holds.

Example

struct WeightedSquare <: SeparableFunction            # Σᵢ wᵢ xᵢ² / 2
    w::Vector{Float64}
end
ProximalBase.value_entry(f::WeightedSquare, v, i) = f.w[i] * v^2 / 2
ProximalBase.prox_entry(f::WeightedSquare, v, i, γ) = v / (1 + γ * f.w[i])
ProximalBase.iszeropreserving(::WeightedSquare) = true
ProximalBase.isconvex(::Type{WeightedSquare}) = true

value_entry(f, v, i)

For a SeparableFunction f(x) = Σᵢ fᵢ(xᵢ), the value fᵢ(v) of the term of entry i at the scalar v. i is a linear index of the array, as in prox_entry.

iszeropreserving(f)

Whether the separable function f has fᵢ(0) = 0 and maps 0 to 0 by its proximal operator, prox_entry(f, 0, i, γ) = 0, for every entry i and step γ, as norms and most penalties do. Then the value and the proximal operator on a sparse iterate touch only its stored entries. It is false unless a function declares it, and it is a property of the function, not of its type: the bounds of an IndBox may exclude zero.

opnorm2(X)

Estimate the squared spectral norm ‖X‖₂² by Lanczos iteration, using only products with X and its adjoint, without forming a Gram matrix or densifying a sparse matrix. The result has floating element type matching X.

This is a numerical estimate, not a certified upper bound: Lanczos can miss an eigenspace orthogonal to its deterministic starting vector. Algorithms requiring a guaranteed bound should use an independently known bound (for example opnorm(X, 1) * opnorm(X, Inf)) or a suitable line search.