Estimating Hierarchically Rank Structured Covariance Matrices

1 Introduction

• Spatial localization yields data-sparse, positive-definite estimates but requires multiscale modifications and sensitive bandwidth tuning; correlation-based localization is more general but difficult to store sparsely and rarely positive definite.The paper frames estimator design as balancing prior flexibility against computational constraints.Source: §1 · pp. 1–2
• The paper uses hierarchical rank structure to regularize covariance estimation, allowing long-range correlations to vary more smoothly than short-range ones rather than assuming they vanish.Unlike post-estimation compression, the structure is imposed within the estimation procedure.Source: §1 · p. 2
• HCov comes with an error bound for asymptotically smooth covariance structures; RSCov uses a factorization-based representation to guarantee positivity.The paper also reports experiments recovering these structures from far fewer samples than traditional methods.Source: §1 · p. 2

1.1 Code Availability

• The experiments can be reproduced using the authors’ publicly available code.The passage provides a release link for the code.Source: §1.1 · p. 2

2.1 The Need for Small-Sample Covariance Estimators

• Model reduction uses snapshot covariance to identify dominant variability directions, while ensemble data assimilation uses forecast covariance to propagate observations to unobserved state variables.Both settings rely on costly simulation snapshots or ensemble members.Source: §2.1 · p. 3
• Sample covariance estimates require enough samples to represent the covariance’s effective dimension, but simulation costs often leave applications below the O(D(C) log n) threshold.The passage defines effective dimension in terms of the trace and spectral norm of the covariance.Source: §2.1 · pp. 3–4

2.2 Regularizing Small-Sample Covariance Estimates

• Shrinkage regularizes the sample covariance spectrum toward a target, while thresholding removes small covariance entries and may fail to preserve positive semidefiniteness.Linear shrinkage toward the identity is one simple form; soft thresholding can be formulated to ensure positivity.Source: §2.2 · pp. 4–5
• Covariance sparsity assumes each component correlates with few others, whereas precision-matrix sparsity encodes conditional dependence and, for Gaussian variables, a graph of neighbors.The graphical lasso and row-wise regression methods are examples of precision-matrix estimators.Source: §2.2 · p. 5
• Spatial localization suppresses long-range correlations, while correlation-based localization adjusts sample correlations and can accommodate more general structures.Spatial localization bandwidth must be tuned; power-law correction is a correlation-based variant, and an even integer power guarantees positive semidefiniteness.Source: §2.2 · pp. 5–6

2.3 Hierarchical Rank Structure

• Hierarchical rank structure represents interactions between well-separated spatial domains with low-rank submatrices, reflecting interactions that vary more smoothly at long distances.The paper proposes using this structure as a covariance-estimation regularizer rather than only for later compression.Source: §2.3 · p. 6
• The error analysis is limited to low-dimensional settings: its dimension-dependent constants grow exponentially as d increases, and the authors know of no realistic applications with d ≥4.The restriction to d ≤3 is not technically necessary.Source: §2.3 · p. 6
• H-matrices organize admissible blocks in a tree, with rank-k representations on those blocks; their tree-based operations can have low asymptotic complexity.The block tree splits until blocks are admissible or reach the cluster-tree depth.Source: §2.3 · pp. 7–8

3.1 Problem Formulation

• Asymptotic smoothness of the covariance kernel is the key assumption that gives it hierarchical rank structure and permits close H-matrix approximation.The kernel is assumed analytic away from x = y and to satisfy the stated derivative bounds.Source: §3.1 · p. 9

3.2 General Framework

• HCov estimates covariance matrices as H-matrices by handling admissible cross-covariance blocks with low-rank factors and non-admissible blocks as unrestricted submatrices.The block-tree procedure runs O(n) estimators in parallel and enforces symmetry between paired blocks.Source: §3.2 · p. 9

3.3 Estimation via Polynomial Approximation

• Polynomial bases represent admissible cross-covariance blocks with few coefficients, and RegCC estimates each block by projecting its unregularized estimate into the chosen subspaces.The resulting estimate has rank at most the smaller of the two subspace dimensions.Source: §3.3 · p. 11

3.4 Error Analysis

• The error bound separates polynomial-truncation bias, which decreases exponentially with rank k, from sampling error, which scales with polynomial-subspace rank rather than total dimension n.The exponential rate depends on η; smaller η imposes greater smoothness on admissible blocks.Source: §3.4 · p. 12
• Sampling error also reflects unregularized non-admissible blocks, whose contribution depends on their total size, while η^-d and tree depth L further affect the bound.The η^-d factor reflects the number of cross-covariance leaves requiring separation of within-domain and between-domain variation; L reflects accumulated error across tree levels.Source: §3.4 · p. 13

4 Enforcing Positive Definiteness

• HCov may have negative eigenvalues, so it cannot guarantee exact Cholesky factorization or support the intended efficient sample generation.The section develops an SPD estimator while seeking to retain HCov's accuracy.Source: §4 · p. 13

4.1 Recursive Skeletonization

• Recursive skeletonization uses sparse, nonorthogonal transforms to reduce the matrix approximately to block-diagonal form, with small blocks and a residual matrix.The resulting factorization is SPD exactly when each diagonal block and the residual are SPD.Source: §4.1 · pp. 13, 15

4.2 Estimation Algorithm

• RSCov modifies HCov's recursive skeletonization and produces a nearby SPD covariance estimate in factorized form.The construction enforces positive definiteness on the small blocks and the final residual.Source: §4.2 · p. 17

5.1 Synthetic Covariance Kernel

• HCov and RSCov substantially outperform shrinkage, spatial localization, PLC, and the unregularized estimator on the synthetic multiscale kernel; RSCov remains close to HCov despite a small positive-definiteness correction cost.The kernel combines shortwave variations with longwave correlations, which the hierarchical methods filter without erasing the long-distance structure.Source: §5.1 · p. 20, Figure 5.3 · p. 21
• For the synthetic example, HCov and RSCov perform best at smaller subspace dimensions, while PLC overtakes HCov in VFIS at larger dimensions; RSCov’s spectrum has a plateau likely caused by diagonal modification.The unregularized estimator’s rank is limited by the sample size, which is 55 in this comparison.Source: §5.1 · p. 20, Figure 5.4 · p. 21

5.2 Climatological Covariance of Quasigeostrophic Flow

• On QG climatology, HCov and RSCov preserve long-distance correlations while filtering stochastic-forcing noise; RSCov storage is over an order of magnitude cheaper than the full matrix at small sample sizes.At larger sample sizes, RSCov storage is also cheaper than storing the samples; HCov shows comparable but less dramatic storage gains.Source: §5.2 · p. 23, Figure 5.6 · p. 23
• For QG climatology, HCov trails spatial localization in VFIS at moderate-to-large subspace dimensions, while RSCov outperforms HCov; the covariance-estimation comparison favors HCov and RSCov over localization.The eigenspace comparisons use 55 samples.Source: §5.2 · p. 23, Figure 5.7 · p. 24

5.3 Multiscale Lorenz Forecast Covariance

• For the multiscale Lorenz forecast covariance, algorithm comparisons follow the earlier examples, but the HCov–RSCov performance gap is more substantial and grows moderately at large sample sizes.The test system has short-wavelength x dynamics and longer-wavelength y dynamics.Source: §5.3 · p. 26

5.4 Surface Quasigeostrophic Forecast Covariance

• The SQG setup generates 8,192 nearly-identical states from a noise-perturbed, 100-day seed in a two-layer model with 64 × 64 buoyancy grids, then extracts first-layer vectors of dimension 4,096.The ensemble is integrated forward until its variance reaches 10% of its steady-state value.Source: §5.4 · p. 26
• The covariance is relatively sparse off the diagonal, favoring spatial localization, while its nonzero off-diagonal structures are more oscillatory and may degrade HCov and RSCov performance.The matrix is not completely zero off the diagonal.Source: §5.4 · p. 27
• HCov and RSCov use relatively small optimal approximation degrees at extremely low sample sizes, but storage grows quickly as more samples require higher degrees.Storage can exceed the O(n^2) cost of the full covariance matrix; the observed plateau likely reflects the experiment's maximum allowed approximation degree.Source: §5.4 · p. 27
• Spatial localization and PLC outperform HCov and RSCov in estimated-subspace quality, while spatial localization also outperforms them in covariance estimation.The subspace result is consistent with the covariance-estimation comparison.Source: §5.4 · p. 27

• Hierarchical rank regularization reduces estimation error for covariance matrices arising from asymptotically smooth covariance kernels.Source: §6 · p. 27
• HCov uses low-degree polynomial approximation with theoretical error analysis, while a modified recursive skeletonization factorization can enforce positive definiteness.Source: §6 · p. 27
• HCov and RSCov perform well on shortwave–longwave spatial processes, including examples without analytically verified asymptotic smoothness, but are not universally superior when simpler localization is appropriate.Source: §6 · p. 27

• The comparison uses a 50-point uniform-grid discretization of the 1D process, with regularization strength tuned optimally for each sample size.Source: A Comparison With the Graphical Lasso · p. 33
• Except at high sample sizes, every other estimator outperforms the graphical lasso in covariance accuracy, although it beats the unregularized estimator.Source: A Comparison With the Graphical Lasso · p. 33
• The graphical lasso assumes a sparse precision matrix, whereas this example has a non-sparse precision matrix associated with genuine long-range couplings between distant wave crests and troughs.Source: A Comparison With the Graphical Lasso · p. 33
• Despite weaker overall covariance accuracy, the graphical lasso performs comparatively well for spectrum and eigenspace estimation.Source: A Comparison With the Graphical Lasso · p. 33

• A thin SVD A = UΣV^T supplies orthonormal factors and a diagonal matrix of singular values.Source: B Proof of Frobenius Norm Identity · p. 33
• Orthogonal invariance makes z = V^Tx standard Gaussian, enabling the transformed expectation that completes the claim.Source: B Proof of Frobenius Norm Identity · pp. 33, 35
• The excerpt equates the expectation of tr(xTATAx) for x∈N(0,I) with the expectation of tr(zTΣ2z) for z∈N(0,I), without stating a starting identity for F.Source: B Proof of Frobenius Norm Identity · p. 35

• Differentiating k times yields sine terms for odd k and cosine terms for even k, with parity-dependent signs.Source: C Asymptotic Smoothness Proof for… · p. 35
• The derivative bound uses the kth Hermite polynomial and [41, eq. 8.22.8].Source: C Asymptotic Smoothness Proof for… · p. 35
• Equation (5.2) defines an asymptotically smooth kernel with growth parameter γ = 4π.Source: C Asymptotic Smoothness Proof for… · p. 35

• The proof is adapted nearly detail-for-detail from [6, lemma 3.15], with the scale factor changed to ∥C∥∞.Source: D Proof of Lemma 3.1 · p. 35
• Analyticity of C* supports a Taylor expansion of C(·, y) about a fixed point ξ in the cluster.Source: D Proof of Lemma 3.1 · p. 35
• Asymptotic smoothness bounds r(k), with the estimate applied to η-admissible cluster pairs.Source: D Proof of Lemma 3.1 · p. 36

• Lemma E.2 bounds spectral or Frobenius norm deviations for iid Gaussian vectors with covariance R, including singular R via factorization.Source: E Proof of Theorem 3.1 · pp. 36–37
• Theorem E.1 gives a local error bound for admissible leaf blocks under its stated restriction on η and γ.Source: E Proof of Theorem 3.1 · p. 37
• To control nearby cluster counts, the proof encloses nearby sidelength-ℓ cubes in a ball and bounds their number.Source: E Proof of Theorem 3.1 · p. 39
• Theorem E.2 states an error bound for estimator bCHM in terms of tree depth and the maximum bottom-level cluster size p.Source: E Proof of Theorem 3.1 · pp. 39–40
• Bottom-level blocks use a standard outer-product estimator, while higher-level admissible leaves use the local error bound.Source: E Proof of Theorem 3.1 · pp. 37, 40, 41
• The proof establishes eq. (3.5), bounding its first term by ε²∥C∥².Source: E Proof of Theorem 3.1 · p. 42

F Interpolative Decomposition Subroutine

• The decomposition selects a well-conditioned column subset whose span approximates the range of A.This is the column-selection task used in RSCov.Source: F Interpolative Decomposition Subroutine · p. 42
• Column-pivoted QR chooses the permutation through a greedy strategy for maximizing the minimum singular value of the selected columns.The procedure is specified in Algorithm F.1.Source: F Interpolative Decomposition Subroutine · p. 42
• The algorithm selects k as the largest index whose diagonal entry in R exceeds ϵ times the magnitude of R(1,1).The input is a relative error tolerance ϵ > 0.Source: F Interpolative Decomposition Subroutine · p. 42
• The interpolation matrix is computed as T = R(1:k,1:k)^-1 R(1:k,k+1:n), and the algorithm returns k, Π, and T.The output dimensions are k ≥ 1, Π ∈ {0,1}^(n×n), and T ∈ R^(k×n).Source: F Interpolative Decomposition Subroutine · p. 42

G Tuning Procedure for PLC Localization

• The PLC exponent β is tuned using a procedure inspired by NICE and a discrepancy rule comparing the sample covariance with the PLC estimate.The cited rule selects β using a tolerance δξ based on an a priori estimate of covariance error.Source: G Tuning Procedure for PLC Localization · pp. 42–43
• The search replaces integer β values with an infimum over [0, 10], approximated by fixed-iteration bisection.The original rule takes a minimum over positive integers.Source: G Tuning Procedure for PLC Localization · p. 43
• Instead of lookup tables, the method estimates ξ on the fly with a Monte Carlo procedure based on Fisher’s z-transformation.The estimate targets the Frobenius norm discrepancy on the sampled indices.Source: G Tuning Procedure for PLC Localization · p. 43
• The authors set δ = 1 and approximate the covariance discrepancy using 100 randomly chosen indices rather than evaluating the full norm.This replaces the full Frobenius norm evaluation with a sampled submatrix calculation.Source: G Tuning Procedure for PLC Localization · p. 43