PaperlayerView exact revision on arXiv

Source-linked AI summary

Flexible covariance estimation in graphical Gaussian models

Bala Rajaratnam, Hélène Massam, Carlos M. Carvalho

arXiv:0901.3267v1math.STstat.TH

TL;DR

Covariance estimation is difficult when the dimension is large relative to the sample size, and graphical models offer a way to encode structure. The paper derives closed-form Bayes estimators using flexible conjugate priors; the estimators show good frequentist risk properties and eigenvalue shrinkage in the reported examples.

  • Problem

    Covariance estimation is challenging in high dimensions: when n < r the sample covariance is not positive definite, and its eigenstructure can be distorted even when n > r.

    Source: §1 · p. 1

  • Method

    For decomposable-graph Gaussian models, the paper uses the conjugate WPG family to derive closed-form Bayes estimators under entropy and graph-adapted squared-error losses.

    Source: §1 · p. 2, §3 · p. 6, §5.1 · p. 16

  • Results

    The authors report that their Bayes estimators have good frequentist risk properties and shrink eigenvalues in numerical examples.

    Source: §8 · p. 26

  • Takeaways & Limitations

    The WPG family provides k + 1 shape parameters, allowing more flexible prior specification than the single-shape inverse hyper inverse Wishart.

    Source: §1 · p. 2, §3 · p. 6

    Under the graph-adapted squared-error loss, nonfree covariance entries must not contribute to the loss for the posterior mean to be Bayes optimal.

    Source: §5.1 · p. 16

Abstract

from arXiv · show

In this paper, we propose a class of Bayes estimators for the covariance matrix of graphical Gaussian models Markov with respect to a decomposable graph $G$. Working with the $W_{P_G}$ family defined by Letac and Massam [Ann. Statist. 35 (2007) 1278--1323] we derive closed-form expressions for Bayes estimators under the entropy and squared-error losses. The $W_{P_G}$ family includes the classical inverse of the hyper inverse Wishart but has many more shape parameters, thus allowing for flexibility in differentially shrinking various parts of the covariance matrix. Moreover, using this family avoids recourse to MCMC, often infeasible in high-dimensional problems. We illustrate the performance of our estimators through a collection of numerical examples where we explore frequentist risk properties and the efficacy of graphs in the estimation of high-dimensional covariance structures.

1. Introduction.

Covariance estimation becomes especially challenging when the dimension is large relative to the sample size, motivating structured alternatives to the sample covariance. The paper combines graphical-model regularization with Bayesian shrinkage and develops closed-form estimation using flexible conjugate priors.

1. Introduction.

  • When n < r, the sample covariance is not positive definite, and even for n > r its eigenstructure can be distorted unless r/n is extremely small.The introduction frames these difficulties in high-dimensional settings where the number of samples is often much smaller than the number of parameters.Source: §1 · p. 1
  • The proposed approach combines graphical-model regularization with Bayesian shrinkage for covariance and precision estimation under decomposable graphs.The resulting precision matrices have zeros in the pattern prescribed by the graph.Source: §1 · p. 2
  • The inverse WPG family is conjugate and includes the Wishart for complete graphs and the inverse hyper inverse Wishart for decomposable graphs, while offering k + 1 shape parameters.The extra shape flexibility addresses the single-shape-parameter restriction of the Wishart and inverse hyper inverse Wishart.Source: §1 · p. 2
  • The paper derives a noninformative reference prior for covariance and precision, with closed-form posterior estimation because the covariance posterior remains in the IWPG family.It also derives Bayes estimators under two graph-adapted loss functions and compares estimator performance in high-dimensional examples.Source: §1 · p. 3

2. Preliminaries.

The preliminaries describe graphical Gaussian models through conditional independences encoded by precision-matrix zeros and set notation for decomposable graphs. They also establish the covariance completion and likelihood-estimation conventions used later.

2. Preliminaries.

  • For a decomposable graph, a perfect clique order defines each new layer R_j and its separator S_j as the intersection with the preceding cliques.The separators are minimal, and repeated separators are counted by their multiplicity ν(S).Source: §2 · p. 3
  • In a graphical Gaussian model, variables at non-neighboring vertices are conditionally independent given the others, represented by zeros in the precision matrix.The covariance parameter is represented by the entries of the inverse precision matrix on the graph's specified positions.Source: §2 · pp. 3–4
  • A partial covariance matrix on the graph has a unique positive-definite completion when the graph is decomposable.The passages also describe a bijection between the covariance and precision parameter spaces.Source: §2 · p. 4
  • The graphical-model maximum likelihood estimate of the precision matrix exists when n exceeds the largest clique size; in the saturated model, the maximum likelihood estimate is S^-1.The paper distinguishes mle for the saturated model from mleg for the graphical model and uses tilde notation for estimators.Source: §2 · pp. 4–5

3. Flexible conjugate priors for Σ and Ω.

For decomposable graphs, the WPG/IWPG conjugate family extends the hyper inverse Wishart’s one-dimensional shape parameterization. Explicit posterior means support layerwise shrinkage, with different intensities across parts of the covariance matrix.

3.1. The WPG distribution and its inverse.

  • The WPG-Wishart family gives precision-matrix priors a k + 1-dimensional shape parameter, versus the one-dimensional shape parameter of the hyper inverse Wishart.The hyper inverse Wishart is a special case of the IWPG family.Source: §3 · pp. 5–6, §3.1 · p. 7
  • The explicit IWPG mean formulas let the posterior mean of the covariance matrix be reconstructed layer by layer.The posterior mean follows by updating the scale parameter with the observed data and adjusting the shape parameters.Source: §3.1 · pp. 8–9

3.2. Shrinkage by layers and the choice of the scale parameter θ.

  • Choosing the scale parameter θ sets the shrinkage target, while the shape parameters allow different shrinkage intensities across matrix layers.The posterior mean of each layer is a linear combination of its sample block and the corresponding scale block.Source: §3.2 · pp. 9, 11
  • With D equal to the average sample eigenvalue times the identity, posterior eigenvalues shrink toward that average when ν > r − 1, reducing their span.For eigenvalues below the average, the posterior values increase; above-average values shrink when ν exceeds the stated threshold.Source: §3.2 · p. 10

4. Reference prior.

The paper derives a reference prior by reparametrizing the graphical Gaussian model using an ordered Choleski parameterization. The induced prior is IWPG, and the posterior mean approaches the precision maximum-likelihood estimate as sample size grows.

4. Reference prior.

  • The reference prior is derived for an ordered Choleski parameterization induced by a perfect clique order, then induced on the covariance parameter.The reference prior for the ordered parameter is independent of the order of its components.Source: §4 · pp. 11–12
  • The induced prior on the covariance parameter has IWPG form, and its posterior remains a proper IWPG when the stated moment-existence conditions hold.The additional posterior hyperparameter conditions ensure that the moments used in the formulas exist.Source: §4 · pp. 12–13
  • As n grows, the posterior mean expression becomes close to the maximum-likelihood estimate of the precision matrix..Source: §4 · p. 13

5. Decision-theoretic results.

The paper adapts decision-theoretic covariance and precision estimation to decomposable graphical models by modifying losses to reflect their structural zeros. It derives Bayes estimators under these losses and establishes a risk property for the graphical maximum-likelihood estimator.

5. Decision-theoretic results.

  • The paper derives Bayes estimators of covariance and precision under entropy and graph-adapted squared-error losses for HIW, IWPG, and reference priors.The entropy-loss covariance estimator is the inverse of the posterior expected precision.Source: §5 · p. 13, §5.1 · pp. 15–16

5.1. Bayes estimators for Σ and Ω.

  • The precision matrix is the graphical model’s canonical parameter and has fewer parameters than the full positive-definite matrix space.Its structural zeros mean the adapted losses depend only on nonzero precision elements.Source: §5.1 · p. 14
  • For decomposable graphs, the squared-error loss is modified to depend on free graphical parameters, and posterior means are Bayes-optimal only under this graph-adapted loss.The traditional saturated-model losses must be reconsidered because the precision matrix has fixed zeros.Source: §5.1 · pp. 14–16
  • The graphical maximum-likelihood covariance estimator is the best entropy-loss estimator among its scalar multiples.The risk is minimized when the scalar multiplier is 1/n.Source: §5.1 · p. 17

6. Risk comparisons and numerical properties.

The section compares Bayes covariance and precision estimators by simulation, examining flexible shrinkage across cliques and selection among candidate graphs. The examples report lower risk for flexible-prior estimates in the two-clique setting and more parsimonious graph choices at smaller sample sizes.

6. Risk comparisons and numerical properties.

  • Across simulations, the study compares frequentist risks under L1 and L2 losses for covariance and precision estimates, alongside predictive and eigenvalue properties.In the two-clique example, the authors report that δ_i = 1/(4c_i) performs best and is a compromise among shrinkage choices.Source: §6 · p. 18, §6.1 · p. 19

6.1. Example 1: “Two Cliques.”

  • The two-clique simulation uses 100 variables and clique sizes 70 and 40, with different shrinkage levels intended for the different clique dimensions.The covariance is constructed from call-center sample covariance entries retained by the graph and completed to preserve its conditional independences.Source: §6.1 · p. 18
  • For L1(Ω) at n = 75, IW PG(1/4ci,D) reduces risk by 78% versus HIW (3,I) and 97% versus mleg.The reductions are 76% and 87% at n = 100, 20% and 30% at n = 500, and 9% and 16% at n = 1000, respectively.Source: §6.1 · p. 19

6.2. Example 3: choosing the graph.

  • The graph-selection procedure chooses among banded decomposable graphs by marginal likelihood, using the selected graph for covariance estimation.With a uniform graph prior, posterior graph probabilities are equivalent to marginal likelihood comparisons; the example uses HIW (3,I).Source: §6.2 · pp. 21–22
  • For n = 100, the selected average band size is 4.36, rising to 7.38 at n = 500 and 18.80 at n = 1000.At n = 100 and n = 500, oracle losses are always larger than those of the estimators; the oracle is relatively competitive only at n = 1000.Source: §6.2 · p. 22

7. Call center data.

The call-center analysis compares Bayesian graphical-model estimators with standard and banded alternatives, then tests whether allowing different band widths improves prediction. The results favor differential banding over the best uniformly banded model and show different selected structures for the two halves of the day.

7. Call center data.

  • The analysis uses 239 days of call-center records, with 102 ten-minute intervals per day, and compares graphical estimators for prediction.The data cover 7:00 am to midnight, with weekends, public holidays, and days with equipment malfunction excluded.Source: §7 · p. 22

7.1. Analysis with banded precision matrices.

  • Cross-validation selected k = 4 for all eight estimators, while Bayesian model selection selected k = 5; both procedures yielded similar parsimonious models.The forecast task predicted the second half of the day from the first.Source: §7.1 · p. 23
  • All Bayesian estimators had lower prediction errors than the mleg and mle, and their eigenvalue estimates shifted inward relative to the sample covariance estimator S.The Bayesian estimates were lower for the largest eigenvalues and higher for the smallest.Source: §7.1 · p. 24

7.2. Analysis with differential banding.

  • Cross-validation predicting the first half of the day from the second favored AR(16) (k = 16), suggesting one k may not suffice and that correlation structures differ across the day.The differentially banded model (k1 = 14, k2 = 4, r = 58) had the lowest prediction error when using IW PG priors.Source: §7.2 · p. 25
  • 16% improvement over the best k-banded model was reported for differential banding, while k-banding reduced prediction error 23% over the sample covariance matrix S.The differential-band model also predicted better than k-banded models for almost all time points.Source: Figure 2 · p. 24, §7.2 · p. 26
  • For predicting the second half from the first, the IW PG L1 estimator selected (k1 = 14, k2 = 1, r = 55) as the lowest-error differentially banded model.The authors frame graphical models as a more flexible alternative to simple banding.Source: §7.2 · p. 26

8. Discussion.

The discussion presents flexible-prior Bayes estimators for graphical covariance and precision matrices, along with a reference prior. It reports favorable frequentist risk properties and eigenvalue shrinkage, while identifying unresolved asymptotic and eigenvalue questions.

8. Discussion.

  • The approach uses flexible conjugate priors with multiple shape parameters to address limitations of one-parameter conjugacy and MCMC in very high-dimensional settings.The authors describe those existing approaches as potentially inadequate or infeasible in such settings.Source: §8 · p. 26
  • The paper derives Bayes estimators under two loss functions adapted to graphical models and a noninformative reference prior for Σ and Ω.The authors also report good frequentist risk properties and shrinkage in the eigenvalues.Source: §8 · p. 26
  • The authors identify the estimators’ asymptotic properties when both r and n grow, and their eigenvalue behavior, as questions for future work.They describe further assessment of estimator properties as still needed.Source: §8 · p. 26

A.1. Proof of Theorem 3.1.

The appendix develops calculations used to prove Theorem 3.1, including a matrix-normal lemma and blockwise expectation identities. It uses these results in deriving the mean of the IW PG distribution.

A.1. Proof of Theorem 3.1.

  • The appendix sets up the proof of Theorem 3.1 by recalling the matrix-normal distribution and establishing a lemma for matrix expectations.The stated matrix-normal representation uses covariance σ1 ⊗ σ2.Source: §A.1 · pp. 26–27
  • The lemma relates XaX^t to a sum over columns and uses the vectorized matrix-normal distribution to derive block expectation identities.The passages refer to the resulting identities as (A.2) and (A.3).Source: §A.1 · p. 27
  • The proof concludes after verifying the lemma’s stated identities by inspection.Source: §A.1 · p. 27

Proof of Theorem 3.1.

This proof computes entries of the IW PG matrix expectation using distributional properties of its subblocks. It proceeds layer by layer, following the dependence structure among entries.

Proof of Theorem 3.1.

  • The proof computes expected values of the different entries of an IW PG matrix using distributional properties of its subblocks.It relies on Theorem 4.4 of.Source: Proof of Theorem 3.1 · pp. 27–28
  • Because the entries have different dependencies, the expectation is computed in an ordered, layer-by-layer sequence.The text specifies starting from the top layer.Source: Proof of Theorem 3.1 · pp. 28–29
  • The calculation invokes the general expectation formula for an inverse Wishart distribution and applies similar calculations to later layers.Source: Proof of Theorem 3.1 · pp. 28–29

A.2. Proof of Theorem 4.1.

The proof derives the reference prior through likelihood and Fisher-information calculations, then transforms it using Jacobians. The hyper Wishart expectation simplifies one step, and the prior is independent of component order.

A.2. Proof of Theorem 4.1.

  • Differentiating the log-likelihood twice yields a block-diagonal Fisher information matrix, with the indicated blocks also block diagonal.The proof gives the determinant of these blocks for j = 2,...,k.Source: §A.2 · p. 30
  • For the hyper Wishart case, E(X) is proportional to σ, avoiding the general WQG expectation expression in (4.21) of.The proof uses this relation in the subsequent calculation.Source: §A.2 · p. 30
  • The reference prior for φ follows from natural exponential family theory and is independent of the order of its components.The proof invokes and Section 2.3 of.Source: §A.2 · p. 31
  • The induced prior for σ is obtained through Jacobians from φ to σ̂⁻¹ and from σ̂⁻¹ to σ.The proof concludes by combining these Jacobians to establish Theorem 4.1.Source: §A.2 · pp. 31–32
Loading 0901.3267v1…