PaperlayerView exact revision on arXiv

Source-linked AI summary

Efficient Implementations of the Generalized Lasso Dual Path Algorithm

Taylor Arnold, Ryan Tibshirani

arXiv:1405.3222v2stat.COcs.LGstat.ML

TL;DR

The paper develops efficient implementations of the generalized lasso dual path algorithm for general and structured penalty matrices. Specialized fused-lasso and trend-filtering strategies use structure to reduce per-iteration work, while approximate linear-system solutions can risk errors accumulating along the path.

  • Problem

    The generalized lasso dual path algorithm requires repeated linear-system solves, and its per-iteration cost can make computing the path burdensome.

    Source: §2.3 · p. 8, §3 · p. 10

  • Method

    The authors implement the dual path algorithm generically and exploit fused-lasso and trend-filtering structure with specialized solvers.

    Source: §1 · p. 1, §4 · p. 11, §4.2.2 · p. 12

  • Results

    The wide strategy costs O(rn) per subsequent iteration versus O(r^2n) for the naive strategy, and the trend-filtering algorithm requires O(r) operations.

    Source: §2.3 · p. 8, §3 · p. 10

  • Takeaways & Limitations

    The implementations cover general and structured penalties with X = I or general X, and are available in the genlasso R package.

    Source: §8 · p. 23

    Approximate linear-system solutions can accumulate errors along the path, so the proposed implementations use direct solvers.

    Source: §1.3.2 · p. 7

Abstract

from arXiv · show

We consider efficient implementations of the generalized lasso dual path algorithm of Tibshirani and Taylor (2011). We first describe a generic approach that covers any penalty matrix D and any (full column rank) matrix X of predictor variables. We then describe fast implementations for the special cases of trend filtering problems, fused lasso problems, and sparse fused lasso problems, both with X=I and a general matrix X. These specialized implementations offer a considerable improvement over the generic implementation, both in terms of numerical stability and efficiency of the solution path computation. These algorithms are all available for use in the genlasso R package, which can be found in the CRAN repository.

1 Introduction

The paper studies the generalized lasso under a full-column-rank predictor matrix, deriving efficient implementations of its dual path algorithm. It connects established penalty structures to a unified path framework and develops generic and specialized computational strategies.

1 Introduction

  • The generalized lasso includes the standard lasso when D = I, and assuming X has full column rank ensures a unique solution for every λ..Source: §1 · p. 1
  • The authors derive efficient implementations of the generalized lasso dual path algorithm, which computes the solution over the full range of regularization values.They provide a generic implementation and specialized versions for trend filtering and fused lasso problems.Source: §1 · p. 1
  • For trend filtering, the order-(k+1) difference penalty yields piecewise kth-degree polynomial solutions; the 1d fused lasso is its k=0 special case.The difference matrix is banded with bandwidth k + 2.Source: §1 · p. 2

1.1 Related work

  • Prior methods solve either fixed-λ problems or solution paths; the dual path framework permits general D while also supporting specialized implementations.The authors focus on stable, efficient implementations rather than detailed comparisons with alternative methods.Source: §1.1 · p. 3

1.2 The dual path algorithm

  • The dual path algorithm tracks dual coordinates on the constraint boundary as λ decreases, recording when coordinates enter or leave to recover the primal path.For general X, a transformation of the outcome and penalty matrix reduces the dual problem to the same form as the signal-approximator case.Source: §1.2.1 · p. 4, §1.2.2 · p. 5

1.2.1 The signal approximator case, X = I

  • With X = I, the algorithm’s main computational work is repeated least-squares solves as the boundary set changes; the dual solution need not be unique.The path proceeds from λ = ∞ toward λ = 0.Source: §1.2.1 · p. 4

1.2.2 The general X case

  • For general X, the dual can be transformed using the pseudoinverse X+, after which the signal-approximator algorithm applies; a small ridge term can handle rank-deficient X.The paper assumes full column rank for its main algorithm.Source: §1.2.2 · pp. 4–5

1.3 Implementation overview

  • Updating a QR decomposition reduces generic per-iteration work, while trend-filtering and fused-lasso structure enables faster banded or Laplacian solves.For general X, transforming D can destroy that structure; specialized methods nevertheless avoid an initial pseudoinverse and require solving a full system in X^T X only at the path’s end. Total path steps are not precisely understood, so the analysis focuses on per-iteration complexity.Source: §1.3.1 · p. 5, §1.3.2 · p. 6

2 QR-based implementation for a general D

The general-D implementation uses QR decompositions to solve and update the least-squares problems along the dual path. A wide strategy is more efficient when D has full row rank and m ≤ n, while the tall strategy applies regardless of D’s dimensions and rank.

2 QR-based implementation for a general D

  • For general D, the implementation updates QR decompositions as rows are added to or removed from D−B, reducing computational work by essentially an order of magnitude versus recomputing them.A general full-column-rank X adds O(np^2) operations to compute X+, after which the QR approach applies to the transformed penalty matrix.Source: §2 · p. 7

2.1 The wide strategy

  • The wide strategy maintains a QR decomposition of D^T and applies when m ≤ n and rank(D) = m.It checks the rank through the triangular R factor; its initial decomposition costs O(m^2n), and subsequent iterations cost O(rn).Source: §2 · p. 7, §2.1 · pp. 7–8

2.2 The tall strategy

  • The tall strategy maintains a rotated QR decomposition of D and handles any dimensions or rank of D.Its initial work costs O(mn · min{m, n}); updates cost O(max{r^2, n^2}) per iteration.Source: §2 · p. 7, §2.2 · p. 8, §2.3 · p. 8

2.3 Computational complexity comparisons

  • For the tall strategy, the advantage over the naive approach is largest when r = Θ(n); near the path’s end, switching to the naive strategy can be more efficient.The wide strategy uses O(rn) operations per iteration versus O(r^2n) for the naive strategy.Source: §2.3 · pp. 8–9

3 Specialized implementation for trend filtering, X = I

For trend filtering with X = I, the banded structure of the discrete difference operators enables fast dual-path computations. Sparse QR can offer better numerical accuracy, while neither proposed decomposition approach updates its factorization across iterations.

3 Specialized implementation for trend filtering, X = I

  • With fixed k, each trend-filtering path iteration takes O(r) operations, linear in the number of interior coordinates.The banded linear systems have bandwidth 2k + 3 and can be solved in O(rk^2) operations.Source: §3 · pp. 9–10
  • Sparse QR gives better numerical accuracy than the banded Cholesky approach because it operates on a matrix whose condition number is squared relative to D(k+1)−B.The genlasso package uses a sparse QR algorithm.Source: §3 · p. 10
  • Neither banded Cholesky nor sparse QR updates a single matrix decomposition across iterations; sparsity-preserving update schemes remain future work.Dense update rules can create fill-in when used with sparse decompositions.Source: §3 · p. 10

4 Specialized implementation for the fused lasso, X = I

For fused lasso with X = I, the algorithm replaces difficult minimum-norm solves with graph-structured projections and Laplacian systems. The implementation tracks changing connectivity and extends to sparse fused lasso by incorporating coefficient penalties.

4 Specialized implementation for the fused lasso, X = I

  • The key obstacle is computing minimum ℓ2 norm solutions of sparse, underdetermined systems, which is harder than finding arbitrary basic solutions.The specialized strategy targets these solves in the dual path steps.Source: §4 · p. 10

4.1 Alternative form for Steps 1 and 3(a) in Algorithm 1

  • An alternate formulation obtains the minimum ℓ2 norm solution by projecting the right-hand side onto row(D), solving D^T D z = d, and setting x = Dz.For fused lasso, the needed null-space projections can be computed by averaging over connected components.Source: §4.1 · p. 11, §4 · p. 11

4.2.1 Null space of the oriented incidence matrix

  • For an oriented graph incidence matrix, null(D) is spanned by connected-component indicators, so projection onto it is componentwise averaging.After edges in B are removed, the same averaging applies within components of G−B.Source: §4.2.1 · p. 12

4.2.2 Solving Laplacian linear systems

  • Because D^T D is the graph Laplacian, its systems can be split by connected component and solved after removing one row and column per connected subsystem.The reduced system is sparse and nonsingular, permitting a direct sparse Cholesky solve; approximate iterative solvers can cause practical path-algorithm issues.Source: §4.2.2 · pp. 12–13

4.2.3 Tracking graph connectivity across iterations

  • Across iterations, G−B changes by one edge, so a breadth-first or depth-first search from an incident node checks whether connectivity changed.With sparse Cholesky on reduced systems, the reported per-iteration operation count is approximately O(max{r, n}), where r = m − |B| and n is the node count.Source: §4.2.3 · p. 14

4.3 Extension to sparse fused lasso problems

  • Sparse fused lasso appends a constant multiple of the identity to the graph incidence penalty matrix, yielding componentwise systems that are either Laplacian or nonsingular.The nonsingular systems can be factored directly; for X = I, the sparse fused lasso path can also be obtained by soft-thresholding the fused lasso path.Source: §4.3 · pp. 14–15

5 Specialized implementations with a general X

With a general full-column-rank X, directly applying the path algorithm to the transformed penalty can lose the efficiency of structured penalties. The paper instead uses null-space bases and structured solves, with specialized bases for trend filtering and graph penalties.

5 Specialized implementations with a general X

  • The transformed penalty DX+ need not retain D's structure, so directly applying the usual path algorithm can substantially reduce relative efficiency.The proposed approach constructs solutions to the relevant least-squares problems without relying on that structure being preserved.Source: §5 · p. 16

5.1 Alternate form of computations in Algorithm 2

  • When null-space bases and systems in D or D−B are easy to compute, the method avoids generic dense solves in the transformed penalty and avoids explicit X+ calculations.Its projections require systems in (XH)^T XH, whose dimension is the nullity q; q is typically small at early path stages of interest.Source: §5.1 · pp. 17–18

5.2 Trend filtering, general X

  • For order-k trend filtering, null(D) has dimension k + 1 and is spanned by vectors formed from repeated cumulative sums.When rows indexed by B are removed, null(D−B) has dimension k + 1 + d, where d = |B|, with an extended basis.Source: §5.2 · p. 18

5.3 Fused lasso and sparse fused lasso, general X

  • For fused lasso, null(D−B) is spanned by indicators of the connected components of G−B; for sparse fused lasso, padded component indicators provide a basis.These bases feed the structured Laplacian-based solves used in the general-X procedure.Source: §5.3 · pp. 18–19

6 Chicago crime data example

Using reported burglaries per household across Chicago census blocks, the example applies the fused lasso to estimate burglary risk and group adjacent blocks. The resulting solution shows broad spatial regions, while fused lasso's local adaptivity lets cluster size respond to the measurements.

6 Chicago crime data example

  • The measurements are reported burglaries per household for 2005–2009, aggregated within Chicago's 2010 census block groups.They are treated as noisy measurements of underlying household burglary probabilities.Source: §6 · p. 19
  • The fused lasso groups adjacent census blocks by penalizing differences between neighboring block estimates, allowing exact fusions that smooth penalties cannot produce.The example uses 2,167 blocks and 14,060 neighboring-block connections, with X = I; the first 2,500 path steps took a little over a minute on a laptop.Source: §6 · p. 19
  • The 8-degree-of-freedom solution divides Chicago into roughly four regions, with highest risk in the south, lowest in the north, and a small isolated low-risk area.The isolated area is buffered by blocks with no data, so it incurs only a small penalty for separating from the main graph.Source: §6 · p. 19
  • Fused lasso is locally adaptive: it chooses cluster size from nodewise measurements rather than tending toward roughly equal-sized clusters.Source: §6 · p. 19

7 Empirical timings

The paper measures runtime for the first 100 path steps across three problem classes and reports near-linear empirical scaling. The reported runtimes are for default R 3.1, and the authors note that many useful grid-graph solutions may require prohibitively many steps.

7 Empirical timings

  • The authors time the first 100 steps because full paths can require super-linear iteration counts, and many applications need the more-regularized solutions near the path's start.For 1d fused lasso with X = I, the number of steps to termination is an exception and is always n − 1.Source: §7 · p. 22
  • Empirical scaling is approximately linear for 1d fused lasso and cubic trend filtering, and slightly slower but still near-linear for 2d fused lasso; the authors project 100 steps on problems in the millions to take under an hour.Figure 3 reports first-100-step runtimes for n = 1,000 to 50,000 and presents the three problem classes in separate panels.Source: Figure 3 · p. 22, §7 · p. 22
  • Interesting solutions on geographic-region graphs may require fewer steps, while grid graphs can require prohibitively many, especially for image denoising with many connected components.Source: §7 · p. 22
  • The runtime measurements were obtained with default R, specifically version 3.1.Source: §7 · p. 23

8 Discussion

The discussion summarizes efficient generalized-lasso dual-path implementations for general penalties, trend filtering, and fused lasso, covering both X = I and general predictor matrices.

8 Discussion

  • The implementations cover general penalty matrices, discrete-difference trend-filtering penalties, and graph-incidence fused-lasso penalties, for both X = I and general X, and are included in genlasso.Source: §8 · p. 23

A The QR decomposition and least squares problems

The appendix briefly reviews QR decomposition and its application to least-squares problems.

A The QR decomposition and least squares problems

A.1 The QR decomposition of a full column rank matrix

The appendix develops QR decompositions for full-rank and rank-deficient least-squares problems, including minimum-norm solutions. It also explains how Givens rotations create and maintain the triangular structure used in these decompositions.

A.1 The QR decomposition of a full column rank matrix

  • After computing a QR decomposition, a full-column-rank least-squares solution is found by forming c = Q1^T b and back-solving R1x = c in O(mn) operations.The initial QR decomposition costs O(mn^2); it can be reused to solve for multiple right-hand sides.Source: §A.1 · pp. 23–24

A.2 The QR decomposition of a column rank deficient matrix

  • For rank-deficient A, a basic least-squares solution sets the free coordinates to zero and back-solves for the remaining coordinates, costing O(mn).The least-squares problem can have infinitely many solutions when rank(A) < n.Source: §A.2 · p. 24

A.3 The minimum ℓ2 norm least squares solution

  • The minimum-ℓ2-norm least-squares solution uses Givens rotations to transform the rank-deficient QR factorization before solving for the required coordinates.The stated total cost is O(n · max{m, n}).Source: §A.3 · pp. 24–25

A.4 The minimum ℓ2 norm least squares solution and the transposed QR

  • A QR decomposition of A^T can also be used to compute the minimum-ℓ2-norm least-squares solution, with total cost O(m · max{m, n}).The method uses an orthogonal change of variables and forward-solves for the nonzero coordinates.Source: §A.4 · p. 25

B Givens rotations

  • Givens rotations are orthogonal transformations used to create or maintain upper-triangular structure, including when a matrix gains or loses a row or column.A single rotation can be applied to a matrix in O(m) operations when it affects two columns, or O(n) when it affects two rows.Source: B Givens rotations · p. 26, §B.2 · p. 27

B.1 Simple Givens rotations in two dimensions

  • In two dimensions, a Givens rotation can rotate a vector onto the first coordinate axis without computing the angle explicitly.The rotation parameters are chosen directly from the vector components.Source: §B.1 · p. 26

B.2 Givens rotations in higher dimensions

  • In higher dimensions, a Givens rotation changes only two components and can zero a selected entry while preserving existing zeros.Applying rotations to rows or columns can build an upper-triangular structure.Source: §B.2 · pp. 26–27

C Updating the QR decomposition in the full rank case

The section updates an existing QR decomposition when a matrix changes by one row or one column. Updating the decomposition saves an order of magnitude in computational time compared with recomputing it from scratch.

C Updating the QR decomposition in the full rank case

  • Updating the QR decomposition after a one-row or one-column change saves an order of magnitude in computational time over computing it from scratch.The section treats row and column updates separately.Source: C Updating the QR decomposition in the… · p. 27

C.1 Adding or removing a row

The section uses Givens rotations to update a full-rank QR decomposition when a row is added or removed. For a row addition, the updated factorization is obtained with O(mn) operations.

C.1 Adding or removing a row

  • Adding a row uses Givens rotations to obtain the updated QR decomposition in O(mn) operations.The described procedure applies n rotations; the passage also gives a specific count of m − 1 rotations for a row-addition construction.Source: §C.1 · pp. 28–29
  • Removing a row uses Givens rotations to transform the factorization into the QR decomposition of the matrix with that row deleted.The derivation assumes m > n so that removing a row does not change the rank.Source: §C.1 · pp. 28–29

C.2 Adding or removing a column

For full-rank matrices, adding or removing a column updates the QR decomposition in O(mn) operations. The rank-deficient treatment handles row changes directly and column changes through a QR decomposition of the transpose.

C.2 Adding or removing a column

  • In the full-rank case, adding or removing a column updates the QR decomposition in O(mn) operations.The addition assumes the new column is outside the span of the existing columns.Source: §C.2 · p. 29

D Updating the QR decomposition in the rank deficient case

  • For rank-deficient matrices, the update problem is to compute minimum-ℓ2-norm least-squares solutions after a row or column is added or removed.The initial QR decomposition must be chosen to suit whether the update is to rows or columns.Source: D Updating the QR decomposition in the… · pp. 29–30

D.1 Adding or removing a row

  • Rank-deficient row additions and removals use Givens rotations, with costs O(n · max{m, n}) and O(m · max{m, n}), respectively.The addition procedure applies O(n) rotations; the removal procedure restores the required QR form when rank decreases.Source: §D.1 · pp. 31–32

D.2 Adding or removing a column

  • For rank-deficient column changes, the method updates a QR decomposition of A^T; addition costs O(m · max{m, n}) and removal costs O(n · max{m, n}).The transposed decomposition is first computed for the initial least-squares problem and then updated for the changed matrix.Source: §D.2 · p. 32

E More plots from the Chicago crime data example

The section presents five additional fused lasso solutions along the path for the Chicago crime data example, at different degrees of freedom.

E More plots from the Chicago crime data example

Loading 1405.3222v2…