Modular proximal optimization for multidimensional total-variation regularization

1 Introduction

• Prox-operators localize the complexity of 1D and 2D TV regularization within a broader sparsity-inducing optimization problem.The paper focuses on weighted anisotropic TV, a structured-sparsity regularizer.Source: §1 · pp. 1–2

1.1 Contributions

• Geometric analysis yields an efficient taut-string variant and a hybrid connected to Condat’s method, including a new weighted ℓ1-norm 1D-TV solver.The authors report that their taut-string version performs better than most recently developed TV proximity methods.Source: §1.1 · p. 2
• For ℓp-TV beyond ℓ1, the paper presents a specialized Newton method for p = 2 and projection-free and projection-based first-order approaches for general p ≥ 1.The introduction also describes the work as reviewing, analyzing, implementing, and experimentally comparing fast algorithms.Source: §1.1 · p. 2
• For multidimensional TV, the paper reviews proximal-splitting methods and reports alternating reflections as most effective in practice, with modular reuse of 1D solvers and parallelization.Source: §1.1 · p. 2
• The paper provides multi-threaded C++, Matlab, and Python implementations of the reviewed and developed methods, with code, datasets, and scripts intended to support independent verification.Source: §1.1 · p. 3

1.2 Related work

• Anisotropic TV appears in denoising, deconvolution, fused-lasso, and change-point detection,,; fast prox-operators also benefit higher-dimensional TV and trend-filtering solvers, [87].Source: §1.2 · p. 3

1.3 Summary of the paper

• The paper moves from 1D ℓ1-TV and ℓp-TV proximity to multidimensional splitting methods, then presents experiments and applications.Source: §1.3 · pp. 3–4

2 TV-L1: Fast prox-operators for Tv1D

• The section focuses on carefully tuned 1D TV-L1 prox-solvers as fast, robust, low-memory building blocks for higher-dimensional TV problems.It introduces taut-string methods for unweighted TV-L1 before discussing elementwise weighted TV, which is important in applications such as segmentation and submodular optimization [43].Source: §2 · p. 4

2.1 The taut-string method for Tv1D

• The classic taut-string method computes the TV-L1 solution as the discrete gradient of a shortest polyline constrained to a λ-wide tube around the cumulative input.Each tube segment enters and leaves the convex minorant or concave majorant once, giving an overall O(n) runtime.Source: §2.1 · pp. 5–6

2.1.1 Efficient implementation of taut-strings

• A fixed-length array deque makes classic taut-string majorant and minorant operations constant-time, with one initial allocation and no later memory management.The queue tracks valid segments with pointers; the algorithm also precomputes segment slopes for repeated comparisons.Source: §2.1.1 · pp. 6–7

2.2 Linearized taut-string method for Tv1D

• The linearized taut-string method replaces stored convex and concave bounds with affine approximations and is equivalent, nearly line by line, to Condat’s algorithm [28].It uses a constant number of bookkeeping variables and supports in-place updates, but restarting at breakpoints permits O(n^2) worst-case runtime.Source: §2.2 · pp. 8–10

2.3 Comparison of taut-string methods and a hybrid strategy

• The classic method guarantees linear runtime, while the linearized method equivalent to Condat uses less memory and has faster per-step traversal but can take O(n^2) in the worst case.The authors report that a deque-based classic implementation can be very efficient in practice, while reserving empirical review of the linearized method’s practical runtime for later sections.Source: §2.3 · p. 11
• Capping the linearized solver at n^S steps (1 < S < 2) before switching to the classic solver gives a worst-case O(n^S) bound while retaining the classic method’s O(n) guarantee.The switch occurs if the linearized method has not found the solution by the cap; the classic solver handles the remaining signal.Source: §2.3 · p. 12

2.4 Taut-string methods for weighted Tv1D

• Weighted TV-L1 can be solved by taut-string methods using per-component bounds |s_i − r_i| ≤ w_i, which create a variable-width tube rather than the unweighted fixed-width tube.The classic method sets the tube limits to r_i − w_i and r_i + w_i.Source: §2.4 · p. 12
• The linearized taut-string method extends to weighted TV by adapting its tube-violation checks, slope recomputations, and restarts; its equivalence with the classic method still holds.The weighted variant retains the benefits and limitations described for the unweighted solvers.Source: §2.4 · p. 12

3.1 TV-L2: Proximity for Tv1D

• A hybrid MSN–GP approach is recommended for TV-L2, with GP competitive only over a limited range of λ values.The specialized MSN method exploits the tridiagonal Hessian structure.Source: §3.1 · pp. 13, 15
• The TV-L2 MSN iteration runs in O(n) time by exploiting the tridiagonal structure of DDT + αI.Source: §3.1 · pp. 13–14

3.2 TV-Lp: Proximity for Tv1D

• For 1 < p < ∞, the dual problem favors a hybrid of GP and FW because Newton-type methods do not apply easily and the GP projection is involved.Source: §3.2 · p. 15

3.2.1 Efficient projection onto the ℓq-ball

• The projected-Newton direction for the ℓq-ball projection can be computed in O(n) time by exploiting the reduced-Hessian structure.Source: §3.2.1 · p. 16

3.2.2 Frank-Wolfe algorithm for TV-Lp proximity

• Frank–Wolfe obtains a descent direction by solving a linear problem over the constraint set and uses a closed-form stepsize for the quadratic objective.Source: §3.2.2 · pp. 16–17

3.3 Prox operator for TV-L∞

• TV-L∞ proximity uses GP, whose nontrivial step is projection onto the ℓ1-ball, available through O(n)-time routines.Source: §3.3 · p. 17

4.1 Proximity stacking

• Proximity stacking combines individual regularizer prox operators to compute the prox for their sum while preserving modularity, enabling a framework for convex machine-learning problems with multiple regularizers.Source: §4.1 · p. 18

4.1.1 Proximal Dykstra (PD)

• Proximal Dykstra combines two regularizers, while its parallel generalization handles more than two and permits parallel prox calls.Source: §4.1.1 · pp. 18–19

4.1.2 Alternating reflections – Douglas-Rachford (DR)

• The alternating-reflections Douglas–Rachford variant computes the solution through proximity steps and does not require hyperparameter tuning.Source: §4.1.2 · p. 20

4.1.3 Alternating-Direction Method of Multipliers (ADMM)

• ADMM applies consensus constraints and alternates minimization over x and each z_i with gradient steps over each u_i; its prox calls can run in parallel.Source: §4.1.3 · p. 21

4.1.4 Dual proximity methods

• Dual proximity methods can be advantageous when dual prox operators are easier to compute than primal ones; isotropic TV is given as an instance motivating anisotropic variants.Source: §4.1.4 · p. 21

4.2 Two-dimensional TV

• Anisotropic 2D-TV splits into row and column regularizers, each composed of independent 1D-TV terms on distinct rows or columns.Source: §4.2 · p. 24

4.3 Higher-dimensional TV

• Higher-dimensional TV decomposes into dimension-wise sums of 1D-fiber TV terms, which can be combined through the m-term strategies and parallelized.Source: §4.3 · p. 24

• The experiments examine 1D-TV methods, other one-dimensional variants, and multidimensional TV.Source: §5 · p. 24
• The solvers are implemented in C++ and distributed through proxTV with Matlab and Python bindings.Source: §5 · p. 24
• The library’s matrix operations use LAPACK.Source: §5 · p. 24

5.1.1 Running time results for synthetic data

• On synthetic signals, direct methods are much faster than iterative methods, and Condat’s and linearized taut-string methods show practical O(n) behavior despite a theoretical O(n^2) worst case.Johnson’s and Kolmogorov’s methods have guaranteed O(n) runtime but are slower than those methods in these experiments.Source: §5.1.1 · p. 26

5.1.2 Worst case scenario

• The proposed hybrid avoids quadratic runtimes in a constructed worst-case signal where Condat’s and linearized taut-string methods perform much worse than other direct methods.The experiment was designed from Condat’s worst-case description and tested increasing signal lengths.Source: §5.1.2 · p. 26

5.1.3 Running times on natural images

• Across natural-image runs, the classic taut-string method tends to perform best, while the hybrid tracks Condat at low penalties and Kolmogorov at large penalties.Condat’s method performs well at low penalties but scales poorly as penalties increase; the classic taut-string, Johnson, and Kolmogorov methods degrade less.Source: §5.1.3 · pp. 26–27

5.1.4 Running time results for weighted TV-L1

• Weighted and unweighted TV-L1 problems have almost identical performance for the tested methods, although adapting other reviewed algorithms to weighted problems while maintaining performance is not straightforward.The comparison includes Projected Newton and linearized taut string on both standard and weighted TV-L1 proximity problems.Source: §5.1.4 · p. 28

5.2.1 Running time results for TV-L2

• For TV-L2, MSN meets the stopping criterion across tested settings, whereas GP can fail after 10,000 iterations at moderately large penalties; the hybrid mimics the better method’s behavior.The hybrid uses GP when λ < ∥y∥₂ and otherwise MSN, switching from GP to MSN if it fails within 50 iterations.Source: §5.2.1 · pp. 29–30

5.2.2 Running time results for TV-Lp

• For TV-Lp, the Projected Newton projection is more efficient than SLEP’s epp within GP, while Frank–Wolfe is significantly faster than Projected Newton; a GP+FW hybrid attains the target duality gap faster than GP.The projection subroutine must be accurate to about 10^-12 in duality-gap terms for GP to meet the desired stopping threshold.Source: §5.2.2 · pp. 31–32

5.2.3 Running time results for TV-L∞

• For TV-L∞, GP with a standard ℓ1-projection subroutine converges faster than the general TV-Lp GP solver.The reported explanation is that ℓ1 projection is easier than general ℓq projection.Source: §5.2.3 · p. 32

5.2.4 Application: Proximal optimization for Fused-Lasso

• The fused-lasso variants can be solved by plugging their gradients and proximal operators into proximal splitting, and the combined ℓ1-plus-TV-L1 prox can be computed directly for FL and LFL.That direct decomposition removes the need for a combiner algorithm for those two models; the Variable Fusion variant is described as new and uses a more conservative TVp penalty that does not oversmooth estimates.Source: §5.2.4 · p. 33

5.2.5 Fused-Lasso experiments: simulation

• Replacing SLEP’s default FLSA subroutine with the fast taut-string TV-L1 solver speeds optimization when synthetic data matrices have very large column dimensions, but yields negligible gains for more balanced matrices.For large row counts, gradient computation dominates the cost and limits the benefit of faster proximity computations.Source: §5.2.5 · p. 34

5.2.6 Fused-Lasso Experiments: Microarray classification

• On microarray classification, logistic-loss fused-lasso models generally yield better accuracies than least-squares models, except on Ovarian; TV-L1 generally performs better than TV-L2, which is competitive in some settings.The reported comparison covers ArrayCGH, Leukemias, Colon, Ovarian, and Rat datasets.Source: §5.2.6 · pp. 34–35

5.3.1 Image denoising through anisotropic filtering

• Condat and Chambolle-Pock perform similarly and slowly; Split Bregman performs significantly better, while Goldfarb-Yin is slower and appears suboptimal.The authors attribute Split Bregman’s performance to its implementation and fast approximate matrix inversions, and Goldfarb-Yin’s suboptimality to its discrete approximation.Source: §5.3.1 · p. 37
• Mid-quality solutions suffice for denoising as ISNR saturates rapidly; Douglas-Rachford is faster than Proximal Dykstra, while Kolmogorov and ADMM suit high accuracy.The authors conclude their proposed methods are best suited to the considered image-processing task.Source: §5.3.1 · pp. 37–38, Figure 19 · p. 39

5.3.2 Parallelization experiments

• Parallel speedups are substantial, though gains flatten near eight cores for smaller images; Split Bregman is harder to parallelize.The authors attribute a speedup bump at 14 processors to a quirk of their multicore machine.Source: §5.3.2 · p. 40, Figure 20 · p. 41

5.3.3 Anisotropic image deconvolution

• In deconvolution, Douglas-Rachford is slightly better at finding accurate solutions and faster to converge to the final ISNR than Kolmogorov et al.The comparison uses five Douglas-Rachford or ten Kolmogorov inner iterations as the configurations that reach relative distance below 10^-3.Source: §5.3.3 · p. 42, Figure 21 · p. 43, Figure 22 · p. 44

5.3.4 2D Fused-Lasso Signal Approximator

• For 2D fused-lasso signal approximation, the solver’s performance depends on 2D-TV proximity, and Douglas-Rachford performs best when reconstruction ISNR is the primary concern.The method decomposes the problem into a 2D-TV prox operator followed by closed-form soft-thresholding.Source: §5.3.4 · pp. 42–43

5.4 Application of higher-dimensional TV

• The video filter computes a 3D-TV prox on a tensor whose dimensions represent image height, width, and frames, using Parallel Proximal Dykstra.The experiments use grayscale sequences with zero-mean Gaussian noise of variance 0.01.Source: §5.4 · pp. 44–45
• Parallel Proximal Dykstra is fastest to reach mid-quality solutions, while ADMM is preferable for high accuracy; both converge to equivalent solutions in practice.Goldfarb-Yin again provides suboptimal solutions because it uses a discrete approximation.Source: Figure 24 · p. 46, §5.4 · p. 46
• Intermediate PPD solutions yield better ISNR than convergence on the coastguard and bicycle sequences, so the authors claim efficiency improvements, not novel ISNR quality.They note that more advanced denoising models exist for these signal types.Source: §5.4 · p. 46

Prox operators and primal-dual relations

• A prox operator maps an input to the unique point defined by the proximity problem, and Moreau decomposition relates a function’s prox to its conjugate’s prox.The appendix says these primal-dual relations provide tools used in the algorithms.Source: A Mathematical background · p. 51

Submodular functions and Lovász extensions

• Submodular functions have diminishing returns, and their Lovász extensions are convex with support functions equal to those of their base polytopes.The appendix also states that submodular minimization can be cast as convex optimization.Source: A Mathematical background · pp. 52–53
• A base polytope is specified by inequalities over subsets and an equality for the complete set, but the paper notes it is not used directly.The constraints are combinatorial in number.Source: A Mathematical background · p. 52

Geometric formulation and projections

• For two functions, the resulting geometric problem is finding closest points between convex sets; projections onto base polytopes can be computed through prox operators.Projection onto the negated base polytope follows by negating the input and output of the corresponding projection.Source: A Mathematical background · pp. 53–54

B proxTV toolbox

• The proxTV toolbox provides TV proximity solvers in C++, Matlab, and Python, with a user-friendly top-level Matlab function.The example function TV accepts a signal and regularization value.Source: B proxTV toolbox · p. 54
• The TV function handles weighted 1D proximity by accepting a vector of weights as its regularization parameter.It also solves the unweighted problem for signals of any dimension.Source: B proxTV toolbox · p. 54
• Optional norm and regularizer settings allow different choices by signal dimension or combinations of TV regularizers along one dimension.Providing a third parameter selects the general Tvp proximity problem and the library chooses a solver.Source: B proxTV toolbox · p. 54

C Proof on the equality of taut-string problems

• For a nonzero vector w, problems C.1 and C.2 share the same minimum s*=ŝ*.The proof establishes optimality in both directions by transferring the KKT conditions between the problems.Source: C Proof on the equality of taut-string… · pp. 54, 56
• The two problems differ only in their objective functions, so their KKT conditions have the same form apart from the first condition.The proof uses this correspondence to transfer feasibility and complementary-slackness conditions between the formulations.Source: C Proof on the equality of taut-string… · p. 55
• Because f is contractive, the proof satisfies condition C.8 by setting (ˆαi, ˆβi) = (kαi, kβi) for some 0 ≤ k < 1; this works because αi and βi cannot both be zero.Source: C Proof on the equality of taut-string… · p. 55

D Proof on the equivalence of linearized taut-string method

• Affine approximations to the greatest convex minorant and smallest concave majorant do not change the taut-string solution.The proof establishes this by showing that the affine and convex/concave versions detect overlaps equivalently.Source: D Proof on the equivalence of… · pp. 56–57
• An overlap occurs when the greatest convex minorant of f+λ is no higher than the smallest concave majorant of f−λ at some point.The proof compares this condition with overlap detected by the affine approximations.Source: D Proof on the equivalence of… · p. 56
• The common origin clamp in the taut-string method is key to proving that affine overlap detection also holds for the convex/concave counterparts.The clamp is at the start of the method or at the last fixed segment after a restart.Source: D Proof on the equivalence of… · p. 56

E Projected-newton for weighted Tv1D

• Taut-string approaches are empirically superior to projected Newton, though the derivation may help develop subroutines for ℓp-TV prox-operators and related problems.The passage names group total-variation and ℓ1-trend filtering [87] as examples of related structures.Source: E Projected-newton for weighted Tv1D · p. 57
• Projected Newton solves the weighted-TV dual, a bound-constrained quadratic program, by exploiting its structure rather than using generic solvers off the shelf.The paper notes that TRON, L-BFGS-B, or projected Newton could solve the QP, but structure exploitation is needed for efficiency.Source: E Projected-newton for weighted Tv1D · p. 57
• At each iteration, projected Newton identifies active constraints, forms a reduced Hessian, scales the gradient, takes a damped step, and projects onto the constraints.Variables at a bound whose gradient points inside the feasible region are excluded from the update.Source: E Projected-newton for weighted Tv1D · pp. 57–58
• The reduced Hessian remains symmetric tridiagonal, enabling linear-time Cholesky decomposition and linear-time solution of the resulting systems.This structure holds for any active set.Source: E Projected-newton for weighted Tv1D · p. 58
• The stepsize uses backtracking with quadratic interpolation because the original Armijo search along the projection arc is inordinately expensive for this problem.The gradient is modified at boundary variables when the proposed direction points outside the feasible set.Source: E Projected-newton for weighted Tv1D · p. 58
• Each full projected-Newton iteration costs O(n), reducing the general method’s O(n^3) per-iteration cost by exploiting problem structure.When the unconstrained optimum is feasible, it is obtained directly in O(n); the implementation uses a duality-gap tolerance of 10^-5.Source: E Projected-newton for weighted Tv1D · p. 59

F Testing images and videos, and experimental results

• The image set includes public shape and phantom images, QR barcodes generated from random text, trollface and comic, a high-resolution Gaudi photograph, and author-created images.Gaudi’s Casa Batlló image is 3197 × 3361 pixels and is used in multicore experiments.Source: F Testing images and videos, and… · p. 59, Figure 25 · p. 60
• The video experiments use the publicly available salesman, coastguard, and bicycle sequences, with selected salesman frames displayed.The appendix captions show noisy and denoised frames from the salesman sequence.Source: F Testing images and videos, and… · p. 59, Figure 30 · p. 65, Figure 31 · p. 66
• The appendix displays test-image denoising results and deconvolution results, alongside noisy and convoluted image versions.The captions identify the figure contents but do not state comparative outcomes.Source: Figure 26 · p. 61, Figure 27 · p. 62, Figure 28 · p. 63