Source-linked AI summary
Implementation of an Optimal First-Order Method for Strongly Convex Total Variation Regularization
Tobias Lindstrøm Jensen, Jakob Heide Jørgensen, Per Christian Hansen, Søren Holdt Jensen
TL;DR
Large tomography problems restrict algorithm choices, while Nesterov’s method requires parameters that are often unavailable in practice. The paper estimates those parameters during optimization and tests the resulting method on 3D TV reconstruction. UPN is competitive with other first-order methods and superior on ill-conditioned problems solved to high accuracy.
Problem
General linear operators and large 3D tomography systems limit algorithm choices, while Nesterov’s method requires strong-convexity and Lipschitz parameters that are often unavailable in practice.
Source: §1 · p. 2Method
UPN estimates the Lipschitz constant by backtracking and adjusts the strong-convexity estimate using a local inequality, with restarts when an estimate is inadequate.
Source: §1 · p. 2Results
UPN is competitive with other first-order algorithms and superior for ill-conditioned problems, including high-accuracy solutions.
Source: §1 · p. 2, §9 · p. 17Takeaways & Limitations
Simulations report favorable convergence for UPN even on problems that are not strongly convex.
Source: §9 · p. 17UPN is theoretically suboptimal relative to the stated optimal complexity of O(√Q log ε^-1).
Source: §6 · p. 12
Abstract
from arXiv · showhide
We present a practical implementation of an optimal first-order method, due to Nesterov, for large-scale total variation regularization in tomographic reconstruction, image deblurring, etc. The algorithm applies to $\mu$-strongly convex objective functions with $L$-Lipschitz continuous gradient. In the framework of Nesterov both $\mu$ and $L$ are assumed known -- an assumption that is seldom satisfied in practice. We propose to incorporate mechanisms to estimate locally sufficient $\mu$ and $L$ during the iterations. The mechanisms also allow for the application to non-strongly convex functions. We discuss the iteration complexity of several first-order methods, including the proposed algorithm, and we use a 3D tomography problem to compare the performance of these methods. The results show that for ill-conditioned problems solved to high accuracy, the proposed method significantly outperforms state-of-the-art first-order methods, as also suggested by theoretical results.
1 Introduction
The paper adapts Nesterov’s method for total-variation reconstruction by estimating its otherwise unknown parameters, then evaluates it on tomography problems. It also formulates a smoothed, bound-constrained discrete TV problem for 3D reconstructions.
1 Introduction
- Large 3D tomography systems with general linear operators limit algorithm choice, and their memory requirements preclude second-order methods with quadratic convergence.TV regularization is used when reconstructions need sharp edges.Source: §1 · p. 2
- Nesterov’s method requires known strong-convexity and Lipschitz parameters, which are often unavailable in practice.The paper estimates the Lipschitz constant by backtracking and adjusts the strong-convexity estimate using a local inequality, with restarts to ensure convergence when an estimate is inadequate.Source: §1 · p. 2
- The paper compares UPN with gradient-projection variants GP and GPBB, and with UPN0, which does not enforce strong-convexity information.GPBB uses Barzilai–Borwein acceleration; UPN0 is described as similar to FISTA [4].Source: §1 · p. 2
- UPN is significantly faster than GP and, for ill-conditioned problems, significantly faster than GPBB; it is consistently faster than UPN0 at high accuracy.For moderately ill-conditioned problems, UPN is as fast as GPBB.Source: §1 · p. 2
2 The Discrete Total Variation Reconstruction Problem
- The reconstruction represents the image as a 3D voxel array and approximates its gradient with forward differences, using periodic boundary conditions.The paper smooths the TV functional using the Huber function, described as prox-function smoothing.Source: §2 · p. 3
- The discrete 3D reconstruction assumes voxel values in [0,1], imposing the bound-constrained feasible set Q = {x | 0 ≤ x_j ≤ 1}.The problem uses a linear system Ax ≈ b and TV regularization.Source: §2 · p. 4
3 Smooth and Strongly Convex Functions
This section defines smoothness and strong convexity and relates them to the condition number governing convergence. For the smoothed TV reconstruction objective, it gives parameter bounds and notes that they need not be tight.
3 Smooth and Strongly Convex Functions
- The section frames the analysis around minimizing a convex, continuously differentiable objective over a convex feasible set.Its smoothness condition is also described as Lipschitz continuity of the gradient.Source: §3 · p. 4
- The condition number Q = L/µ governs convergence rates and, when A has full column rank, equals κ(A)^2.The paper defines the function class F_µ,L using strong-convexity parameter µ and Lipschitz constant L.Source: §3 · pp. 4–5
- For the reconstruction objective, the stated strong-convexity parameter is λ_min(A^T A), and the Lipschitz constant is ||A||_2^2 + α||D||_2^2/τ.If A has full column rank, µ = σ_min(A)^2 > 0; if rank(A) < N, the stated µ is zero.Source: §3 · p. 5
- The derived µ and L are not guaranteed to be the tightest possible parameters for the objective.For rank-deficient A, some problems have a singular Hessian and µ = 0, though the paper does not claim this always occurs.Source: §3 · p. 5
4 Some Basic First-Order Methods
The section contrasts gradient projection and Barzilai–Borwein approaches with optimal first-order methods. For strongly convex objectives, the latter have a better dependence on the condition number, but Nesterov’s method requires known parameters.
4 Some Basic First-Order Methods
- The Barzilai–Borwein step cheaply approximates a Newton step, but precise iteration complexity is difficult to give for the resulting general algorithms.The section says not to expect better than gradient-method complexity from BB for the stated function classes.Source: §4 · p. 6
- The paper defines optimal complexity as worst-case complexity over a problem class and reviews methods for both F_0,L and F_µ,L.It also notes O(L/ε) complexity for optimal methods on the non-strongly-convex class.Source: §4 · p. 6
- Nesterov’s optimal method for strongly convex objectives requires known µ and L, with an initial parameter θ_0 ≥ √(µ/L).The method maintains two vector sequences, x^(k) and y^(k).Source: §4 · p. 7
- For ill-conditioned problems, optimal methods’ √Q dependence offers an advantage over the gradient method’s Q dependence.The section motivates specialized optimal first-order methods on this basis.Source: §4 · p. 8
5 First-Order Inequalities for the Gradient Map
The section develops inequalities for the gradient map to support parameter estimation and stopping decisions. It relates the map to optimality and to objective error for constrained problems.
5 First-Order Inequalities for the Gradient Map
- The gradient map generalizes the gradient to constrained optimization and vanishes at a minimizer as a necessary and sufficient optimality condition.When Q = R^N, the map equals ∇f(x).Source: §5 · p. 8
- Local inequalities can permit estimates with µ̄ ≥ µ and L̄ ≤ L, which the proposed algorithm uses when the parameters are unknown.The inequalities hold for particular iterates even though µ and L remain valid general choices.Source: §5 · pp. 8–9
- The gradient-map norm is used in a stopping criterion because the section relates it to distance from the minimizer and objective error.The criterion uses the current iterate and a user-specified tolerance based on requested accuracy.Source: §5 · p. 9
6 Nesterov’s Method With Parameter Estimation
UPN adapts Nesterov’s method by estimating the unknown Lipschitz and strong-convexity parameters during iteration. Backtracking, a local strong-convexity heuristic, and restarts support its use when those parameters are unavailable.
6 Nesterov’s Method With Parameter Estimation
- UPN estimates Lk by backtracking and µk with a heuristic because computing the required parameters explicitly is too expensive.The µ estimate is chosen from known iterates and is intended to decrease toward a better estimate.Source: §6 · pp. 9–10
- The µk heuristic selects the largest estimate satisfying the local strong-convexity inequality, without exceeding the previous estimate.Because an estimate may be too large, the restart procedure detects inadequacy and restarts with a reduced estimate from the current iterate.Source: §6 · pp. 10–11
- UPN adds an initial projected-gradient step for a restart bound, a stopping-criterion step, parameter-estimation logic, and a restart for inadequate µ estimates.The complete method combines Nesterov’s method with backtracking on Lk and the µk heuristic.Source: §6 · pp. 11–12
- In the worst case, one UPN iteration may require twice the computational work of a GPBB iteration because it makes two backtracking calls.The authors note the calls may be implemented more efficiently than simply running backtracking twice.Source: §6 · p. 12
- Under the stated parameter conditions, UPN stops using a gradient-map criterion and returns an ε-suboptimal solution, but its iteration complexity is suboptimal to the optimal bound.The method can nevertheless be applied when µ and L are unknown.Source: §6 · p. 12
7 The 3D Tomography Test Problem
The test problem models 3D tomography as a sparse linear system built from ray paths through voxels. It uses a 3D Shepp–Logan phantom with orthogonal projections distributed over the sphere.
7 The 3D Tomography Test Problem
- Tomographic reconstruction estimates an object from projections produced by rays such as X-rays or seismic waves, whose intensities are partially absorbed by the object.The absorption provides information about the object.Source: §7 · pp. 12–13
- Figure 7.1 shows orthogonal phantom slices, a central horizontal slice, and a solution for α = 1 and τ = 10^-4; smaller α yields a less smooth solution.Display colors are clipped to [0.1, 0.4] for contrast, while original voxel values are 0.0, 0.2, 0.3, and 1.0.Source: Figure 7.1 · p. 13
- The generic tomography model reconstructs an approximation of the object’s attenuation function from recorded ray measurements.Ray intensity decays according to a line integral through the object.Source: §7 · p. 13
- Each ray’s line integral is discretized as a sum of its path-length contributions through voxels, yielding a sparse system Ax = b.The matrix has M rows for rays or measurements and N columns for voxels.Source: §7 · p. 13
- The phantom is a 43×43×43-voxel 3D Shepp–Logan object, reconstructed using parallel-beam projections in directions distributed over the unit sphere.The projection directions are Lebedev quadrature points.Source: §7 · pp. 13–14
8 Numerical Experiments
Experiments compare UPN with gradient-projection and FISTA-like variants on TV-regularized tomography problems. UPN is especially effective at high accuracy on ill-conditioned problems, and it also converges linearly in the tested non-strongly convex case.
8 Numerical Experiments
- The experiments use two representative test problems with noisy data at 1% relative noise, and initialize TV algorithms from the fifth conjugate-gradient iteration for least squares.Accuracy comparisons use a reference solution computed at accuracy ε̄·10^-4.Source: §8 · p. 14
8.1 Influence of α and τ on the convergence
- For larger condition numbers, UPN converges significantly faster at high accuracy; at low accuracy, all four methods perform more or less equally well.For low α/τ ratios, GPBB and GP can require comparable or fewer iterations; when UPN is not best, it requires only slightly more iterations than the best alternative.Source: §8.1 · p. 15
8.2 Restarts and µk and Lk histories
- For T1 with α = 100 and τ = 10^-4, restarts occur in the first eight iterations, reducing µk by a factor of 0.7 each time; afterward, µk and Lk remain constant.At α = 1 and τ = 10^-4, no restarts occur; µk declines gradually except for a final jump, while Lk remains almost constant.Source: §8.2 · pp. 15–16
8.3 A non-strongly convex example
- On non-strongly convex T2, UPN still converges linearly, slightly slower than on T1; its linear convergence makes it faster than the alternatives for high-accuracy solutions.At low accuracy UPN0 is comparable to UPN, while GP and GPBB converge much more slowly.Source: §8.3 · p. 16
8.4 Influence of the heuristic
- On T1, UPN with heuristic parameters converges faster than UPN using the true µ and L, which the authors attribute to the large strong-convexity modulus implied by the true parameters.The heuristic chooses locally sufficient µk values rather than restricting each step to a globally valid µ.Source: §8.4 · pp. 16–17
9 Conclusion
The paper implements Nesterov’s optimal first-order method for large-scale optimization by estimating smoothness and strong-convexity parameters during iterations. In 3D TV-regularized tomography, it is competitive overall and superior on ill-conditioned problems requiring high accuracy.
9 Conclusion
- The implementation estimates smoothness and strong-convexity parameters during iterations, making Nesterov’s method practical when those parameters are unknown.The method is intended for large-scale problems with smooth, strongly convex objectives.Source: §9 · p. 17
- In 3D TV-regularized tomography, UPN is competitive in iteration count and superior on ill-conditioned problems solved to high accuracy.The experiments also report favorable strong-convexity-associated convergence rates in practice for problems that are not strongly convex.Source: §9 · p. 17
A The Optimal Convergence Rate
The appendix analyzes an optimal first-order method using iteration-dependent estimates of smoothness and strong convexity. It concludes that the method attains the optimal convergence rate for both cited function classes.
A The Optimal Convergence Rate
- The analysis introduces iteration-dependent Lk and µk estimates to support reasoning about how those parameters should be selected.It uses relations involving the supporting variables z(k) and γk.Source: A The Optimal Convergence Rate · p. 18
- The appendix examines the method’s recurrences and inequalities, including an update relation for γk+1z(k+1).The relation combines the previous γkz(k), µky(k), and gradient-map term GLk(y(k)).Source: A The Optimal Convergence Rate · p. 18
- The resulting convergence rate is stated to be optimal simultaneously for the classes F0,L and Fµ,L.The paper identifies this as the optimal rate for both classes.Source: A The Optimal Convergence Rate · p. 19
B Complexity Analysis
The appendix derives UPN’s worst-case complexity from the number of restarts and the cost of successful termination. In the worst case, the heuristic never reduces the strong-convexity estimate within a restart stage.
B Complexity Analysis
- UPN’s total worst-case complexity combines the worst-case number of restarts with the worst-case cost of successful termination.The appendix identifies these as the two components of the complexity bound.Source: B Complexity Analysis · p. 19
- In the worst case, the heuristic never reduces μ_k, so μ_k,r = μ_0,r throughout restart stage r.Source: B Complexity Analysis · p. 19
- The restart count must satisfy R ≥ log(μ_0,0/μ)/log(1/ρ_μ) for the initial estimate to reach the target bound.Source: B Complexity Analysis · p. 19
B.1 Termination Complexity
This section bounds UPN’s termination and restart costs, then combines them into a total iteration complexity. It also relates the stopping tolerance to ε-suboptimality.
B.1 Termination Complexity
- After at most R restarts, a sufficiently small μ_0,r makes condition (6.9) hold, allowing termination complexity to be bounded.Source: §B.1 · p. 20
B.2 Restart Complexity
- Restart detection uses rule (6.9); before restarting, the gradient-mapping norm remains at least ε̄, and the analysis bounds the iterations across R restarts.Source: §B.2 · pp. 20–21
B.3 Total Complexity
- UPN’s total iteration complexity is the sum of the restart and successful-termination bounds.Source: §B.3 · p. 21
- Setting (2/μ)ε̄² = ε relates the stopping criterion to ε-suboptimality and yields the iteration complexity in terms of ε.Source: §B.3 · p. 21