Source-linked AI summary
Exact and Stable Covariance Estimation from Quadratic Sampling via Convex Programming
Yuxin Chen, Yuejie Chi, Andrea Goldsmith
TL;DR
The paper studies how to estimate structured covariance matrices from high-dimensional data when acquisition devices have limited memory and computation. It uses quadratic measurements with structure-specific convex recovery and establishes universal, stable, and near-minimal recovery guarantees across several covariance models. The framework also supports streaming and phaseless measurement applications while remaining robust to noise and imperfect structural assumptions.
Problem
High-dimensional covariance estimation must operate with a single pass, minimal storage, and low computational complexity, motivating recovery from a small number of measurements under structural assumptions.
Method
The paper combines quadratic rank-one measurements with convex relaxations tailored to low-rank, Toeplitz low-rank, sparse, and jointly sparse rank-one covariance structures.
Results
The proposed algorithms achieve exact universal recovery, stable recovery under noise and imperfect structure, and near-minimal measurement guarantees for most considered covariance structures.
Takeaways & Limitations
Quadratic covariance sketches provide a practical framework for streaming and phaseless applications, with low storage and computational costs and no additional incoherence conditions.
Abstract
from arXiv · showhide
Statistical inference and information processing of high-dimensional data often require efficient and accurate estimation of their second-order statistics. With rapidly changing data, limited processing power and storage at the acquisition devices, it is desirable to extract the covariance structure from a single pass over the data and a small number of stored measurements. In this paper, we explore a quadratic (or rank-one) measurement model which imposes minimal memory requirements and low computational complexity during the sampling process, and is shown to be optimal in preserving various low-dimensional covariance structures. Specifically, four popular structural assumptions of covariance matrices, namely low rank, Toeplitz low rank, sparsity, jointly rank-one and sparse structure, are investigated, while recovery is achieved via convex relaxation paradigms for the respective structure. The proposed quadratic sampling framework has a variety of potential applications including streaming data processing, high-frequency wireless communication, phase space tomography and phase retrieval in optics, and non-coherent subspace detection. Our method admits universally accurate covariance estimation in the absence of noise, as soon as the number of measurements exceeds the information theoretic limits. We also demonstrate the robustness of this approach against noise and imperfect structural assumptions. Our analysis is established upon a novel notion called the mixed-norm restricted isometry property (RIP-$\ell_{2}/\ell_{1}$), as well as the conventional RIP-$\ell_{2}/\ell_{2}$ for near-isotropic and bounded measurements. In addition, our results improve upon the best-known phase retrieval (including both dense and sparse signals) guarantees using PhaseLift with a significantly simpler approach.
I. INTRODUCTION
The paper addresses covariance estimation for high-dimensional, rapidly changing data under limited memory and computation by using quadratic measurements and structure-specific convex recovery. It develops guarantees for several covariance structures, including universal, stable, and near-minimal recovery, with applications to streaming and phaseless measurement settings.
- Motivation: High-dimensional data streams require covariance estimation in a single pass with minimal storage and low computational complexity.The paper motivates structural assumptions because acquisition devices may have limited memory and processing power relative to the data volume and rate.
- Applications and implementation: Quadratic covariance sketches use storage complexity m, linear per-instance sketching cost, and support distributed and asynchronous aggregation.Each randomized sketch is described as a compressive snapshot of second-order statistics, unlike uncompressed measurements that do not capture correlation information.
- Applications and implementation: The measurement scheme applies to streaming covariance estimation and compressive phase space tomography, where approximately low-rank correlation matrices can be recovered from few measurements.The paper also identifies broader phaseless-measurement applications for high-dimensional covariance or correlation structures.
- Measurement model and structures: The framework reconstructs low-rank, Toeplitz low-rank, sparse, and jointly sparse rank-one covariance matrices from a small number of quadratic measurements.Each quadratic measurement has the form a_i^⊤Σa_i plus noise, and the recovery algorithms use convex relaxations tailored to the presumed structure.
- Guarantees: Exact and universal recovery holds with high probability, while stable recovery remains accurate under imperfect structural assumptions and measurement noise.The stated stability guarantee bounds estimation error by a constant multiple of the noise level, including possibly adversarial noise.
- Guarantees: Near-minimal measurement guarantees are obtained for most structures, using RIP-ℓ2/ℓ1 for general models and RIP-ℓ2/ℓ2 for Toeplitz low-rank matrices.The paper states that the Toeplitz result follows by showing that linear combinations of quadratic measurements satisfy RIP-ℓ2/ℓ2 on the relevant class.
- Applications and implementation: Rank-one measurements are easier and cheaper to implement than full-rank random measurement matrices and require no additional covariance incoherence conditions.The claimed universality contrasts with standard matrix completion requirements for incoherence.
D. Organization
The paper develops convex recovery methods for structured covariance matrices from random quadratic measurements and establishes exact, universal, and stable guarantees for low-rank and Toeplitz low-rank cases.
- Problem and model: The analysis uses quadratic measurements from i.i.d. sub-Gaussian sensing vectors, with bounded measurement noise modeled in ℓ1 or ℓ2 norm.The measurements are represented by a linear operator acting on covariance matrices.
- Convex relaxation: Trace minimization replaces NP-hard rank minimization and provides a convex surrogate for PSD covariance matrices.The relaxation remains effective for approximately low-rank matrices and bounded measurement noise.
- Low-rank covariance recovery: m = Θ(nr) measurements yield exact recovery of rank-r PSD covariance matrices without noise, matching the intrinsic degrees of freedom order.The guarantee holds with exponentially high probability.
- Low-rank covariance recovery: The same program gives universal recovery for all low-rank covariance matrices once the sensing vectors are fixed, across a broad class of sub-Gaussian measurements.The result holds near the information-theoretic measurement limit.
- Low-rank covariance recovery: Approximately low-rank covariance matrices admit almost accurate estimates under natural power-law spectral decay without prior signal knowledge beyond that decay assumption.The noiseless reconstruction guarantee applies when m is about the order of nr.
- Stationary covariance recovery: Toeplitz rank-r covariance matrices have stable and universal recovery guarantees when m > c0r log10 n, under sub-Gaussian sampling, µ4 ≤ 3, and bounded ℓ2 noise.The Toeplitz requirement is substantially smaller than the general low-rank requirement and general Toeplitz degrees of freedom.
C. Recovery of Sparse Covariance Matrices
The paper recovers sparse covariance and jointly sparse rank-one matrices from quadratic measurements using convex relaxations, with exact, universal, and robust guarantees under sub-Gaussian sampling.
- Sparse covariance recovery: ℓ1 minimization replaces intractable support-size minimization for approximately sparse covariance matrices.The relaxation remains stable under bounded measurement noise and approximate sparsity.
- Sparse covariance recovery: m > c1k log(n2/k) measurements suffice for simultaneous recovery guarantees over all sparse covariance matrices.The bound uses the best k-sparse approximation and universal positive constants.
- Sparse covariance recovery: Exact k-sparse covariance matrices are recovered without noise with exponentially high probability at measurement order k log(n2/k).The guarantee is stated as optimal within a constant factor.
- Sparse covariance recovery: The same sensing mechanism universally recovers all sparse covariance matrices and remains robust for approximately sparse matrices.The paper states that quadratic measurements are order-wise at least as good as linear measurements in this setting.
- Jointly sparse and rank-one recovery: For jointly sparse rank-one matrices, a trace-plus-ℓ1 convex program balances low-rank and sparse structure while accommodating bounded noise.The regularization parameter controls the two convex surrogates.
- Jointly sparse and rank-one recovery: O(k2 log n) noise-free measurements universally recover all k-sparse signals with exponentially high probability.This guarantee extends prior Gaussian-sensing results to a broad class of sub-Gaussian sensing vectors using a simpler proof.
- Jointly sparse and rank-one recovery: Under imperfect sparsity or noisy samples, the recovered signal error is controlled by the structural tail and is at most proportional to per-entry noise.The signal estimate is obtained from the top normalized eigenvector of the lifted estimate, with Davis–Kahan perturbation analysis.
- Analysis: The analysis introduces RIP-ℓ2/ℓ1, measuring input strength by the Frobenius norm and output strength by the ℓ1 norm.This mixed-norm property supports universal recovery of low-rank, sparse, and sparse rank-one covariance matrices without dual-certificate constructions.
B. RIP-ℓ2/ℓ1 of Quadratic Measurements for Low-rank and Sparse Matrices
The paper removes measurement bias with an auxiliary operator and establishes mixed-norm restricted isometry for low-rank, sparse, and low-rank-plus-sparse matrix classes under sub-Gaussian sampling.
- Auxiliary operator: The original quadratic sampling operator fails RIP-ℓ2/ℓ1 because its measurement matrices have non-zero mean.The paper introduces debiased auxiliary measurement matrices to remove this bias.
- Measurement representation: The auxiliary formulation doubles the measurement count for notational simplicity without changing the order-wise results.This representation uses rank-2 measurements B_i instead of the original rank-one notation.
- Proof strategy: The auxiliary operator B exhibits RIP-ℓ2/ℓ1 with minimal measurements, established through a pointwise proposition and a covering argument.The proof combines Proposition 1 with standard covering techniques.
- Low-rank matrices: B satisfies RIP-ℓ2/ℓ1 for rank-at-most-r matrices when m > c4nr.The guarantee holds with probability exceeding 1 − C3 exp(−c3m).
- Sparse matrices: B satisfies RIP-ℓ2/ℓ1 for sparsity-at-most-k matrices when m > c4k log(n2/k).The result holds uniformly over all matrices in the sparse class with exponentially high probability.
- Low-rank-plus-sparse matrices: For low-rank-plus-sparse matrices, RIP-ℓ2/ℓ1 holds when the measurement count exceeds a universal-constant multiple of the class-specific maximum term.The class combines a rank-r matrix in a k × k subspace with an ℓ-sparse matrix.
C. Proof of Theorems 1, 3 and 4 via RIP-ℓ2/ℓ1
The proofs establish recovery guarantees through RIP-ℓ2/ℓ1 for low-rank, sparse, and jointly structured covariance matrices, while a separate RIP-ℓ2/ℓ2 route handles Toeplitz low-rank matrices.
- RIP-ℓ2/ℓ1 with sufficiently small constants is the common condition used to prove Theorems 1, 3, and 4.
- m > c4 (K1 + 2r) n measurements establishes the low-rank recovery guarantee in Theorem 1.
- m > c4 (K2 + 2k) log(n2/k) measurements establishes the sparse recovery guarantee in Theorem 3.
- Theorem 4 follows from a best k-term approximation condition and RIP-ℓ2/ℓ1 constants satisfying the stated inequality.
- General Toeplitz low-rank matrices require a separate RIP-ℓ2/ℓ2 analysis because RIP-ℓ2/ℓ1 covering-number arguments are not rigorously characterized for that case.
- For bounded and near-isotropic operators, RIP-ℓ2/ℓ2 yields universal stable recovery with m = Ω(nrpolylogn), strengthening guarantees for Fourier-type measurements.
B. Construction of RIP-ℓ2/ℓ2 Operators for Toeplitz Low-rank Matrices
The construction converts generally non-isotropic and poorly bounded quadratic measurements into truncated operators that are near-isotropic and well-bounded on Toeplitz matrices. This enables exact and stable recovery of rank-r matrices once the measurement count exceeds O(nrpolylogn).
- Construction: The original measurement matrices are generally non-isotropic and not well-bounded, motivating a new measurement construction.The paper introduces new matrices to facilitate application of the RIP-ℓ2/ℓ2 theorem.
- Construction: The construction defines rank-at-most-3 matrices, adds independently generated Gaussian matrices, and truncates the resulting matrices.The Gaussian matrices have independent standard Gaussian entries.
- Isotropy: A linear combination of the measurement matrices can be made isotropic when restricted to matrices with identical diagonal entries, including Toeplitz matrices.Lemma 4 supplies the isotropy property used for the associated operator.
- Boundedness: The operators formed before truncation are generally not well-bounded, but probabilistic norm bounds control the constituent matrices and their combinations.The stated bounds hold with probabilities exceeding 1−n^-10, 1−3n^-8, and 1−n^-7 in the cited steps.
- Recovery: The truncated operator is near-isotropic, and exact and stable recovery holds for all rank-r matrices when M exceeds O(nrpolylogn).This result establishes Theorem 2 through an equivalence argument.
C. Proof of Theorem 2
The paper develops and evaluates convex recovery of structured covariance matrices from quadratic measurements, with experiments indicating near-information-theoretic recovery and robustness to noise. It also identifies computational trade-offs and open sampling-model questions.
- Low-rank recovery: 20-trial experiments with n = 50 compare Gaussian and symmetric Bernoulli sensing vectors using empirical success probabilities over (m, r) pairs.Figure 2 overlays the information-theoretic limit as a red line.
- Low-rank recovery: The practical phase transition curve is very close to the theoretical sampling limit for low-rank covariance recovery.Recovery is declared when relative Frobenius error is below 10^-3.
- Low-rank recovery: Trace minimization experiments vary rank, measurement count, and bounded noise, reporting NMSE for n = 40.The noiseless and noisy cases are shown separately, with r = 5 in the noise experiments.
- Computational comparison: POCS requires more measurements than trace minimization to succeed, but has substantially lower computational cost.The comparison uses POCS results from 2000 iterations.
- Other structures: The framework extends across Toeplitz low-rank, sparse, and jointly rank-one and sparse covariance structures, with phase transitions near information-theoretic limits.The paper also reports stability under noise and imperfect structural assumptions.
- Proof and guarantees: The proof framework establishes recovery guarantees for low-rank covariance matrices from sub-Gaussian quadratic measurements using convex programming.The broader analysis uses RIP-based arguments for structured covariance models.
APPENDIX A PROOF OF PROPOSITION 1
The proof of Proposition 1 bounds the expected absolute quadratic measurement and establishes concentration by combining moment estimates, sub-exponential control, and Bernstein-type inequalities.
- Moment bounds: The proof first derives upper and lower bounds for E[|⟨B_i, X⟩|] before applying a Bernstein-type inequality.This establishes a large-deviation bound for the measurement process.
- Concentration: Hanson-Wright concentration shows that the quadratic form ⟨B_i, X⟩ is sub-exponential for sub-Gaussian sensing variables.Its sub-exponential norm is controlled by the Frobenius norm of X.
- Concentration: Bernstein’s inequality yields probability bounds exceeding 1 − 2 exp(−cmϵ) for deviations of averaged absolute measurements.The constants depend only on the sub-Gaussian norm of the sensing vectors.
- Uniform control: The proof uses symmetrization, Gaussian-process bounds, and entropy arguments to control RIP quantities uniformly over low-rank tangent spaces.The argument extends prior Pauli-measurement analysis to general near-isotropic measurements.
- Recovery implication: A small RIP-ℓ2/ℓ2 constant is then converted into recovery for all matrices of rank at most r.The proof concludes after establishing the required concentration and uniform bounds.
APPENDIX C PROOF OF LEMMA 1
The proof of Lemma 1 introduces the tangent space associated with a rank-r matrix and uses decompositions and norm inequalities to control recovery error components.
- Tangent-space setup: The proof defines the tangent space T from the singular-value decomposition of a rank-r matrix and decomposes perturbations into T and T⊥ components.PT and PT⊥ denote the corresponding orthogonal projections.
- Error decomposition: For an approximately low-rank matrix, the error is written as H around the target and constrained using trace-minimization optimality.The decomposition separates the best rank-r approximation from the residual.
- Norm control: Orthogonal components are organized into rank-controlled blocks to bound their nuclear norms relative to the tangent-space component.The proof uses singular-value ordering across successive blocks.
- Recovery bound: Feasibility bounds the measurement norm of the error, which is combined with the preceding decomposition to obtain the stated recovery inequality.The final bound depends on universal constants.
- Conclusion: The argument concludes by combining the derived inequalities with universal constants C1, C2, and C3.
APPENDIX D PROOF OF LEMMA 2
The proof of Lemma 2 analyzes perturbations of sparse rank-one matrices by combining support decompositions, tangent-space geometry, and nuclear- and ℓ1-norm control.
- Sparse decomposition: The proof separates the k largest support entries from the residual and decomposes the error across support and tangent-space components.The support projection is denoted PΩ.
- Block decomposition: The proof partitions off-support errors into blocks with controlled support size and decreasing entry magnitudes.This block structure enables norm comparisons for sparse recovery.
- Rank-one geometry: For rank-one structure, X is written as xx⊤ and its best k-term approximation defines the principal tangent space.The residual is Xc = X − XΩ.
- Recovery bound: Combining feasibility, tangent-space compatibility, and norm inequalities yields the claimed bound with universal constants.The final steps explicitly combine inequalities (82)–(86).
- Convex optimality: Optimality is handled with a subgradient combining nuclear and ℓ1 norms, including sign information outside the principal support.The construction uses matrices W and Y constrained in operator and infinity norms.
APPENDIX F PROOF OF LEMMA 4
The proof determines coefficients a, b, and c so that B is isotropic, derives the resulting form, and verifies the required discriminant condition.
- The proof seeks coefficients a, b, and c that make B isotropic.
- The isotropy constraint imposes E[B] = a + b + c = ϵ √n.
- Solving the resulting quadratic equation determines the coefficient values.
- The discriminant is positive when ϵ2 > 1.5 · |3 − µ4|.
- Using β = bα and γ = cα, the proof derives the form of Bi introduced in (39).
APPENDIX G PROOF OF LEMMA 5
The proof analyzes a symmetric Toeplitz matrix through its diagonal averages and a circulant embedding, then bounds the resulting quadratic-form eigenvalues with high probability.
- For symmetric Toeplitz M, each entry M_k equals the average of the corresponding constant descending diagonal.
- The diagonal averages satisfy E[M_0] = 1 and E[M_k] = 0 for 1 ≤ k < n.
- The proof embeds M into a (2n − 1) × (2n − 1) circulant matrix C_M and bounds C_M's spectral norm.
- The eigenvalues λ_i of C_M are quadratic forms in z_1, z_2, ···, z_n, with Eλ_i = EM_0 = 1.
- The resulting bound holds with probability exceeding 1 − 1/n^10, completing the proof together with (95).
APPENDIX H PROOF OF LEMMA 6
The proof combines quadratic-form tail bounds, Gaussian-process control, and nuclear-norm metric-entropy estimates to bound the relevant supremum over low-rank matrices.
- The argument begins by introducing auxiliary events while using that the restriction of B_i to Toeplitz matrices is isotropic.
- A tail inequality for quadratic forms supplies bounds beyond the threshold t > (20n log n)^2.
- The proof controls a Gaussian-process supremum using Dudley’s inequality.
- For normalized matrices with rank at most r, the pseudo metric is bounded through max_i |B_i(X − Y)|.
- The covering number of the nuclear-norm ball D1 is bounded using a standard metric-entropy procedure and the containment 2rD2^2r ⊆ D1^2r ⊆ D1.