PaperlayerView exact revision on arXiv

Source-linked AI summary

Sparse Approximation via Polynomial Equations

Matija Tomić, Raphaël Widdershoven, Lieven De Lathauwer

arXiv:2609.11215v1eess.SPcs.NAmath.NA

TL;DR

Sparse approximation asks for sparse solutions to underdetermined linear systems, a problem for which brute-force search is intractable. This paper recasts sparsity as a polynomial system and develops tensor-based recovery methods; experiments report stronger recovery than BP and OMP, alongside higher computational cost.

  • Problem

    Sparse approximation seeks sparse solutions to underdetermined linear systems, while brute-force search is intractable and practical methods offer approximate solutions with guarantees under specific conditions.

    Source: §I · p. 1

  • Method

    The paper formulates sparse approximation as polynomial equations, uses a Macaulay construction to recover all roots, and proposes two optimization-based alternatives.

    Source: §I · p. 1, C. The Macaulay Construction · p. 4

  • Results

    Across experiments, SESP-D and SESP-P recovered signals from fewer measurements than BP and OMP and succeeded at higher sparsity, while with m = 32 measurements they recovered a 12-sparse signal in approximately 90% of trials versus roughly 50% for BP and OMP.

    Source: E. Recovery and Runtime in the… · p. 7, §VII · p. 8

  • Takeaways & Limitations

    The methods recover multiple equivalent sparse representations and let users specify target sparsity directly as an upper bound rather than an exact requirement.

    Source: §VII · p. 8

    The proposed methods have higher computational cost, most notably in the noisy setting; BP or BPDN remains cheaper when it succeeds.

    Source: §VII · p. 8

Abstract

from arXiv · show

We consider the problem of finding sparse solutions of an underdetermined linear system $Ax=b$. In contrast to conventional approaches based on greedy algorithms or convex relaxation, we reformulate sparse approximation as a structured system of polynomial equations and connect with the literature on tensor methods. We develop an eigenvalue decomposition based method that formally guarantees recovery of all sparse solutions if there is more than one. We also develop two optimization-based methods achieving favorable computational complexity. The new methods allow explicit control of the target sparsity. Numerical experiments illustrate the performance and compare to basis pursuit (denoising) and orthogonal matching pursuit.

I. INTRODUCTION

Sparse approximation seeks sparse solutions to underdetermined linear systems, but brute-force search is intractable at scale. The paper formulates sparsity as polynomial constraints and connects the resulting system to tensor methods.

I. INTRODUCTION

  • Brute-force search is computationally intractable for large n and sparse approximation is NP-hard; greedy and convex-relaxation methods provide approximate solutions with guarantees under specific conditions.Examples include Orthogonal Matching Pursuit (OMP) – and Basis Pursuit (BP) and Basis Pursuit Denoising (BPDN).Source: §I · p. 1

II. SPARSE APPROXIMATION AND POLYNOMIALS

  • Products of every s + 1 distinct entries of an s-sparse vector vanish, so sparse approximation can be written as Ax = b plus sparsity-inducing monomial equations.Solutions sparser than the specified level s are also included.Source: §II · pp. 1–2
  • The polynomial formulation represents all sparse solutions, including cases where the sparsest solution is not unique.Geometrically, the sparsity constraints define a union of hyperplanes, whose intersections with Ax = b give the sparse solutions.Source: §II · pp. 1–2

A. Explicit Canonical Polyadic Decomposition

  • A tensor can be represented by a Canonical Polyadic Decomposition (CPD), and the paper uses nonlinear least squares to compute a target-rank CPD.The CPD is unique up to permutation and scaling under mild conditions.Source: A. Explicit Canonical Polyadic… · p. 2

B. Implicit CPD

The paper uses constrained CPD to handle tensors specified implicitly, then encodes sparsity and data fit in a structured multilinear system. For multiple polynomial roots, coefficient matrices place root-specific monomial vectors in a shared null space.

B. Implicit CPD

  • When a linear system defining a tensor is underdetermined, a CPD-constrained solution can be found by minimizing a nonlinear least-squares objective.The implicit CPD computation requires the objective, its gradient, and the Jacobian Gramian.Source: B. Implicit CPD · p. 2

C. Monomials and Rank-1 Tensors

  • The sparsity-promoting monomials form a rank-1 tensor system, and adding Ax = b extends it to enforce both sparsity and data fit.A selection matrix picks the distinct-index terms, which the system constrains to zero.Source: C. Monomials and Rank-1 Tensors · p. 2

IV. POLYNOMIAL EQUATIONS AND MULTIPLE SOLUTIONS

  • Unlike a formulation for one solution, the null space of a coefficient matrix jointly describes all solutions of a polynomial system.The paper frames recovery of multiple roots as separating rank-1 terms within that null space.Source: §IV · p. 3, A. Roots as Structured Null Space… · p. 3

A. Roots as Structured Null Space Vectors

  • For a system with Q distinct roots, each root yields a Vandermonde vector in the null space of M0, whose columns can be collected into a shared matrix.The vectors contain distinct monomials evaluated at each root; the setup assumes Q distinct roots.Source: A. Roots as Structured Null Space… · p. 3

B. Recovering the Roots from the Null Space

A numerical null-space basis can be reshaped into a tensor whose CPD encodes the roots. Under the stated conditions, an eigenvalue decomposition uniquely recovers the factors.

B. Recovering the Roots from the Null Space

  • A numerical null-space basis can be reshaped into a tensor whose CPD rank-1 terms encode the roots.The second factor matrix contains the roots.Source: B. Recovering the Roots from the Null… · p. 3
  • Under condition (10), the decomposition is unique, and its factors can be computed by a simple matrix eigenvalue decomposition.The condition ensures full column rank for two factors and pairwise non-collinear columns in the third because the roots are distinct.Source: B. Recovering the Roots from the Null… · p. 3

C. The Macaulay Construction

The Macaulay construction represents polynomial-system roots through matrix null spaces, while the paper also develops optimization-based sparse-approximation methods. The construction can recover all sparse solutions, but its matrix size makes computation costly.

C. The Macaulay Construction

  • When assumption (10) fails for a general polynomial system, extra null-space directions can invalidate the basis condition in (11).Those directions do not correspond to roots.Source: C. The Macaulay Construction · p. 3
  • Adding rows from shifted polynomials preserves every root while removing null-space directions that do not correspond to roots.The Macaulay matrix includes the original coefficient rows and coefficients of the shifted polynomials.Source: C. The Macaulay Construction · p. 3
  • At the degree of regularity, the Macaulay null space is spanned only by vectors corresponding to roots, and its stabilized dimension gives the root count.The degree is increased until the null-space dimension stabilizes; a variant of the root-recovery method is then applied.Source: C. The Macaulay Construction · pp. 3–4
  • Under the stated column-independence condition, sparsity-inducing equations turn the underdetermined linear system into a finite set with at most one solution per support of size at most s.The Macaulay procedure can recover all such roots at once, including solutions sparser than s.Source: C. The Macaulay Construction · p. 4
  • The Macaulay matrix can be costly to compute, and dropping columns for monomials constrained to zero provides a smaller matrix for the null-space calculation.The passage gives the matrix size as Nd but does not provide the completed expression.Source: C. The Macaulay Construction · p. 4

V. EFFICIENT OPTIMIZATION-BASED ALGORITHMS

  • Exploiting the selection-matrix structure of Z expresses Gauss–Newton quantities through elementary symmetric polynomials, which can be evaluated cheaply.The stated structure supports computationally efficient implementations.Source: §V · p. 4

A. Algorithms

  • SESP-D jointly optimizes data fitting and sparsity promotion, whereas SESP-P seeks the sparsest exact fit by optimizing over the affine solution set x = xp + Vy.Each returns one solution per initialization; reinitializing may yield different solutions when multiple sparse solutions exist.Source: A. Algorithms · p. 4

B. Computational Complexity of the LS-CPD-Based Methods

The paper reduces the computational cost of the optimization-based methods by exploiting recursions for elementary symmetric polynomials. Its complexity accounting separates sparsity-enforcement costs from the costs of enforcing Ax = b.

B. Computational Complexity of the LS-CPD-Based Methods

  • The recursion reduces evaluation of the SESP-D and SESP-P objective functions from naive O(ns+1) to O(ns) operations.The stated elementary-symmetric-polynomial computation is O(nd) for degree d.Source: B. Computational Complexity of the… · p. 4
  • The recursion reduces the costs of J⊤r and J⊤J from O(n2s) and O(n3s) to O(ns) and O(n2s), respectively.These are per-iteration quantities; iteration-independent terms such as V and A⊤A are precomputed once.Source: B. Computational Complexity of the… · pp. 4–5
  • Since s ≪ m, the sparsity-enforcing terms in Table I do not dominate the overall per-iteration complexity.The terms ns and n2s enforce sparsity; nm, n(n −m), and n2(n −m) enforce Ax = b.Source: B. Computational Complexity of the… · p. 5, Table I · p. 5

VI. EXPERIMENTS

  • The experiments compare the proposed methods with BP/BPDN and OMP, evaluating recovery accuracy and average runtime in noiseless and noisy settings.The algebraic framework also allows explicit target-sparsity control and recovery of multiple valid sparse representations.Source: §VI · p. 5

1) SESP-D and SESP-P:

The experimental section specifies implementation and evaluation procedures for SESP-D, SESP-P, BP/BPDN, and OMP. It also describes initialization and reports that restarts may help SESP-D at higher sparsity.

1) SESP-D and SESP-P:

  • The proposed methods use a Gauss–Newton-NLS dogleg trust-region implementation, with a maximum of 2000 iterations and support-restricted least-squares stopping.The stopping threshold is 10−10 in noiseless trials and σ√m in noisy trials.Source: 1) SESP-D and SESP-P · p. 5

2) BP/BPDN:

  • Noiseless BP uses MATLAB’s linprog interior-point solver, while noisy BPDN uses SPGL1.The noiseless BP solver uses optimality tolerance 10−8 and constraint tolerance 10−6, with a maximum of 2000 iterations.Source: 2) BP/BPDN · p. 5

3) OMP:

  • OMP greedily selects the column most correlated with the current residual and stops when the residual norm falls below its prescribed threshold.The thresholds are 10−10 in noiseless trials and σ√m in noisy trials.Source: 3) OMP · p. 5

B. General Experimental Setup

  • Trials use normalized Gaussian sensing-matrix columns and randomly supported Gaussian signals; sparsity is known, and noisy-case performance is reported as relative error.Noiseless success requires relative recovery error below 10−5, with success rate measured over 100 trials.Source: B. General Experimental Setup · p. 5

C. Initialization

  • Because the methods are nonconvex and Gauss–Newton is local, initialization may affect the outcome.For SESP-D, the initialization scale is chosen to equalize the expected initial data-fitting and sparsity-promotion terms.Source: C. Initialization · p. 6, 1) SESP-D · p. 6

1) SESP-D:

  • At low sparsity, one SESP-D initialization suffices; at higher sparsity, reinitializations may improve success rate.Figure 2 plots SESP-D success rate against sparsity for different numbers of random initializations.Source: Figure 2 · p. 6

2) SESP-P:

The experiments compare SESP-D and SESP-P with BP/BPDN and OMP across initialization strategies, regularization choices, and noiseless and noisy settings. Recovery gains are accompanied by a runtime trade-off, especially with noise.

2) SESP-P:

  • For SESP-P, minimizing ∥xp + Vy0∥2 gives the initialization y0 = −V⊤xp, using the orthonormal null-space basis V.A BP/BPDN solution can instead warm-start SESP-P with y0 = V⊤(xBP/BPDN − xp).Source: 2) SESP-P · p. 6, 3) Warm Start from BP/BPDN · p. 6

3) Warm Start from BP/BPDN:

  • Warm-starting SESP-D and SESP-P from BP further improves recovery beyond BP alone, and SESP-P outperforms BP as sparsity increases.The polynomial-based refinement further sparsifies BP's ℓ1-optimal solution.Source: 3) Warm Start from BP/BPDN · p. 6

D. Regularization Parameter

  • For low-to-medium sparsity, noiseless success is insensitive to λ across several orders of magnitude; at higher sparsity, quite small λ values tend to be preferable.Smaller λ places relatively more weight on data fitting than on sparsity promotion; noisy results show a similar effect.Source: D. Regularization Parameter · p. 6, Figure 4 · p. 6

E. Recovery and Runtime in the Noiseless Case

F. Recovery and Runtime in the Noisy Case

G. Sparse Modeling with Equivalent Representations

The equivalent-representations experiment tests recovery of multiple sparse solutions while controlling target sparsity; the conclusion reports favorable recovery comparisons alongside higher computational cost. The paper also contrasts the methods’ scope and practical trade-offs, including Macaulay’s small-problem limitation and open directions.

G. Sparse Modeling with Equivalent Representations

  • SESP-D and SESP-P recover all four equivalent sparse representations, whereas BP and OMP are limited to a single solution.The example evaluates multiple equivalent representations of the same signal.Source: G. Sparse Modeling with Equivalent… · p. 8
  • With target sparsity s=2, SESP-D and SESP-P always recover x̂1; with s=3, they may recover x̂2, x̂3, or x̂4, and may also recover x̂1 depending on initialization.Across 1000 random initializations per γ, the lower-sparsity x̂1 is most often recovered when target sparsity overestimates the minimum.Source: G. Sparse Modeling with Equivalent… · p. 8, Figure 10 · p. 8
  • Macaulay formally recovers all solutions at once in this example: x̂2, x̂3, and x̂4 have multiplicity one, while x̂1 has multiplicity three.The limited problem size makes this evaluation feasible, unlike the prior experiments.Source: G. Sparse Modeling with Equivalent… · p. 8

VII. CONCLUSION AND FUTURE WORK

  • Macaulay is tractable only for small problems, while SESP-D and SESP-P exploit sparsity-constraint structure through inexpensive ESP estimation and have per-iteration costs comparable to BP.The conclusion identifies these as optimization-based alternatives.Source: §VII · p. 8
  • The methods recovered sparse signals from fewer measurements than BP and OMP, succeeded at higher sparsity levels, and remained close to the oracle bound under additive noise, at higher computational cost.The cost increase is most pronounced in the noisy setting.Source: §VII · p. 8
  • Unlike BP and BPDN, the formulation specifies target sparsity directly as an upper bound, and can recover several equivalent sparse representations that BP and OMP do not provide.Source: §VII · p. 8
  • Open directions include error bounds for noisy estimation, polynomially expressible constraints, polynomially parameterized dictionaries, and faster Macaulay solvers for larger problems.Source: §VII · p. 8

APPENDIX A STRUCTURE OF Z⊤Z

Appendix A characterizes Z⊤Z as a projector onto basis vectors with pairwise-distinct indices. Appendix B uses the resulting structure to derive the SESP-D objective and Jacobian Gramian.

APPENDIX A STRUCTURE OF Z⊤Z

  • Because Z has one-hot rows, Z⊤Z is diagonal with entries in {0, 1} and projects onto basis vectors whose indices are pairwise distinct.The basis is indexed by multi-indices and formed from Kronecker products of standard basis vectors.Source: Appendix A · pp. 8–9

APPENDIX B SESP-D DERIVATION

  • The SESP-D objective combines the squared linear residual ∥Ax−b∥₂² with an E_{s+1} penalty weighted by λ²(s+1)!/2.Source: Appendix B · p. 9
  • The SESP-D Jacobian Gramian splits into a linear-system term A⊤A and a sparsity-residual Jacobian term.J(1)=A, while J(2) is the Jacobian of r(2).Source: Appendix B · p. 9
  • Each subset of s+1 distinct indices contributes (s+1)! equal entries to r(2), which share the same gradient.Source: Appendix B · p. 9

APPENDIX C SESP-P DERIVATION

The SESP-P derivation uses the same objective-building approach as Appendix B while optimizing in y coordinates, where x is represented as xp+Vy.

APPENDIX C SESP-P DERIVATION

  • The SESP-P objective is obtained using the same derivation approach as Appendix B.Source: Appendix C · p. 9
  • SESP-P optimizes over y rather than x, using the coordinate relation x=xp+Vy.Source: Appendix C · p. 9
  • The Jacobian of x with respect to y is determined through the map x=xp+Vy.Source: Appendix C · p. 9

APPENDIX D PROOF OF LEMMA 1

The appendices give a recursion for elementary symmetric polynomials and use left-right partial computations to obtain the stated costs for the subsequent quantities.

APPENDIX D PROOF OF LEMMA 1

  • The recursion e_j(z_1,…,z_k)=e_j(z_1,…,z_(k−1))+z_k e_(j−1)(z_1,…,z_(k−1)), initialized by e_0=1, computes e_d in O(nd) operations.Source: Appendix D · p. 10

APPENDIX E PROOF OF LEMMA 2

  • Left-to-right and right-to-left partial ESP tables each cost O(ns); partitioning subsets around each index then computes the needed values in O(ns) total.Source: Appendix E · p. 10

APPENDIX F PROOF OF LEMMA 3

  • Applying the partial-ESP computation to each index yields all E(i,j) values in O(n^2s) time.Source: Appendix F · p. 10
Loading 2609.11215v1…