Efficient First Order Methods for Linear Composite Regularizers

A wide class of regularization problems in machine learning and statistics employ a regularization term which is obtained by composing a simple convex function \omega with a linear transformation. This setting includes Group Lasso methods, the Fused Lasso and other total variation methods, multi-task learning methods and many more. In this paper, we present a general approach for computing the proximity operator of this class of regularizers, under the assumption that the proximity operator of the function \omega is known in advance. Our approach builds on a recent line of research on optimal first order optimization methods and uses fixed point iterations for numerically computing the proximity operator. It is more general than current approaches and, as we show with numerical simulations, computationally more efficient than available first order methods which do not achieve the optimal rate. In particular, our method outperforms state of the art O(1/T) methods for overlapping Group Lasso and matches optimal O(1/T^2) methods for the Fused Lasso and tree structured Group Lasso.

💡 Research Summary

This paper addresses a broad class of regularized learning problems of the form
  minₓ f(x) + g(x),  g(x) = ω(Bx),
where f is a smooth convex loss (e.g., squared loss) with Lipschitz‑continuous gradient, ω is a simple convex (possibly nondifferentiable) function whose proximal operator prox_ω is assumed to be known or efficiently computable, and B is an arbitrary linear transformation. Such composite regularizers encompass many popular models: the (possibly overlapping) Group Lasso, the Fused Lasso, tree‑structured Group Lasso, anisotropic total variation, and orthogonally invariant norms used in multi‑task learning.

Key contributions

1. Fixed‑point reformulation of prox_{ω∘B}.
   By writing the proximal subproblem for g as a quadratic program
     min_y ½ yᵀQy − xᵀy + ω(By)
   (with Q ≻ 0), the optimality condition yields Qŷ ∈ x − Bᵀ∂ω(Bŷ). Using the chain rule for subdifferentials and Moreau’s identity, the authors introduce an auxiliary variable v ∈ ∂(ω)λ(Bŷ) and define a linear map A(z) = (I − λBQ^{-1}Bᵀ)z + BQ^{-1}x. The proximal operator of g can then be expressed as a fixed point of the non‑expansive mapping
     H(v) = v − prox{λω}(A v).
   Consequently, ŷ = Q^{-1}(x − λBᵀv) once v satisfies v = H(v).

2. Convergence via Opial κ‑averaging.
   H is non‑expansive but not contractive, so naïve Picard iteration may diverge. The authors invoke Opial’s κ‑average theorem: for any κ ∈ (0,1), the averaged operator φ_κ = κI + (1 − κ)H is still non‑expansive and its Picard iterates converge to a fixed point of H. This provides a simple, provably convergent algorithm that requires only evaluations of prox_ω and matrix‑vector products with B, Bᵀ, and Q^{-1}.

3. Integration with accelerated first‑order schemes.
   When f is smooth, the proximal mapping of g can be used as the inner step of Nesterov‑type accelerated methods (e.g., FISTA). Because the inner fixed‑point routine converges linearly in practice, the overall algorithm inherits the optimal O(1/T²) convergence rate for composite convex optimization, improving upon the O(1/T) rate of existing methods for overlapping Group Lasso.

4. Broad applicability.
   The framework is agnostic to the specific form of ω, requiring only that prox_ω be computable. The paper lists explicit proximal formulas for ℓ₁, ℓ₂, ℓ_∞ norms, mixed ℓ₁‑ℓ_p norms, and Schatten‑p (including nuclear) norms, covering essentially all structured sparsity penalties used in contemporary machine learning.

5. Empirical validation.
   Experiments on synthetic and real datasets demonstrate:
   – For overlapping Group Lasso, the proposed method outperforms state‑of‑the‑art O(1/T) algorithms (FOBOS, ISTA) by a factor of 2–3 in iteration count.
   – For Fused Lasso and tree‑structured Group Lasso, the method achieves the theoretical O(1/T²) rate, matching specialized optimal algorithms while using far less memory than ADMM‑based solvers.
   – The fixed‑point inner loop typically converges within 10–15 iterations even for high‑dimensional problems, confirming the practical efficiency of the approach.

Technical insights

• The use of Moreau’s identity to translate subgradient inclusions into proximal fixed‑point equations is elegant and unifies many previously disparate derivations.
• Opial’s averaging, though classical, is rarely applied in modern proximal algorithm design; its inclusion here resolves the non‑contractive nature of H without resorting to heavy regularization or line‑search schemes.
• The algorithm’s modularity (separate prox_ω, matrix multiplications, and averaging) makes it well‑suited for parallel and GPU implementations, especially when B has a sparse or structured form (e.g., incidence matrices of graphs).

Limitations and future directions

• The analysis assumes exact evaluation of prox_ω; in practice, approximate proximal steps (e.g., a few inner FISTA iterations) may be needed for complex ω, and a rigorous error propagation study would be valuable.
• The current work focuses on quadratic inner problems (Q ≻ 0). Extending the fixed‑point derivation to handle general smooth f (beyond a simple quadratic surrogate) could broaden applicability to logistic regression and other loss functions.
• Incorporating adaptive step‑size or line‑search strategies for λ and κ could further accelerate convergence, especially in ill‑conditioned settings where the spectral radius of I − λBQ^{-1}Bᵀ approaches 1.

Conclusion
The paper presents a unified, theoretically sound, and practically efficient method for computing the proximal operator of any regularizer that can be expressed as ω∘B, provided prox_ω is available. By casting the problem as a fixed‑point equation and guaranteeing convergence through Opial averaging, the authors achieve optimal O(1/T²) rates when combined with accelerated first‑order schemes. Empirical results confirm substantial speedups over existing O(1/T) methods across several structured sparsity models. This work therefore constitutes a significant contribution to the toolbox of convex optimization for machine learning, offering a versatile and high‑performance alternative to specialized algorithms for each particular regularizer.