Multidimensional Rational Covariance Extension

• Multidimensional RCEP is a problem that reconstructs a nonnegative, rational spectral density from a finite set of covariance moments.
• It uses convex optimization and duality frameworks to model and identify multidimensional stationary processes, particularly for ARMA and shaping-filter applications.
• Recent advances, including efficient algorithms and matrix-valued extensions, have enhanced its application in image compression, radar, and medical imaging.

The Multidimensional Rational Covariance Extension Problem (RCEP) is a central inverse problem in systems theory, signal processing, and multidimensional spectral analysis. Its objective is to reconstruct a nonnegative power spectral density, with a rational absolutely continuous part, from a finite set of covariance moments. By enforcing a rational structure on the spectrum, the RCEP encapsulates both the modeling and identification of multidimensional stationary stochastic processes, especially for ARMA and shaping-filter models. The multidimensional setting introduces significant new mathematical and computational phenomena compared to the well-understood one-dimensional case, including unique geometric, convexity, and identifiability challenges. The RCEP is formalized as a truncated trigonometric moment problem (TMTMP) under rationality constraints, seeking measures or spectral densities rational in several complex variables, with applications ranging from image compression and texture generation to radar, sonar, and medical imaging. Recent advances have achieved well-posedness, effective parametrization, and robust solution algorithms, including connections to information divergence, convex optimization, and Riccati-type equations.

1. Formal Definition and Fundamental Structure

Let $\Lambda \subset \mathbb{Z}^d$ be a finite, symmetric index set ($0 \in \Lambda$, $-\Lambda = \Lambda$), and $c = (c_k)_{k\in\Lambda}$ be a Hermitian-symmetric vector of moments ($c_{-k} = \overline{c_k}$). Define the $d$-torus $\mathbb{T}^d = [-\pi,\pi]^d$. The moment-matching condition is

$$c_k = \int_{\mathbb{T}^d} e^{i\langle k, \theta\rangle} \Phi(\theta) \, d\theta, \quad k\in\Lambda,$$

where $d\theta = (2\pi)^{-d} d\theta_1\cdots d\theta_d$, and $\Phi \geq 0$ is the (unknown) spectral density.

The RCEP seeks to represent $\Phi$ as a rational function of the form

$$\Phi(\theta) = \frac{P(e^{i\theta})}{Q(e^{i\theta})}$$

where $P, Q$ are real Hermitian trigonometric polynomials, strictly positive on $\mathbb{T}^d$ except possibly on a set of measure zero where a singular part may be supported. In the multivariate and multidimensional case, analogously, the matrix-valued spectral density $\Phi$ is represented as $P(z)Q(z)^{-1}$, with $z_j = e^{i\theta_j}$.

The moment problem thus becomes: Given $c = (c_k)_{k\in\Lambda}$, find $\Phi = P/Q$ (with $P, Q \in \mathcal{P}_+$, the cone of positive trigonometric polynomials) and an optional singular measure $d\nu$ supported where $Q = 0$, such that the moment equations are satisfied exactly or approximately.

2. Existence, Uniqueness, and Parametric Structure

Existence and uniqueness of solutions critically depend on both the data $c$ and the dimension $d$. The dual cone

$$\mathcal{C}_+ = \{ c : \langle c, p \rangle > 0 \ \forall P \in \overline{\mathcal{P}}_+ \setminus \{0\} \}$$

characterizes feasible moment sequences. For $c \in \mathcal{C}_+$ and $P \in \overline{\mathcal{P}}_+ \setminus \{0\}$, the strictly convex dual problem

$$\min_{Q>0} \langle c, q \rangle - \int_{\mathbb{T}^d} P \log Q$$

admits a unique minimizer $Q^\ast$, and the rational spectrum $P/Q^\ast$ (plus possibly a singular part) solves the RCEP (Ringh et al., 2015).

For dimensions $d \leq 2$ and strictly positive $P$, the minimizer $Q^\ast$ is strictly positive, and the solution is purely absolutely continuous with $\Phi = P/Q^\ast$ (Ringh et al., 2015). In higher dimensions, the positivity of $Q$ can fail on a set of positive codimension; uniqueness of the singular part is generally lost, and only the matched moments are uniquely attained.

A complete smooth parametrization of all rational solutions is possible (Ringh et al., 2015). Fixing $P$ within the interior of the positive cone, the mapping $P \mapsto Q^\ast$ is a global diffeomorphism. The parameter-to-moment map is also a diffeomorphism, as shown for specific basis choices and convex optimization schemes (Wu et al., 13 Jan 2026).

3. Convex Optimization and Duality Frameworks

Modern approaches uniformly recast the RCEP as a convex optimization problem over the coefficients of $P$ and $Q$:

• Primal Formulation: Maximize a generalized entropy subject to moment constraints,

$$\max_{\Phi \geq 0} \int_{\mathbb{T}^d} P(e^{i\theta}) \log \Phi(e^{i\theta})\, d\theta, \quad \text{subject to} \quad c_k = \int e^{i\langle k,\theta \rangle} \Phi(\theta)\, d\theta.$$

• Dual Formulation: Minimize a convex cost in the denominator parameters,

$$J_P(q) = \langle c, q \rangle - \int_{\mathbb{T}^d} P \log Q, \quad Q(e^{i\theta}) = \sum_{k\in\Lambda} q_k e^{-i\langle k,\theta \rangle}.$$

Extensions to generalized moment problems, including cepstral moment constraints, yield joint optimization over $(P, Q)$, with dual objectives of the form

$$J(P, Q) = \langle c, q \rangle - \langle m, p \rangle + \int_{\mathbb{T}^d} P \log \frac{P}{Q}$$

with $p_0 = 1$ fixing normalization (Zhu et al., 2021).

Convexity and strict positivity of the Hessian ensure well-posedness and strong duality. Regularization (e.g., penalizing $\int P^{1-\nu}$ or $\int Q^{1-\nu}$) enforces interiority and enables robust solution under approximate or noisy moment sequences (Zhu et al., 2021).

4. Approximate Covariance Matching and Statistical Properties

In the presence of noise or when empirical covariances fall outside the feasible moment cone, exact moment matching is infeasible. Two principled approximate schemes are formulated (Ringh et al., 2017):

• Soft-constraint formulation: Minimizes a weighted sum of Kullback–Leibler divergence and moment deviation,

$$\min_{\Phi \geq 0,\, d\nu \geq 0} D(P\|\Phi) + \tfrac{1}{2}(r-c)^* W^{-1} (r-c),$$

with $D(P\|\Phi)$ the extended Kullback–Leibler divergence.

• Hard-constraint formulation: Minimizes entropy subject to a hard bound on moment deviation,

$$\min_{\Phi \geq 0,\, d\nu \geq 0,\, (r-c)^* W^{-1} (r-c) \leq 1} D(P\|\Phi).$$

These variants are linked by a homeomorphism between their regularization (weight) matrices. Both variants admit unique solutions, continuous dependence on data, and are solved by finite-dimensional strictly convex minimization (Ringh et al., 2017).

Statistically, with an explicit basis and block-Toeplitz moment representation, estimators can achieve consistency, (asymptotic) unbiasedness, optimal $O_P(N_k^{-1/2})$ convergence rates, and (under Gaussianity) efficiency (Wu et al., 13 Jan 2026). Efficient computation exploits the Kronecker structure of multidimensional trigonometric bases.

5. Matrix-Valued and Multivariate Extensions

The matrix-valued and multivariate versions of RCEP generalize the scalar problem to nonnegative Hermitian matrix-valued spectra. The essential solution methods include:

• τ-divergence approach: For prior $\Psi(e^{i\theta})$ and integer parameter $\nu\geq 2$, the minimum of

$D_{1-1/\nu}(\Phi || \Psi)$

under moment constraints, always yields rational optimal spectra with the closed-form (Zhu et al., 2020)

$$\Phi_\nu^*(\theta) = [\Psi^{-1}(\theta) + Q(\theta)/\nu ]^{-\nu}\Psi(\theta)^{\nu - 1}.$$

Rationality is assured for sufficiently large $\nu$, and the dual problem is strictly convex. This framework bridges from Itakura–Saito (ν = 1) to weighted Kullback–Leibler (ν → ∞) solutions.

• Shaping-filter parameterization: Once the rational spectrum is computed, multidimensional spectral factorization (e.g., Fejér–Riesz, Wilson's algorithm on block-circulants) produces shaping filters $W(z)$ such that $\Phi = W(z) W(z^{-1})^*$ for system identification and simulation (Zhu et al., 2021).

6. Riccati-Type and Analytic Interpolation Perspectives

An alternative solution leverages rational analytic interpolation via the Covariance Extension Equation (CEE) (Cui et al., 2020, Cui et al., 2024). Here, the RCEP is recast as finding a rational Carathéodory function $F(z)$ positive-real on the unit disk, matching prescribed Taylor or moment coefficients, with a McMillan degree bound. The spectral density is then $F(z) + F(z)^*$.

CEE is a Riccati-type matrix equation parameterized by spectral-zero parameters or shaping matrix polynomials. In the scalar and block cases, all rational positive-real interpolants are generated via the family

$P = \Gamma (P - P H^*H P) \Gamma^* + G(P)G(P)^*,$

where the problem data enter through structured state-space matrices. The CEE parametrization is a (local) diffeomorphism under mild conditions, and homotopy continuation algorithms yield all solutions, supporting model reduction via the rank–degree relationship ($\mathrm{rank}\,P = \deg \Phi_+$) (Cui et al., 2020, Cui et al., 2024).

7. Computational Algorithms and Applications

Solution of the multidimensional RCEP involves the following workflow (Wu et al., 13 Jan 2026, Ringh et al., 2017):

1. Form sample covariance moments (biased or unbiased), construct block-Toeplitz moment matrices, and select basis (often the multilinear Kronecker product basis).
2. Choose a prior $P$ and formulate the dual convex minimization problem.
3. Employ Newton-type, BFGS, or interior-point methods, leveraging multidimensional FFTs for integrals and gradients.
4. Recover the optimal spectral form, $P/Q$, and, if desired, spectrally factorize for ARMA/system identification.
5. For approximate data, select between soft/hard constraint formulations and tune regularization or weight parameters.

Applications include multidimensional spectral estimation in ARMA fields, texture synthesis via Wiener-type systems, robust multi-dimensional system identification, and image compression through sparse cepstral/covariance representations (Ringh et al., 2015, Ringh et al., 2017, Wu et al., 13 Jan 2026).

8. Open Problems, Limitations, and Generalizations

Several open technical and conceptual issues remain:

• In $d > 2$, rational solutions need not be purely absolutely continuous; support of singular measures may have nontrivial topology (Ringh et al., 2015).
• Uniqueness failures have been demonstrated in the multivariable THREE-parametrization with non-scalar priors, as the moment map can have bifurcation points and lose injectivity (Zhu, 2018).
• Spectral factorization in higher dimensions can be obstructed by rank conditions (Geronimo–Woerdeman); practical algorithms rely on block-circulant/FFT methods.
• The parameter–moment map's injectivity and global diffeomorphic properties depend intricately on the basis and regularity of the chosen parameter domain (Wu et al., 13 Jan 2026, Zhu et al., 2020).
• Model reduction via Riccati rank-deficiency and the precise characterization of minimal positive-real degree given a covariance sequence (Kalman's question) remain active research directions (Cui et al., 2020).

Continued development focuses on robust regularization, computational scalability, statistical optimality, and expansion to broader classes of moment/covariance constraints, with a deep interplay between convex geometry, spectral theory, and information divergences across all generalizations of the Multidimensional Rational Covariance Extension Problem.