Joint Mean–Covariance Estimation

• JMCE is a framework for simultaneous estimation of mean and covariance structures in multivariate, functional, and high-dimensional data, ensuring enhanced efficiency and interpretability.
• It leverages methodologies such as hierarchical Bayesian modeling, constrained optimization, and penalized likelihood to improve estimation accuracy and robustness.
• Applications span functional data smoothing, matrix-variate genomics, differential privacy in streaming data, and panel data analysis, providing practical solutions in complex settings.

Joint Mean–Covariance Estimation (JMCE) refers to a family of methodologies for simultaneous inference of both mean and covariance (or scatter) structures in multivariate, functional, or structured data settings. This paradigm departs from traditional approaches that estimate means and covariances (or precision matrices) sequentially or in isolation, enabling improved efficiency, robustness, and interpretability by properly accounting for their interdependence.

1. Hierarchical and Bayesian JMCE for Functional Data

JMCE for functional data—where each observational unit is a function or curve defined over a domain—addresses the challenge that existing methods can overly smooth individual curves and miss systematic cross-curve structure. A hierarchical Bayesian formulation posits observed noisy curves

$$y_i(t_{ij}) = Z_i(t_{ij}) + \epsilon_{ij},\quad \epsilon_{ij}\sim N(0,\sigma_\epsilon^2)$$

with latent true curves $Z_i(\cdot)$ modeled as i.i.d. GPs:

$$Z_i(\cdot) \sim GP(\mu(\cdot), \Sigma(\cdot,\cdot))$$

The mean function $\mu(\cdot)$ and covariance kernel $\Sigma(\cdot,\cdot)$ are inferred nonparametrically via second-level process priors:

• $$\mu(\cdot)\mid\Sigma \sim GP(\mu_0(\cdot),\Sigma(\cdot,\cdot)/c)$$
• $$\Sigma(\cdot,\cdot)\sim \mathrm{IWP}(\delta,\Psi(\cdot,\cdot))$$ (Inverse-Wishart process)

Posterior computation is performed by Gibbs sampling with all priors chosen to be (conditionally) conjugate, enabling cycling between Gaussian and inverse–Wishart full conditionals for all parameters and latent values. Critically, this model "borrows strength" across all curves via the shared mean and covariance, yielding improved estimation particularly in low-SNR or sparse settings. The automatic simultaneous smoothing also preserves systematic features instead of discarding them as noise. This framework accommodates data measured on irregular or uncommon grids without requiring data imputation or uniform sampling, and admits both stationary and nonstationary covariance structures depending on the choice of $\Psi$ (e.g., stationary Matérn kernel or empirical covariance). The joint posterior is

$$\begin{aligned} & p(Z,Z^*,\mu,\Sigma,\sigma_\epsilon^2,\sigma_s^2 \mid Y) \ &\quad\propto \Biggl[\prod_{i=1}^n \mathcal{N}\bigl(y_i(t_i)\mid Z_i(t_i),\,\sigma_\epsilon^2\,I_{p_i}\bigr)\Biggr] \;\times\; \mathcal{N}\bigl([Z;Z^*]\mid\mu,\Sigma\otimes I_n\bigr) \ &\quad\quad\times \mathcal{N}\bigl(\mu\mid\mu_0,\Sigma/c\bigr) \;\times\; \mathrm{IW}\bigl(\Sigma\mid\delta,\Psi(t,t)\bigr) \;\times\; \mathrm{IG}(\sigma_\epsilon^2\mid a_\epsilon,b_\epsilon) \;\times\; \mathrm{Ga}(\sigma_s^2\mid a_s,b_s). \end{aligned}$$

This formulation achieves lower root-integrated-MSE for both $Z_i$ and $\mu$ compared to conventional methods, and provides competitive or superior covariance estimation relative to established approaches such as PACE (Yang et al., 2014).

2. Mean–Covariance Estimation Under Constraints

In certain settings, mean and covariance must satisfy explicit joint constraints, notably $\Sigma \mu = \mu$ as arises in certain multivariate normal models. Two principal strands exist:

• Lagrange-multiplier-based maximum likelihood, with the Lagrangian

$$\mathcal{L}(\mu,\Sigma,\lambda) = \ell(\mu,\Sigma|X) + \lambda^\top(\Sigma\mu-\mu)$$

and non-linear score equations solved by Newton–Raphson or fixed-point iterations.

• Spectral reparameterization, reducing the parameter dimension by encoding the constraint into the covariance eigenstructure: one eigenvalue/unit-direction fixed by the mean, remainder as free parameters. This approach yields positive-definite, continuous estimators and enables efficient profile likelihood or Bayesian posterior computation via low-dimensional eigendecomposition.

Bayesian formulations use normal priors on the mean and inverse-gamma priors on eigenvalues, with tractable posterior factorization enabling Gibbs–Metropolis sampling or a fast approximate estimation via lower-bound concave maximization. These approaches drastically reduce the parameter count (from $p(p+1)/2$ to $O(p)$), enforce positivity, and attain competitive mean and covariance risk (Kundu et al., 2020, Kundu et al., 2021).

3. Joint Mean–Covariance Estimation for Matrix-Variate and High-Dimensional Data

For high-dimensional unreplicated matrix-variate data, such as genomics applications, JMCE involves modeling the matrix $X\in\mathbb R^{n\times m}$ as $X = M + \varepsilon$, with ${\rm vec}(\varepsilon) \sim N(0, A \otimes B)$. Simultaneous estimation is achieved via a penalized likelihood: $$\min_{M, \Theta\succ0, \Omega\succ0} \sum_{k=1}^m \operatorname{tr}\left[\Theta (X_k - M_k)\Omega (X_k - M_k)^\top\right] - m\log|\Theta| - n\log|\Omega| + \lambda_1\|\Theta\|_{1,\mathrm{off}} + \lambda_2\|\Omega\|_{1,\mathrm{off}}$$ Estimation proceeds by alternating blockwise: generalized least squares on the mean given precisions, and penalized graphical lasso on $\Theta$ and $\Omega$ based on current residuals. Theory provides consistent estimation and rates, with correct calibration of mean parameter inferences and improved power over conventional methods. The joint approach efficiently accounts for dependence graph structure among both samples and features, leading to more accurate mean estimation and properly controlled type I error in, e.g., large-scale differential expression studies (Hornstein et al., 2016).

4. Robust and Semiparametric Approaches for Elliptical and Non-Gaussian Distributions

JMCE has been adapted to semiparametric elliptical models beyond the multivariate Gaussian case. Under complex elliptically symmetric (CES) laws, the joint density admits a nuisance density generator, formally handled via semiparametric inference.

A robust, distribution-free procedure utilizes:

• Tyler's $M$-estimator for location and shape:

$$\sum_{i=1}^{N} q_i^{-1/2}(x_i - \mu) = 0, \quad \frac{p}{N} \sum_{i=1}^{N} (x_i - \mu)(x_i - \mu)^H / q_i = \Sigma$$

• Enhancement of efficiency via a rank-based $R$-estimator for the shape matrix, using ranks of squared Mahalanobis distances and projectile scores.

The semiparametric Cramér–Rao bound is block-diagonal: location and scatter can be estimated essentially independently, asymptotically achieving maximal efficiency in the semiparametric model. These methods are robust, $\sqrt{N}$-consistent, nearly attain the SCRB, and are computationally efficient (matrix operations of order $p^2$ per step) (Fortunati et al., 2021).

5. Nonparametric and Structured JMCE in High-Dimensional and Panel Data

JMCE methods for unbalanced panels (large-scale asset returns, time-varying panels) have been developed using RKHS and kernel-based learning. The COCO estimator formulates the problem as a convex empirical risk minimization using structured kernels to guarantee positive semidefiniteness and handle missing data naturally. Moment-matching loss includes both squared mean and cross-sectional covariance components, penalized by RKHS norms for regularization.

Low-rank approximations via Nyström or pivoted Cholesky factorization reduce computational burden, and convex quadratic objectives over positive semidefinite matrices enable efficient optimization. This estimator achieves $O(1/T)$ rates and exponential tail bounds, is robust to missingness and unbalanced cross-sections, and has delivered high out-of-sample Sharpe ratios and explanatory power for systematic risk in empirical financial data (Filipovic et al., 2024).

6. Joint Moment Estimation Under Differential Privacy and Streaming Settings

In privacy-constrained settings, JMCE has been extended to continuous data streams where both first and second moments must be estimated under $(\epsilon,\delta)$-differential privacy. The Joint Moment Estimation (JME) algorithm exploits the matrix mechanism to release both mean and covariance estimates at no extra privacy cost ("privacy for free" of the second moment). Key features:

• Workload-matrix representations allow general weighting (moving averages, sliding window),
• Shape-calibrated noise is added via Gaussian mechanism with carefully tuned joint sensitivity,
• Theoretical guarantees cover unbiasedness, privacy, and explicit Frobenius-norm error rates for both mean and moment estimates,
• Computational cost is competitive with or superior to baseline approaches.

Empirical gains are documented in both private Gaussian estimation and DP-Adam optimization, with reduced errors compared to postprocessing or baseline DP mechanisms (Kalinin et al., 10 Feb 2025).

7. Model-Specific and Application-Driven JMCE

Additional specialized JMCE methodologies address settings such as:

• Seemingly unrelated regression models under joint sparse mean and covariance estimation, leveraging horseshoe priors for improved shrinkage and support recovery in high dimensions (Li et al., 2019).
• Functional data with "snippets," where only short and irregular functional intervals are observed. Here, the covariance is decomposed into a nonparametrically estimated variance function and a parametric or basis-driven correlation surface, with measurement noise addressed through a novel local constant estimator (Lin et al., 2020).
• Linear models with uncertain signal and noise covariances, combining MAP and variational Bayes strategies with inverse-Wishart conjugacy for tractable fixed-point algorithms and joint estimation (Zachariah et al., 2014).

---

In summary, modern JMCE methodology encompasses hierarchical Bayesian, semiparametric, penalized likelihood, nonparametric kernel, robust, and privacy-aware approaches, targeting high-dimensional, functional, matrix-variate, streaming, and panel data. These methods deliver simultaneous inference for mean and covariance structures, often achieving both statistical optimality and computational scalability in challenging applied domains.