Fast and Flexible ADMM Algorithms for Trend Filtering

Aaditya Ramdas; Ryan J. Tibshirani

Source-linked AI summary

Fast and Flexible ADMM Algorithms for Trend Filtering

Aaditya Ramdas, Ryan J. Tibshirani

arXiv:1406.2082v4Jun 9, 2014stat.MLcs.LGcs.NAmath.NAmath.OCstat.AP

TL;DR

Trend filtering has minimax-optimal error rates for functions with bounded-variation derivatives, but fitting it requires scalable and numerically robust algorithms. This paper introduces a specialized ADMM routine and reports more reliable convergence than PDIP across a broad range of settings, alongside visually acceptable estimates after few iterations.

Problem
Trend filtering is minimax optimal for functions with bounded-variation derivatives, but the paper identifies a need for scalable and numerically stable algorithms to fit its estimates.
Source: §1 · p. 2, §1.2 · p. 3
Method
The paper develops a specialized ADMM algorithm that uses a recursive parametrization to leverage fast exact solvers for the k=0 fused lasso case.
Source: §6 · p. 18
Results
Specialized ADMM converges more reliably than PDIP across broad ranges of n and λ, and produces visually acceptable estimates after few iterations even where criterion convergence is slower.
Source: §1.1 · pp. 2–3
Takeaways & Limitations
The specialized ADMM implementation is simple and can be extended to variants including sparse and isotonic trend filtering.
Source: §6 · p. 18
In particular settings—small problems and small λ for moderate or large problems—specialized ADMM does not match PDIP's convergence rate; more uneven input spacings also slow convergence.
Source: §1.1 · p. 3, §4 · p. 15

Abstract

from arXiv · show

This paper presents a fast and robust algorithm for trend filtering, a recently developed nonparametric regression tool. It has been shown that, for estimating functions whose derivatives are of bounded variation, trend filtering achieves the minimax optimal error rate, while other popular methods like smoothing splines and kernels do not. Standing in the way of a more widespread practical adoption, however, is a lack of scalable and numerically stable algorithms for fitting trend filtering estimates. This paper presents a highly efficient, specialized ADMM routine for trend filtering. Our algorithm is competitive with the specialized interior point methods that are currently in use, and yet is far more numerically robust. Furthermore, the proposed ADMM implementation is very simple, and importantly, it is flexible enough to extend to many interesting related problems, such as sparse trend filtering and isotonic trend filtering. Software for our method is freely available, in both the C and R languages.

1 Introduction

Trend filtering estimates nonparametric functions with adaptively placed knots and, for functions with bounded-variation derivatives, is described as minimax optimal while computationally comparable to smoothing splines. The paper develops a specialized ADMM algorithm to make fitting trend filtering efficient and robust, and compares it with existing optimization methods.

1 Introduction

Order-k trend filtering estimates piecewise-polynomial functions with adaptively selected knots; at k=0, the estimate is piecewise constant.At k=0, trend filtering is the one-dimensional fused lasso, also called total variation denoising.Source: §1 · pp. 1–2
For estimating functions whose derivatives are of bounded variation, trend filtering is minimax optimal and computationally comparable to smoothing splines, which are not minimax optimal.The comparison also describes locally adaptive regression splines as minimax optimal but relatively computationally inefficient.Source: §1 · p. 2

1.1 Overview of contributions

The paper addresses solving trend filtering for a given λ, not choosing λ, though fast solves can speed cross-validation.Source: §1.1 · p. 2
Specialized ADMM converges more reliably than PDIP across broad ranges of n and λ, while producing visually acceptable estimates after few iterations even when its criterion converges more slowly.PDIP is faster in particular settings, including small problems and small λ for moderate or large problems.Source: §1.1 · pp. 2–3

1.2 A motivating example

Specialized ADMM greatly improves convergence over standard ADMM and, in a motivating simulation, reached in 20 iterations what first-order methods and coordinate descent did not reach in thousands.The paper attributes iterative difficulty to ill-conditioned difference operators, whose conditioning worsens as order k increases.Source: §1.1 · p. 3, §1.2 · pp. 3–4

1.1 Overview of contributions

Each specialized ADMM iteration has linear complexity in n and runs empirically about 10 times faster than a PDIP iteration.The paper reports linear complexity for PDIP iterations as well.Source: §1.1 · p. 3
The simple ADMM implementation can extend to sparse, mixed, and isotonic trend filtering, and univariate methods can serve as building blocks for broader models such as generalized additive models.Source: §1.1 · p. 3

2 A specialized ADMM algorithm

The specialized ADMM changes trend filtering’s splitting while retaining linear-time iterations, using direct fused-lasso solvers to improve convergence over standard ADMM. The section also examines why the updates converge differently, warm starts, and the choice of the ADMM parameter ρ.

2 A specialized ADMM algorithm

Both ADMM implementations have O(n) work per iteration, but the specialized method matches the exact fit visually after 20 iterations in the n = 1000 example.The standard ADMM approach in the figure uses 5,000 iterations; the specialized PDIP and ADMM methods use 20.Source: §2 · p. 4, Figure 1 · p. 5
The specialized updates wrap exact linear-time fused-lasso solvers in ADMM, lending their efficiency to higher-order trend filtering.The authors use a dynamic-programming routine for the fused-lasso subproblem.Source: §2 · p. 6, §2.2 · p. 6

2.1 Superiority of specialized over standard ADMM

Across simulation settings, specialized ADMM converges to the optimum faster than standard ADMM, with the gap generally widening as k increases.The reported qualitative pattern holds across the simulation suite and regularization regimes.Source: §2.1 · p. 6

2.2 Some intuition for specialized versus standard ADMM

In basis form, the specialized parametrization has less-correlated regression-matrix blocks than the standard one, which the authors expect to let each update make more progress.The authors describe the specialized updates as decorrelated versions of the original ADMM updates.Source: §2.2 · pp. 7–8

2.3 Superiority of warm over cold starts

Warm starts improve both algorithms’ absolute performance without changing their relative comparison; specialized ADMM gains most over cold starts in the middle λ range.The experiments solve the largest λ first and use each solution to initialize the next.Source: §2.3 · p. 8, Figure 3 · p. 8

2.4 Choice of the augmented Lagrangian parameter ρ

Setting ρ = λ was numerically stable across setups, while adaptively varying ρ was less stable overall and gave neither algorithm consistent benefits.The reported experiments use ρ = λ by default unless stated otherwise.Source: §2.4 · p. 9

3 Comparison of specialized ADMM and PDIP

The specialized ADMM and PDIP routines both take linear time per iteration, but ADMM has a smaller per-iteration cost and is more robust for higher-order trends. The comparison attributes PDIP’s instability to poorly conditioned linear systems, while ADMM’s update includes an eigenvalue buffer.

3 Comparison of specialized ADMM and PDIP

The proposed comparisons use specialized C ADMM and a generalized C++/LAPACK PDIP implementation, with varied PDIP settings failing to yield consistent improvement where it struggled.The PDIP method is applied to the dual problem and solves banded linear systems.Source: §3 · p. 9

3.1 Comparison of cost per iteration

Both routines take O(n) time per iteration, but an ADMM iteration is about 10 times faster than a PDIP iteration.Figure 4 compares time per iteration across problem sizes, trend orders, functions, and λ values.Source: §3.1 · p. 9, Figure 4 · p. 10

3.2 Comparison for k = 1 (piecewise linear fitting)

For k = 1, specialized ADMM and PDIP perform similarly, although PDIP has a relative advantage at small λ values that correspond to overfitting.ADMM’s convergence remains strong in absolute terms.Source: §3.2 · p. 10, Figure 5 · p. 11

3.3 Comparison for k = 2 (piecewise quadratic fitting)

For k = 2, specialized ADMM converges strongly across λ values and problem sizes into the hundreds of thousands, while PDIP struggles at moderate to large λ; k = 3 amplifies PDIP’s problems.The k = 3 issues appear at smaller problem sizes, while ADMM continues to converge steadily.Source: §3.3 · p. 10, Figure 6 · p. 12

3.4 Some intuition on specialized ADMM versus PDIP

ADMM’s system adds an identity-matrix buffer, keeping eigenvalues at least 1, whereas PDIP’s diagonal matrix can lose that buffer and leave its system poorly conditioned.The paper links this PDIP instability particularly to large λ, where many diagonal entries can approach zero.Source: §3.4 · p. 13

4 Arbitrary input points

For arbitrary input locations, the specialized ADMM algorithm replaces the evenly spaced difference operator with a spacing-adjusted one. More irregular spacing worsens conditioning and slows convergence, while an input-adjusted choice of ρ helps the algorithm progress.

4 Arbitrary input points

The arbitrary-input extension replaces the evenly spaced difference operator with a recursively defined, spacing-adjusted operator that retains the banded structure.The adjusted operator is substituted into the specialized ADMM updates.Source: §4 · p. 13
As input spacings become more uneven, conditioning worsens and ADMM converges more slowly, though it remains fairly robust on considerably irregular designs.Figure 7 illustrates convergence and fitted estimates for evenly spaced, random, and mixture-of-Gaussians inputs.Source: Figure 7 · p. 14, §4 · p. 15

4.1 Choice of the augmented Lagrangian parameter ρ

For arbitrary inputs, the authors advocate an input-adjusted ρ that essentially reduces to ρ = λ for unit-spaced inputs.The choice aims to match the first ADMM update of a corresponding unit-spaced routine and makes an important practical difference to algorithm progress.Source: §4.1 · p. 15

5 ADMM algorithms for trend filtering extensions

The paper presents ADMM as a flexible framework for modifying the basic trend-filtering model, while emphasizing that its extension examples are illustrative rather than exhaustive or statistical evaluations.

5 ADMM algorithms for trend filtering extensions

The ADMM framework readily accommodates modifications to trend filtering, but the examples are illustrative rather than an exhaustive catalog or study of statistical properties.The authors frame the framework as a flexible stage for modeling extensions.Source: §5 · p. 15

5.1 Sparse trend filtering

Sparse trend filtering allows a fitted trend to be exactly zero in some regions and vary smoothly elsewhere, using O(n) operations per iteration.The figure’s hand-tuned examples illustrate qualitative differences from ordinary trend filtering, not statistically fair comparisons.Source: §5.1 · p. 15, Figure 8 · p. 16

5.2 Mixed trend filtering

Mixed trend filtering lets either of two polynomial orders dominate at any location; with r mixed trends, each iteration takes O(nr) operations.The formulation generalizes to r polynomial orders.Source: §5.2 · p. 16

5.3 Trend filtering with outlier detection

Trend filtering with outlier detection estimates a trend alongside outliers, identified by nonzero components of ẑ, in O(n) operations per iteration.The paper points to an illustrative example in Figure 8.Source: §5.3 · pp. 16–17

5.4 Isotonic trend filtering

Isotonic trend filtering imposes a monotonicity constraint on fitted values, and its ADMM updates take O(n) time.The updates use an isotonic regression fit.Source: §5.4 · p. 17

5.5 Nearly-isotonic trend filtering

Nearly-isotonic trend filtering penalizes pointwise nonmonotonicities rather than enforcing strict monotonicity, with O(n) time per update round.Its regression subproblem can be computed by modifying the dynamic programming algorithm for the 1d fused lasso.Source: §5.5 · p. 17

6 Conclusion

The paper proposes a simple, specialized ADMM approach that uses exact solvers for the 0th-order fused lasso to solve higher-order trend filtering. It is fast and robust across a wide range of problem sizes and regularization regimes, and can be modified for related problems.

6 Conclusion

The algorithm uses extremely fast, exact solvers for the 0th-order fused lasso to solve higher-order problems.The approach targets trend filtering with k ≥ 1.Source: §6 · p. 18
The specialized ADMM is fast and robust across a wide range of problem sizes and regularization regimes.The paper contrasts its robustness with primal-dual interior point methods.Source: §6 · p. 18
The approach can be modified for related problems, and software is freely available in R and C.The R interface is trendfilter in glmgen.Source: §6 · p. 18

A.1 Algorithm details for the motivating example

The motivating example compares specialized ADMM with several alternatives for a linear trend-filtering fit. Despite ill-conditioning and competing methods with cheap iterations, specialized ADMM reaches a visually accurate fit in few iterations.

A.1 Algorithm details for the motivating example

Condition numbers grow roughly as n^k, and large values can make the criterion very flat around its optimum.The passage notes that this concern can arise even for moderate problem sizes.Source: Figure 9 · p. 19, §A.1 · p. 19
The motivating example uses n = 1,000 points, k = 1, and λ = 1,000, with an exact solution computed by the generalized lasso dual path algorithm.The tuning parameter represents a statistically reasonable regularization level in the example.Source: §A.1 · p. 20
After 10,000 iterations, dual proximal-gradient fits remain wiggly rather than piecewise linear, although acceleration improves the fit.Each iteration takes O(n) operations, and step sizes were hand-selected for the best possible performance.Source: §A.1 · p. 20
After 10,000 iterations, proximal gradient on the lasso formulation remains unsatisfactory, while its accelerated fit has rounded corners.The formulation uses the falling factorial basis matrix H, with matrix products computable in O(n) time.Source: §A.1 · p. 20
After 5,000 iterations, standard ADMM fails to capture piecewise-linear structure overall, while coordinate descent is far from the exact fit.Coordinate descent requires O(n^2) operations per iteration, but delivers a visually perfect piecewise-linear fit after nearly 100,000 iterations.Source: §A.1 · pp. 20–21
After 20 iterations, specialized ADMM and PDIP are visually indistinguishable from the exact solution; specialized ADMM is visually passable after 5.Both algorithms use O(n) operations per iteration.Source: §A.1 · p. 21

A.2 ADMM vs. PDIP for k = 3 (piecewise cubic fitting)

For piecewise cubic fitting, PDIP’s convergence difficulties appear at smaller problem sizes than in the k = 2 case. Specialized ADMM is somewhat slower to converge but remains fast and robust overall.

A.2 ADMM vs. PDIP for k = 3 (piecewise cubic fitting)

For k = 3 piecewise cubic fitting, PDIP’s behavior mirrors the k = 2 case.Its convergence issues begin at problem sizes smaller by an order of magnitude.Source: §A.2 · p. 21
PDIP convergence issues begin at problem sizes an order of magnitude smaller than in the k = 2 case.Source: §A.2 · p. 21
Specialized ADMM converges somewhat more slowly, yet remains fast and robust overall for k = 3.Figure 10 supports this comparison.Source: §A.2 · p. 21

A.3 Prediction at arbitrary points

Trend-filtering fitted values at observed inputs can be viewed as evaluations of an underlying function. The paper describes extending that function to arbitrary locations using falling-factorial basis functions and precomputed coefficients.

A.3 Prediction at arbitrary points

The fitted values at observed inputs are treated as evaluations of an underlying function that can be predicted at arbitrary domain locations.The software package includes this prediction task, although it is not central to the algorithmic proposal.Source: §A.3 · p. 21
The continuous fitted function is represented using falling-factorial basis functions and precomputed coefficients.The coefficients index polynomial and knot-producing basis functions.Source: §A.3 · pp. 21–22
Evaluating the fit at m new locations takes O(mr) operations, where r is the number of nonzero (k + 1)st-order differences.The paper treats k as a constant.Source: §A.3 · p. 22