Iterative Thresholding with Inversion Algorithm

• ITI is a sparse signal recovery method that employs iterative hard-thresholding integrated with inversion and null-space tuning for enhanced performance.
• The algorithm leverages an eigenvalue-based suboptimal feedback approach to significantly reduce computational complexity while maintaining high accuracy.
• Empirical studies show that ITI achieves comparable NMSE to projection-based methods with dramatically lower runtime in large-scale compressive sensing settings.

Iterative Thresholding with Inversion (ITI) is a framework for sparse signal recovery from undersampled linear measurements. The method is fundamentally characterized by the incorporation of an inversion step—either exact or suboptimal—within an iterative hard-thresholding process, coupled with additional null-space correction. ITI, also referred to as Null-Space Tuning with Hard Thresholding and Feedback (NST+HT+FB), is notable for achieving recovery accuracy comparable to projection-based algorithms while offering significantly improved computational scalability, particularly when leveraging an eigenvalue-based suboptimal feedback strategy. This makes ITI applicable to large-scale compressive sensing and related inverse problems (Song et al., 2017).

1. Mathematical Formulation and Algorithmic Structure

The ITI algorithm addresses the canonical sparse recovery model:

$y = A x + e,$

where %%%%1%%%% is $k$-sparse, $A \in \mathbb{R}^{m \times n}$ ($m \ll n$) is the measurement matrix, $y \in \mathbb{R}^m$ are the measurements, and $e \in \mathbb{R}^m$ is noise. Sparsity is enforced by a hard thresholding operator $H_k$, which retains only the $k$ largest (in magnitude) entries of a vector.

Each ITI iteration proceeds as follows:

1. Feedback Step: Compute a proxy

$$g^{(t)} = x^{(t)} + A^\top (A A^\top)^{-1}(y - A x^{(t)}).$$

2. Hard Thresholding: Compute sparse proxy

$\mu^{(t)} = H_k(g^{(t)}).$

3. Feasibility Tuning: Project back to the feasible set using the orthogonal projector onto the null space of $A$,

$$x^{(t+1)} = \mu^{(t)} + A^\top (A A^\top)^{-1}(y - A \mu^{(t)}).$$

Alternatively, a merged update is available:

$$x^{(t+1)} = H_k\left(x^{(t)} + A^\top (A A^\top)^{-1}(y - A x^{(t)})\right).$$

These iterations ensure that recovered vectors gradually approach sparsity and measurement consistency as defined by $A x = y$. The null-space tuning operator is given by

$P = I - A^\top(AA^\top)^{-1} A,$

which maps arbitrary updates back to the feasible set (Song et al., 2017).

2. Suboptimal Feedback and Eigenvalue-Based Approximation

In standard ITI or NST+HT+FB, support-based inversion operations such as inversion of $A_T^\top A_T$ (where $T$ indexes the active support) are required, incurring $O(k^3)$ computational cost per iteration. For large-scale or high-sparsity settings, this cost is prohibitive.

To reduce complexity, the NST+HT+subOptFB variant replaces the explicit inversion with an eigenvalue-based scalar approximation. Specifically, under the restricted isometry property (RIP) and preconditioned RIP (P-RIP), the Gram matrix $A_T^\top A_T$ has eigenvalues in $[1-\delta, 1+\delta]$, with $\delta$ being an appropriate RIP constant. Its inverse is replaced by $\lambda^{-1} I$, where $\lambda \in [1 - \delta, 1 + \delta]$ is chosen to control contraction.

The update becomes:

• Identify $T_t$ as the $k$ largest-magnitude indices of $x^{(t)}$.
• For entries in $T_t$, set

$$\mu^{(t)}_{T_t} = x^{(t)}_{T_t} + \lambda^{-1} A_{T_t}^\top (A x^{(t)} - A_{T_t} x^{(t)}_{T_t}),$$

and $\mu^{(t)}_{T_t^c} = 0$.

• Update as before:

$$x^{(t+1)} = \mu^{(t)} + A^\top (A A^\top)^{-1}(y - A \mu^{(t)}).$$

This version preserves the key properties of ITI while reducing per-iteration complexity from $O(mn + k^3)$ to $O(mn + n)$ (Song et al., 2017).

3. Convergence Guarantees

The convergence analysis employs the preconditioned restricted isometry property (P-RIP). For every $s$-sparse vector $u$,

$$(1 - \gamma_s) \|u\|_2^2 \leq \| (A A^\top)^{-1/2} A u \|_2^2 \leq (1 + \gamma_s) \|u\|_2^2,$$

where $\gamma_s$ is the P-RIP constant.

Theoretical results show that for $y = A x + \eta$ with $x$ $s$-sparse and $\|\eta\|_2 \leq \epsilon$, if $\gamma_{3s} < \sqrt{2}/4$, and $\lambda$ is sufficiently large (dependent on $\gamma_{3s}$ and related constants), then the iterates $\mu^{(t)}$ satisfy

$$\|x - \mu^{(t)}\|_2 \leq \rho_{3s}^t \|x - \mu^{(0)}\|_2 + \kappa_{3s} \frac{1 - \rho_{3s}^t}{1 - \rho_{3s}} \cdot \epsilon,$$

with explicit forms for $\rho_{3s}<1$ and $\kappa_{3s}$ provided in (Song et al., 2017).

In the special case where $A$ is a Parseval frame ($A A^\top = I$), one may guarantee $\rho<1$ by choosing $\lambda > 3.9$ whenever $\delta_{3s} \leq 1/4$.

4. Computational Complexity

The computational cost per iteration is summarized as follows:

Algorithm Variant	Per-Iteration Complexity	Dominant Costs
NST+HT+FB (exact feedback)	$O(mn + k^3)$	Gram matrix inversion on active support ($O(k^3)$)
NST+HT+subOptFB (suboptimal)	$O(mn + n)$	Matrix-vector multiply ($O(mn)$); thresholding ($O(n)$)
HTP (Hard-Threshold Pursuit)	$O(mn + k^3)$	Inversion of $k \times k$ matrix

The eigenvalue-based suboptimal feedback approach eliminates the per-iteration $O(k^3)$ bottleneck, making the method scalable to problem sizes with $n \geq 10^5$ and large $k$. The complexity is ultimately controlled by the cost of matrix-vector multiplication (Song et al., 2017).

5. Adaptive Iterative Thresholding Without Sparsity Knowledge

The adaptive algorithm (AdptNST+HT+subOptFB) addresses scenarios where the true sparsity level $k$ is unknown. The procedure incrementally increases the support size at each iteration. At iteration $t$:

1. Form proxy $$g^{(t)} = x^{(t)} + \lambda^{-1} A^\top(y - A x^{(t)})$$.
2. Let $T_t$ be the indices of the $t$ largest magnitudes of $g^{(t)}$.
3. Set $\mu^{(t)}_{T_t} = g^{(t)}_{T_t}$, others zero.
4. Null-space tune: $$x^{(t+1)} = \mu^{(t)} + A^\top (A A^\top)^{-1}(y - A \mu^{(t)})$$.

Convergence holds under analogous P-RIP conditions, with geometric error decay after support recovery in $O(s)$ iterations, where $s$ is the true sparsity level (Song et al., 2017).

6. Empirical Performance and Comparative Analysis

Extensive simulations reported in (Song et al., 2017) include problem sizes up to $n=10^5$, measurement ratios $m/n$ up to $0.8$, and sparsity ratios up to $0.4$. The ITI variants are benchmarked against:

• NST+HT+FB (exact inversion)
• NST+HT+subOptFB (suboptimal, eigenvalue-based feedback)
• HTP (Hard-Thresholding Pursuit)
• GHTP (graded HTP)
• Adaptive versions: AdptNST+HT+FB, AdptNST+HT+subOptFB

Key findings include:

• The NMSE (normalized mean-square error) achieved by NST+HT+subOptFB is on par with exact feedback and HTP across a broad range of undersampling/sparsity regimes.
• The runtime of NST+HT+subOptFB grows linearly in $n$ and $m$ and is one order of magnitude faster than HTP and exact-feedback NST+HT+FB at $n=10^5$. At $n=10^5$, $m/n=0.5$, $s/m=0.3$, NST+HT+subOptFB is 3–5× faster than exact-feedback NST+HT+FB.
• The adaptive variant AdptNST+HT+subOptFB achieves comparable NMSE to its competitors but with the lowest per-iteration cost and overall runtime.

7. Significance and Applicability

The ITI framework with null-space tuning and eigenvalue-based suboptimal feedback is distinguished by its combination of rapid geometric convergence and low per-iteration complexity. Particularly under mild P-RIP conditions, this makes ITI a compelling choice for sparse recovery in very large-scale systems, such as large-dimensional compressive sensing, where conventional inversion-based schemes become infeasible (Song et al., 2017).