Penalties

Below, $S_t(v) = \operatorname{sign}(v)\max(|v| - t, 0)$ is soft thresholding and $\gamma$ is the step of the proximal operator.

The ℓ1 norm and its relatives

For NormL1, $g(x) = \lambda_0\sum_i w_i|x_i|$, the problem defining the proximal operator separates over the coordinates. Each scalar problem $\min_z \gamma\lambda_0 w_i|z| + \tfrac12(z - v_i)^2$ is solved by soft thresholding at $\gamma\lambda_0 w_i$. A weight $w_i = 0$ leaves the coordinate unpenalized, as for an intercept. The subdifferential of $|\cdot|$ at zero is $[-1, 1]$, so given the gradient $g_i$ of the smooth part, the coordinate $x_i = 0$ is optimal exactly when $|g_i| \le \lambda_0 w_i$. kkt_violation returns $|g_i| - \lambda_0 w_i$: positive when the coordinate should leave zero.

NormL1Nonneg adds the constraint $x \ge 0$, which turns soft thresholding one-sided. ElasticNet adds $\tfrac{\lambda_2}{2}\|x\|^2$. Completing the square, its proximal operator is soft thresholding followed by the shrinkage $1/(1 + \gamma\lambda_2)$.

Nonconvex penalties

ℓ0. The scalar problem $\min_z \gamma\lambda\,\mathbf{1}[z \ne 0] + \tfrac12(z - v)^2$ compares $z = 0$, of value $\tfrac12 v^2$, with $z = v$, of value $\gamma\lambda$. So the proximal operator keeps $v$ when $|v| > \sqrt{2\gamma\lambda}$ and sets it to zero otherwise.

MCP (Zhang, 2010). The penalty is $p(t) = \lambda|t| - t^2/(2a)$ for $|t| \le a\lambda$ and $a\lambda^2/2$ beyond, with $a > 1$. For $\gamma < a$ the scalar problem is convex, and its solution is firm thresholding:

\[\operatorname{prox}_{\gamma p}(v) = \begin{cases} 0 & |v| \le \gamma\lambda,\\ \operatorname{sign}(v)\,\dfrac{|v| - \gamma\lambda}{1 - \gamma/a} & \gamma\lambda < |v| \le a\lambda,\\ v & |v| > a\lambda . \end{cases}\]

Large values pass untouched, so MCP does not shrink large coefficients the way the ℓ1 norm does.

SCAD (Fan and Li, 2001). The penalty is $\lambda|t|$ for $|t| \le \lambda$, $(2a\lambda|t| - t^2 - \lambda^2)/(2(a - 1))$ for $\lambda < |t| \le a\lambda$, and $(a + 1)\lambda^2/2$ beyond, with $a > 2$. For $\gamma < a - 1$,

\[\operatorname{prox}_{\gamma p}(v) = \begin{cases} S_{\gamma\lambda}(v) & |v| \le (1 + \gamma)\lambda,\\ \dfrac{(a - 1)v - \operatorname{sign}(v)\,a\gamma\lambda}{a - 1 - \gamma} & (1 + \gamma)\lambda < |v| \le a\lambda,\\ v & |v| > a\lambda , \end{cases}\]

which for $\gamma = 1$ is Fan and Li's thresholding rule. Each piece comes from setting the derivative of $\gamma p(z) + \tfrac12(z - v)^2$ to zero on the corresponding interval.

For larger steps the scalar problems of MCP and SCAD are not convex. prox_entry still returns an exact minimizer: it compares the candidate minimizers of the pieces. So these penalties work in coordinate descent and proximal gradient methods with any step. At zero their subdifferential is that of $\lambda|\cdot|$, so kkt_violation is $|g_i| - \lambda$.

Capped ℓ1 (Zhang, 2010): $p(t) = \lambda\min(|t|, a)$. On $|z| \le a$ the scalar problem is the lasso's, so its minimizer is soft thresholding clipped at $a$; on $|z| \ge a$ the penalty is constant, so the minimizer is $\max(|v|, a)$. The proximal operator is the better of the two.

Log penalty: $p(t) = \lambda\log(1 + |t|/\varepsilon)$, whose linearization gives reweighted ℓ1 minimization (Candès, Wakin and Boyd, 2008). With $c = \gamma\lambda$ and $t = |v|$, the stationary points of $c\log(1 + z/\varepsilon) + \tfrac12(z - t)^2$ on $z > 0$ are the roots of $z^2 + (\varepsilon - t)z + c - \varepsilon t = 0$. The larger is a local minimum, compared with $z = 0$. Its slope at zero is $\lambda/\varepsilon$.

ℓq, $0 < q < 1$: $p(t) = \lambda|t|^q$. With $\mu = \gamma\lambda$ and $\beta = (2\mu(1 - q))^{1/(2 - q)}$, the proximal operator is zero for $|v| \le \tau = \beta(2 - q)/(2 - 2q)$, where zero and $\beta$ tie, and above it the root $z > \beta$ of $z + \mu q z^{q - 1} = |v|$, 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, which converges monotonically from $|v|$, otherwise. The slope at zero is infinite, so zero is a local minimizer of every problem and kkt_violation is $-\infty$.

Group MCP and group SCAD (Huang, Breheny and Ma, 2012) apply the MCP or SCAD with level $\lambda w_k$ to the Euclidean norm of each group. For any $p$, $p(\lVert z\rVert) + \lVert z - v\rVert^2/(2\gamma)$ is smallest, among the $z$ of a given norm $r$, at $z = r\,v/\lVert v\rVert$, so the proximal operator is the scalar one at $\lVert v_{G_k}\rVert$, along $v_{G_k}$: exact for every step.

Local linear approximation (reweighted ℓ1)

For a nonconvex penalty, the algorithm decides which local solution an estimator is, so recreating a paper means matching both its penalty and its algorithm. The papers that define the capped ℓ1 and log penalties, and the one-step and two-step SCAD and MCP estimators, compute them by local linear approximation (LLA), also called reweighted ℓ1.

The capped ℓ1, log, MCP and SCAD penalties are $f(x) = \sum_i p_i(\lvert x_i\rvert)$ with each $p_i$ concave and nondecreasing on $[0, \infty)$. So $f$ lies below its tangent at any point $\tilde x$:

\[p_i(\lvert t\rvert) \le p_i(\lvert\tilde x_i\rvert) + p_i'(\lvert\tilde x_i\rvert)\,\bigl(\lvert t\rvert - \lvert\tilde x_i\rvert\bigr).\]

A stage of LLA minimizes the loss plus the tangent, the weighted lasso with penalty $\sum_i w_i\lvert x_i\rvert$, $w_i = p_i'(\lvert\tilde x_i\rvert)$. Each stage is convex, so λ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$.

lla_weights!(w, f, x) fills w with these weights, and lla_weight(f, v, i) gives one. They include the level and the entry weights of f, so the penalty of the next stage is NormL1(1, w). For GroupMCP and GroupSCAD, w has one weight per group, $w_k = p_k'(\lVert x_{G_k}\rVert_2)$, and the next stage is GroupNormL2(1, f.groups, w).

At a kink the weight is the slope from the right, a supergradient of the concave $p_i$, so the majorization holds there too. At zero it is $p_i'(0+)$, which makes the first stage from $\tilde x = 0$ the lasso. For the capped ℓ1 penalty at $\lvert t\rvert = a$ it is $0$, the rule of T. Zhang's algorithm: $\lambda w_i$ for $\lvert t\rvert < a$, $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, and lla_weights! throws an error for it.

The papers, with their penalties as they define them:

Breheny and Huang (2011) compute MCP and SCAD paths instead by coordinate descent with the exact proximal operator of the penalty, which can give a different local solution.

With the loss $\tfrac12\lVert x - y\rVert^2$, a stage is weighted soft thresholding, and three stages of LLA for the MCP read:

using ProximalBase
y = [3.0, -0.4, 1.2, 0.1]
f = MCP(0.8, 2.0)
x = zeros(4)
for stage in 1:3
    w = lla_weights(f, x)          # the first stage, from zero, is the lasso
    x .= prox(NormL1(1, w), y)
end
x

4-element Vector{Float64}:
 3.0
 0.0
 0.6999999999999998
 0.0

The stages of a regression run in a solver; ProximalSolvers' reweight! alternates them with the weights.

• Candès, E. J., Wakin, M. B. and Boyd, S. P. (2008). Enhancing sparsity by reweighted ℓ1 minimization. Journal of Fourier Analysis and Applications 14, 877–905.
• Breheny, P. and Huang, J. (2011). Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Annals of Applied Statistics 5, 232–253.
• Fan, J., Xue, L. and Zou, H. (2014). Strong oracle optimality of folded concave penalized estimation. Annals of Statistics 42, 819–849.
• Shen, X., Pan, W. and Zhu, Y. (2012). Likelihood-based selection and sharp parameter estimation. Journal of the American Statistical Association 107, 223–232.
• Zhang, T. (2010). Analysis of multi-stage convex relaxation for sparse regularization. Journal of Machine Learning Research 11, 1081–1107.
• Zou, H. and Li, R. (2008). One-step sparse estimates in nonconcave penalized likelihood models. Annals of Statistics 36, 1509–1533.

Entropy, divergence and barrier functions

The negative entropy, the Kullback–Leibler divergence and the log barrier are separable, convex, and finite only on the nonnegative (or positive) orthant. For the first two, the optimality condition of the scalar problem with $\mu = \gamma\lambda$ is $z/\mu + \log z = c$ for a constant $c$. Writing $z = \mu s$ gives $s + \log s = c - \log\mu$, whose solution is the Wright omega function $\omega(u) = W_0(e^u)$, the Lambert $W$ function at $e^u$. ProximalBase computes $\omega$ without forming $e^u$, which overflows when the step is small and $u = v/\mu - 1 - \log\mu$ is large, so the proximal operators stay accurate there. The log barrier's scalar problem is the quadratic $z^2 - vz - \gamma\lambda = 0$, whose positive root is the one the proximal operator of LogDetLoss takes for each eigenvalue.

The reverse Huber penalty Berhu (Owen, 2007) is the ℓ1 norm up to $M$ and a ridge penalty beyond, with matching slopes at $\pm M$. It is convex and keeps the sparsity of the lasso near zero while shrinking large, correlated coefficients together. Its proximal operator soft-thresholds up to $|v| = M + \gamma\lambda$ and scales by $M/(M + \gamma\lambda)$ beyond.

Norms of vectors and groups

ℓ2. For $g = \lambda\|\cdot\|_2$, the proximal point lies on the segment from $0$ to $v$. Minimizing over its length gives block soft thresholding, $\max(1 - \gamma\lambda/\|v\|_2, 0)\,v$, which sets the whole vector to zero when $\|v\|_2 \le \gamma\lambda$.

ℓ∞. The conjugate of $\lambda\|\cdot\|_\infty$ is the indicator of the ℓ1 ball of radius $\lambda$. So by Moreau's decomposition its proximal operator is $v - P_{\|\cdot\|_1 \le \gamma\lambda}(v)$, through the ℓ1-ball projection of Sets.

Groups. GroupNormL2(λ0, groups, w) is $\lambda_0\sum_k w_k\|x_{G_k}\|_2$ over disjoint groups, the group lasso penalty. By default $w_k = \sqrt{|G_k|}$, and indices in no group are unpenalized. Its proximal operator applies block soft thresholding with threshold $\gamma\lambda_0 w_k$ to each group. A group at zero is optimal when $\|g_{G_k}\|_2 \le \lambda_0 w_k$. GroupNormLinf uses the ℓ∞ norm of each group, with the ℓ1-ball projection, and its optimality test is $\|g_{G_k}\|_1 \le \lambda_0 w_k$. On an AtomIterate that shares their groups, these penalties update one group at a time.

Sparse groups. SparseGroupNormL2(λ1, λ2, groups, w) adds an ℓ1 term, $\lambda_1\|x\|_1 + \lambda_2\sum_k w_k\|x_{G_k}\|_2$, the sparse group lasso of Simon, Friedman, Hastie and Tibshirani (2013): it selects groups, and entries within a selected group. Its proximal operator soft thresholds each entry at $\gamma\lambda_1$ and then applies block soft thresholding at $\gamma\lambda_2 w_k$ to each group,

\[z_{G_k} = \max\Bigl(1 - \frac{\gamma\lambda_2 w_k}{\|S_{\gamma\lambda_1}(v_{G_k})\|_2}, 0\Bigr) S_{\gamma\lambda_1}(v_{G_k}) .\]

Composing them is exact, because soft thresholding preserves the signs and only shrinks, so the subdifferential of the group norm at its output contains that at $z$. A group at zero is optimal when $\|S_{\lambda_1}(g_{G_k})\|_2 \le \lambda_2 w_k$, and an entry at zero inside a nonzero group when $|g_i| \le \lambda_1$. With one group, SparseGroupNormL2(λ1, λ2, [1:K], [1.0]), it is the sparse group lasso of a single vector: as the fiber function of a FiberSum it gives the penalty of the group graphical lasso (see below).

The rows of a $p \times k$ matrix are the groups [i:p:p*k for i in 1:p] of its linear indices. So GroupNormL2 with these groups is the row-sparse penalty of multitask regression; see the multitask example.

Fused pairs. NormL1Fused(λ1, λ2) acts on pairs of arrays $(a, b)$: $\sum_i \lambda_1(|a_i| + |b_i|) + \lambda_2|a_i - b_i|$, the penalty of the fused graphical lasso for two groups. Its proximal operator first solves the fusion part, moving $a_i$ and $b_i$ toward each other by $\gamma\lambda_2$ or to their mean. It then soft thresholds both by $\gamma\lambda_1$. Composing the two proximal operators this way is exact (Friedman, Hastie, Höfling and Tibshirani, 2007).

On a vector $v$ of length $K$, NormL1Fused(λ1, λ2) fuses all pairs of entries: $\lambda_1\lVert v\rVert_1 + \lambda_2\sum_{k<l}|v_k - v_l|$. For $K = 2$ it is the function of one pair. The proximal operator of the fusion part keeps the order of $v$, and on vectors sorted the same way the fusion part is linear: with $v_{(1)} \le \dots \le v_{(K)}$ it is $\lambda_2\sum_i (2i - K - 1) v_{(i)}$. So the proximal operator sorts $v$, subtracts $\gamma\lambda_2(2i - K - 1)$ from the $i$-th smallest entry, pools adjacent violators into a nondecreasing sequence, as for IndIsotonic, restores the order, and soft thresholds by $\gamma\lambda_1$. It is exact and takes $O(K \log K)$ operations.

Penalties across stacked matrices

The joint estimation of $K$ related matrices, such as the precision matrices of $K$ classes or time points, works on an $m \times n \times K$ array $x$. The fiber $x_{ij} = (x^1_{ij}, \dots, x^K_{ij})$ holds entry $(i, j)$ of every matrix. FiberSum(f; diagonal = f) is $\sum_{ij} f(x_{ij})$, and its proximal operator applies that of $f$ to each fiber. diagonal gives the fibers $x_{ii}$ of square slices their own function, and nothing leaves them unpenalized.

SliceSum(f₁, …, f_K) is $\sum_k f_k(x^k)$ for the slices $x^k = x[:, :, k]$, and its value, gradient and proximal operator act slice by slice. With $f_k$ the log-likelihood LogDetLoss(Sₖ, nₖ/n) of class $k$, the joint graphical lasso minimizes SliceSum(f₁, …, f_K) plus a FiberSum. ADMM splits the two: its $x$-update takes $K$ eigendecompositions, and its $z$-update is the proximal operator of the penalty, fiber by fiber. In the fused graphical lasso of Danaher, Wang and Witten (2014), with the diagonal fused but not shrunk, the penalty is FiberSum(NormL1Fused(λ1, λ2); diagonal = NormL1Fused(0, λ2)).

On a matrix, the fibers are the rows and the slices the columns. FiberSum(NormL2(λ)) is then the row-sparse penalty of multitask regression, without the index vectors of GroupNormL2.

Total variation, the fused lasso and SLOPE

These penalties depend on the order of the entries, and their proximal operators pool adjacent entries into blocks of equal values.

Total variation. The proximal operator of $\lambda\sum_i |x_{i+1} - x_i|$ is total-variation denoising. Its solution $z$ is piecewise constant. Its optimality conditions say that the cumulative sums of the residual $v - z$ stay within $\pm\gamma\lambda$ and end at zero, so the cumulative sums of $z$ are a taut string in a tube around those of $v$. Condat's algorithm (2013) scans $v$ once, keeping bounds on the value of the current segment, and goes back only when a segment ends. It is exact, works in place, and takes O(n) time on most inputs, O(n²) at worst.

The fused lasso. For $\lambda_1\|x\|_1 + \lambda_2\sum_i|x_{i+1} - x_i|$, soft thresholding the denoised vector by $\gamma\lambda_1$ gives the proximal operator (Friedman, Hastie, Höfling and Tibshirani, 2007): thresholding keeps the blocks of equal values, so their optimality conditions still hold.

SLOPE. The sorted ℓ1 norm $J_\lambda(x) = \sum_i \lambda_i|x|_{(i)}$, with $\lambda_1 \ge \dots \ge \lambda_n \ge 0$ and $|x|_{(1)} \ge |x|_{(2)} \ge \cdots$, is the ℓ1 norm when the $\lambda_i$ are equal and the ℓ∞ norm when $\lambda = (\lambda_1, 0, \dots, 0)$. With the Benjamini–Hochberg sequence $\lambda_i = \sigma\Phi^{-1}(1 - iq/(2p))$, SLOPE controls the false discovery rate of the selected variables for orthogonal designs (Bogdan, van den Berg, Sabatti, Su and Candès, 2015). The proximal operator keeps the signs and the order of the magnitudes of $v$. With $|v|$ sorted in decreasing order, it is the nonincreasing, nonnegative sequence nearest to $|v|_{(i)} - \gamma\lambda_i$: isotonic regression by pool adjacent violators, clipped at zero, in O(n log n) time for the sort. Entries pooled into a block come out with equal magnitudes, which is how SLOPE groups correlated predictors.

using ProximalBase
v = [3.0, -1.0, 0.5, 2.0]
prox(SLOPE([1.0, 0.8, 0.5, 0.2]), v, 1.0)    # magnitudes 3, 2, 1, 0.5 minus λ

4-element Vector{Float64}:
  2.0
 -0.5
  0.3
  1.2

Losses of residuals and margins

Quantile regression, support vector machines and support vector regression have losses that are nonsmooth but separable, with closed-form proximal operators. So splitting methods such as ADMM apply them to $z = X\beta$ exactly.

The check loss. CheckLoss(τ) is $\sum_i \rho_\tau(x_i)$ with $\rho_\tau(u) = \max(\tau u, (\tau - 1)u)$, for $\tau \in (0, 1)$. Its proximal operator is

\[\operatorname{prox}_{\gamma\rho_\tau}(v) = \begin{cases} v - \gamma\tau & v > \gamma\tau,\\ v + \gamma(1 - \tau) & v < -\gamma(1 - \tau),\\ 0 & \text{otherwise}. \end{cases}\]

Since $\rho_\tau(-u) = \rho_{1-\tau}(u)$, the loss $\sum_i \rho_\tau(y_i - z_i)$ of predictions $z$ is Translate(CheckLoss(1 - τ), y). The conjugate of $\rho_\tau$ is the indicator of $[\tau - 1, \tau]$, the dual feasible set of quantile regression, and Conjugate gives its value.

The hinge loss. HingeLoss(y) is $\sum_i \max(0, 1 - y_i x_i)$ for labels $y_i = \pm 1$. At the margin $u = y_i v$, the proximal operator of $\gamma\max(0, 1 - u)$ is $u + \gamma$ for $u < 1 - \gamma$, $1$ for $1 - \gamma \le u \le 1$, and $u$ above. Without labels, $x$ holds the margins.

The ε-insensitive loss. EpsInsensitiveLoss(ε) is $\sum_i \max(0, |x_i| - \varepsilon)$. Its proximal operator keeps $v$ for $|v| \le \varepsilon$, puts it at $\pm\varepsilon$ for $\varepsilon < |v| \le \varepsilon + \gamma$, and moves it toward zero by $\gamma$ beyond. With $\varepsilon = 0$ it is soft thresholding.

All three take weights $w_i$, which multiply the step. For coordinate descent and Newton steps, the scalar losses ProximalBase.SquaredHinge, ProximalBase.SmoothCheck and ProximalBase.ConvolutionCheck of GLMLoss are smooth versions with bounded curvature; see Smooth functions.

Spectral penalties

For a function of the singular values, $F(X) = f(\sigma(X))$, with $f$ absolutely symmetric, the proximal operator acts on the singular values: $\operatorname{prox}_{\gamma F}(X) = U\operatorname{diag}(\operatorname{prox}_{\gamma f}(\sigma))V^\top$ for $X = U\operatorname{diag}(\sigma)V^\top$ (Lewis, 1995). For NuclearNorm, $f = \lambda\|\cdot\|_1$, which gives singular value thresholding. TracePSD, $\lambda\operatorname{tr}(X)$ on the positive semidefinite matrices, is the analogous function of the eigenvalues: it subtracts $\gamma\lambda$ from each eigenvalue and clips at zero. It is the penalty on the low-rank part of latent-variable graphical models. value_and_prox! gets the value from the thresholded spectrum. The decompositions run in LAPACK workspaces that the functions own. Spectral(f) builds the function of the eigenvalues of a symmetric matrix from any function $f$ of a vector that does not depend on the order of its entries.

SingularValue(f) is the general function of the singular values, for any absolutely symmetric vector function $f$ with a proximal operator: one thin SVD, then the proximal operator of $f$ on the singular values. The nuclear norm is SingularValue(NormL1(λ)). The others:

For a nonconvex separable $f$ whose scalar proximal operator is nondecreasing on $[0, \infty)$, applying it to the singular values still gives a global minimizer (Lu, Zhu, Xu, Yan and Lin, 2015). So MCP and SCAD thresholding of the singular values reduce the bias of singular value thresholding, as they do for the lasso. The conjugate of $f \circ \sigma$ is $f^* \circ \sigma$, so the spectral norm and the nuclear-norm ball are a conjugate pair under Conjugate. The Ky Fan $k$-norm, the sum of the $k$ largest singular values, is the support function of $\{w : 0 \le w_i \le \lambda,\ \sum_i w_i \le \lambda k\}$ at $\sigma$. It needs the halfspace: with the hyperplane, the conjugate of the capped simplex is the sum of the $k$ largest entries, which is not absolutely symmetric.

Partial decompositions

NuclearNorm, TracePSD, IndPSD and IndBallRank keep only the singular triplets or eigenpairs above a threshold, or the $r$ largest, but by default compute the whole decomposition, in $O(p^3)$ for a $p \times p$ matrix. With the keyword partial = true they compute only the pairs they keep, by randomized subspace iteration with a Rayleigh–Ritz step (Halko, Martinsson and Tropp, 2011): for a result of rank $k$, each iteration costs $O(p^2\ell)$ on a block of $\ell \approx k + 10$ vectors.

• Warm starts. A call starts from the block of the previous call, after one more power step, with two fresh random vectors. A solver calls the proximal operator on nearby matrices, so a call takes one to a few iterations, and a direction that enters the leading subspace is amplified from call to call.
• Accuracy. A result is accepted when the residuals $\|Xv_j - \sigma_j u_j\|$ (or $\|Au_j - \theta_j u_j\|$) of the pairs kept are below $\epsilon^{3/4}$ times the largest value, and the next three pairs are certified to lie below the threshold: their value plus their residual is below it. The result then agrees with the full decomposition to about $10^{-12}$ in Float64.
• Unknown ranks. For the thresholds of NuclearNorm, TracePSD and IndPSD(), the block starts from the rank of the previous call and doubles while too few of its values fall below the threshold. For symmetric matrices, it iterates on $A + cI$ with a shift $c$ when negative eigenvalues are large enough to crowd out the ones kept.
• Fallback. The function takes the full decomposition when the block would exceed a quarter of the dimension, or when the iteration does not converge, and that decomposition starts the next call. partial = true is therefore safe to use on any input, and pays when the rank of the result is small.
• What stays full. The value $f(X)$ alone needs every singular value or eigenvalue, and still takes the full decomposition. As for the full decompositions, the workspaces belong to the function object, so do not share one between threads.

On a low-rank matrix plus noise, one call on a slightly changed input, as in a solver, takes (Julia 1.13, one thread of an AMD EPYC 7513; benchmark/partial.jl):

In parentheses, the factor by which this is faster than the full decomposition; at $p = 5000$ the full SVD and eigendecomposition were timed once, for NuclearNorm and TracePSD, and IndBallRank and IndPSD() take the same.

Over ranks $k = 2, 10, 50$ and $p$ from 200 to 5000, a warm call is faster than the full decomposition wherever the partial path applies, and so is a call from a random block, except for TracePSD and IndPSD() at $p = 500$ with $k = 50$, about 7% slower. Told the exact rank, ARPACK and KrylovKit compute the leading pairs of one matrix up to twice as fast for $p \le 2000$, and up to three times as fast for IndPSD(rank = k), but their SVDs are slower than the partial one at $p = 5000$; they allocate at every call, need the rank in advance, and do not start from the previous call. The Lanczos SVD of TSVD is often faster but accurate to about $10^{-7}$, and randomized SVDs without power iterations to a few percent (see benchmark/README.md).