A Weaker Regularity Condition for the Multidimensional $\nu$-Moment Problem

TL;DR

This paper weakens the regularity condition for solving the multidimensional $ u$-moment problem, enabling approximate solutions with smaller $ u$ values, thus broadening the applicability of spectral density reconstruction.

Contribution

The paper demonstrates that the regularity condition for the $ u$-moment problem can be relaxed from $ u \\geq d/2+1$ to $ u \\geq d/2$, allowing for more flexible spectral density solutions.

Findings

Regularity condition can be weakened to \\nu \\geq d/2.

Broader class of spectral densities can be approximated.

Potential for improved spectral density estimation methods.

Abstract

We consider the problem of finding a $d$-dimensional spectral density through a moment problem which is characterized by an integer parameter $\nu$. Previous results showed that there exists an approximate solution under the regularity condition $\nu\geq d/2+1$. To realize the process corresponding to such a spectral density, one would take $\nu$ as small as possible. In this letter we show that this condition can be weaken as $\nu\geq d/2$.

Equations

$${\mathbb{H}}_{\nu }(\Phi ):=\frac{{\nu }^2}{\nu -1}\left({\int }_{{\mathbb{T}}^d}{\Phi }^{\frac{\nu -1}{\nu }}\mathrm{d}\mu -1\right)$$

$$\underset{\Phi \ge 0}{\max }$$

$$c_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i⟨\mathbf{k},\mathbf{\theta }⟩}\Phi \mathrm{d}\mu \forall \mathbf{k}\in \Lambda ,$$

$$m_{\mathbf{k}}=\frac{\nu }{\nu -1}{\int }_{{\mathbb{T}}^d}{e}^{i⟨\mathbf{k},\mathbf{\theta }⟩}{\Phi }^{\frac{\nu -1}{\nu }}\mathrm{d}\mu \forall \mathbf{k}\in {\Lambda }_0.$$

$${\Phi }_{\nu }=(P/Q{)}^{\nu }$$

$$\frac{\lambda }{\nu -1}{\int }_{{\mathbb{T}}^d}\frac{1}{P^{\nu -1}}\mathrm{d}\mu$$

$$\nu \ge d/2+1,$$

$$\frac{\lambda }{\nu -1}{\int }_{{\mathbb{T}}^d}\frac{1}{P^{\nu }}\mathrm{d}\mu$$

$$J_{\nu }(P,Q):=\frac{1}{\nu -1}{\int }_{{\mathbb{T}}^d}\frac{P^{\nu }}{Q^{\nu -1}}\mathrm{d}\mu +⟨\mathbf{q},\mathbf{c}⟩-⟨\mathbf{p},\mathbf{m}⟩,$$

$$Q(e^{i\mathbf{\theta }}):=\sum _{\mathbf{k}\in \Lambda }{q}_{\mathbf{k}}{e}^{-i⟨\mathbf{k},\mathbf{\theta }⟩}$$

$$P(e^{i\mathbf{\theta }})=\sum _{\mathbf{k}\in \Lambda }{p}_{\mathbf{k}}{e}^{-i⟨\mathbf{k},\mathbf{\theta }⟩}$$

$$J_{\nu ,\lambda }(P,Q):=J_{\nu }(P,Q)+\frac{\lambda }{\nu -1}{\int }_{{\mathbb{T}}^d}\frac{1}{P^{\nu }}\mathrm{d}\mu$$

$$=\frac{1}{\nu -1}{\int }_{{\mathbb{T}}^d}g(P(e^{i\mathbf{\theta }}),Q(e^{i\mathbf{\theta }}))\mathrm{d}\mu +⟨\mathbf{q},\mathbf{c}⟩-⟨\mathbf{p},\mathbf{m}⟩$$

$$g(x,y):=x^{\nu }/y^{\nu -1}+\lambda /x^{\nu }$$

$${\mathfrak{P}}_{+}$$

$${\mathfrak{P}}_{+,o}$$

$$\nabla g(x,y)=$$

$${\nabla }^{2}g(x,y)=$$

$$+\left[\begin{array}{cc}\lambda \nu (\nu +1)x^{-\nu -2}& 0\\ 0& 0\end{array}\right],$$

$$\min J_{\nu ,\lambda }(P,Q)\text{s.t.}P\in {\overline{\mathfrak{P}}}_{+,o},Q\in {\overline{\mathfrak{P}}}_{+}.$$

$$c_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i⟨\mathbf{k},\mathbf{\theta }⟩}{\Phi }_0\mathrm{d}\mu \forall \mathbf{k}\in \Lambda ,$$

$$J_{\nu ,\lambda }^{-1}(-\infty ,\beta ]:=\{(P,Q)\in {\overline{\mathfrak{P}}}_{+,o}\times {\overline{\mathfrak{P}}}_{+}:J_{\nu ,\lambda }(P,Q)\le \beta \}.$$

$$f(t):=J_{\nu ,\lambda }(P,Q+t1)$$

$$f^{\prime }(t)=c_0-{\int }_{{\mathbb{T}}^d}{\left[\frac{P}{Q+t1}\right]}^{\nu }\mathrm{d}\mu .$$

$${\int }_{{\mathbb{T}}^d}{\left[\frac{P}{Q+t1}\right]}^{\nu }\mathrm{d}\mu \ge P_{\min }^{\nu }{\int }_{{\mathbb{T}}^d}\frac{1}{(Q+t1{)}^{\nu }}\mathrm{d}\mu \to \infty$$

$$c_{\mathbf{k}}$$

$$m_{\mathbf{k}}$$

$${\epsilon }_{\mathbf{k}}=\frac{\lambda \nu }{\nu -1}{\int }_{{\mathbb{T}}^d}{e}^{i⟨\mathbf{k},\mathbf{\theta }⟩}\frac{1}{{\hat{P}}^{\nu +1}}\mathrm{d}\mu ,\mathbf{k}\in {\Lambda }_0.$$

$$\delta J_{\nu ,\lambda }(P,Q;\delta P,\delta Q)=⟨\mathbf{c},\delta \mathbf{q}⟩-{\int }_{{\mathbb{T}}^d}\delta Q(P/Q{)}^{\nu }\mathrm{d}\mu$$

$$-⟨\mathbf{m},\delta \mathbf{p}⟩+\frac{\nu }{\nu -1}{\int }_{{\mathbb{T}}^d}\delta P\left[(P/Q{)}^{\nu -1}-\lambda /P^{\nu +1}\right]\mathrm{d}\mu$$

$$W(\mathbf{z})={\left[\frac{b(\mathbf{z})}{a(\mathbf{z})}\right]}^{\nu }={\left[\frac{b_0-b_1{z}_1^{-1}-b_2{z}_2^{-1}-b_3{z}_3^{-1}}{a_0-a_1{z}_1^{-1}-a_2{z}_2^{-1}-a_3{z}_3^{-1}}\right]}^{\nu }$$

$${\mathbf{a}}_1={\mathbf{a}}_2=\mathbf{a}$$

$${\tilde{\mathbf{b}}}_1=[1,-0.2,-0.3,-0.4],$$

$${\Lambda }_{+}:=\{(0,0,0),(1,0,0),(0,1,0),(0,0,1)\}.$$

Full text

I Introduction

Multidimensional stationary processes (or stationary random fields) represent a fundamental tool in many applications of signal and image processing. For instance, the data collected from an automotive radar system can be modeled by a multidimensional process of dimension $d=3$, see [1, 2]. In those applications we have to estimate the multidimensional spectral density of the process. This task can be addressed by means of a moment problem, more precisely, a convex optimization problem subject to moment constraints.

In the unidimensional case $(d=1)$ a wide range of spectral estimation paradigms based on moment problems have been proposed, see for instance [3, 4]: in this case the moments correspond to some covariance lags of the process and in the simplest setup the optimal spectrum maximizes the entropy rate. The appealing property of these paradigms is that the optimal spectrum is rational and thus leading to a finite-dimensional linear stochastic system (called “shaping filter” in the literature of signal processing) after spectral factorization.

In the case where the moments include both covariance lags and cepstral coefficients (i.e., logarithmic moments), then it is possible to characterize only an approximate solution to the moment problem: the spectrum maximizing the entropy rate matches the covariance lags and approximately the cepstral coefficients [5]. Such a solution is obtained by considering a regularized version of the dual optimization problem.

These paradigms have been extended also to the multidimensional case, see e.g., [6, 7]. Although spectral factorization is not always possible in the multidimensional setting, rationality still seems to be a key ingredient toward a finite-dimensional realization theory [8, 9]. The main issue is, however, that the solution of the moment problem is not necessarily a spectral density, but rather a spectral measure that may contain a singular part [10]. In particular, if the moments are both the covariance lags and the cepstral coefficients, the existence of an approximate rational solution is only guaranteed when the dimension is $d\leq 2$ [10].

In order to overcome this limitation on the dimension $d$, we have proposed a new moment problem, hereafter called $\nu$-moment problem, in which the entropy rate has been replaced by a more general definition of entropy, called $\nu$-entropy, whose derivation comes from the $\alpha$-divergence [11]. The definition of cepstral coefficients has been generalized accordingly. The $\nu$-moment problem is characterized by the integer parameter $\nu$. In [12, 13] we have shown that for any $d>2$ there exists an approximate rational solution to the $\nu$-moment problem under the regularity condition $\nu\geq d/2+1$. On the other hand, if the solution admits a spectral factorization, then the estimated process can be realized through a cascade of $\nu$ identical linear filters. Therefore, the larger $\nu$ is, the larger the complexity is in order to realize such a process. Accordingly, from a practical perspective one would take $\nu$ as small as possible.

The aim of the present letter is to show that the regularity condition can be weakened as $\nu\geq d/2$ and thus it is possible to estimate a $d$-dimensional process with $\nu=d/2$ whose realization, if admissible, is simpler than the one obtained using the theory in [12]. Such a result is achieved through a new regularization technique. Moreover, the technical proof of the existence of such a solution takes a different route from the one in [12].

The outline of this letter is as follows. In Section II we introduce the multidimensional $\nu$-moment problem and the regularization term for the dual problem, while in Section III we characterize the regularized dual problem. In Section IV we prove the existence of an approximate solution for the $\nu$-moment problem for $\nu\geq d/2$. Some numerical examples are provided in Section V. Finally, in Section VI we draw the conclusions.

II Problem formulation

Consider a $d$-dimensional real111We choose to present the theory for real random fields for simplicity. After suitable adaptation, complex random fields can be handled as well. stationary random field $\{y(\mathbf{t}):\mathbf{t}\in\mathbb{Z}^{d}\}$ with zero mean and a spectral density $\Phi(\boldsymbol{\theta})$. The latter is a nonnegative function on the $d$-dimensional frequency domain $\mathbb{T}^{d}:=(0,2\pi]^{d}$ and $\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{d})\in\mathbb{T}^{d}$ is a frequency vector. In the sequel, we shall use the notation $\Phi(e^{i\boldsymbol{\theta}})$ which is common in Complex Analysis. Here $e^{i\boldsymbol{\theta}}$ is a shorthand for the vector $(e^{i\theta_{1}},\dots,e^{i\theta_{d}})$ representing a point on the $d$-torus which is isomorphic to $\mathbb{T}^{d}$.

The $\nu$-entropy of the random field $y$ with the integer parameter $\nu>1$ is defined as

[TABLE]

where $\mathrm{d}\mu=\frac{1}{(2\pi)^{d}}\prod_{j=1}^{d}\mathrm{d}\theta_{j}$ is the normalized Lebesgue measure on $\mathbb{T}^{d}$. It is worth noting that $\mathbb{H}_{\nu}$ is a more general definition of entropy. Indeed, for the case $\nu=1$ (which is understood in a suitable limit sense), we obtain the usual entropy rate $\mathbb{H}_{1}(\Phi)=\int_{\mathbb{T}^{d}}\log\Phi\mathrm{d}\mu$ [12]. In this letter, we will consider the following multidimensional $\nu$-moment problem

[TABLE]

where $\langle{\mathbf{k}},\,{\boldsymbol{\theta}}\rangle=\sum_{j=1}^{d}k_{j}\theta_{j}$ is the standard inner product in $\mathbb{R}^{d}$ and the multidimensional exponential function is understood as $e^{i\langle{\mathbf{k}},\,{\boldsymbol{\theta}}\rangle}=\prod_{j=1}^{d}e^{ik_{j}\theta_{j}}$. In the moment constraints we have:

• •

$\mathbf{c}=\{c_{\mathbf{k}}\,:\,\mathbf{k}\in\Lambda\}$ which is a covariance multisequence of the random field $y$, namely $c_{\mathbf{k}}=\operatorname{{\mathbb{E}}}[y(\mathbf{t}+\mathbf{k})y(\mathbf{t})]$ where $\operatorname{{\mathbb{E}}}$ denotes the expectation operator. The finite index set $\Lambda\subset\mathbb{Z}^{d}$ contains $\mathbf{0}$ and is symmetric with respect to the origin, i.e., $\mathbf{k}\in\Lambda$ implies $-\mathbf{k}\in\Lambda$. It is well known that the covariances are the Fourier coefficients of the spectral density $\Phi$ which is the meaning of (2b).

• •

$\mathbf{m}=\{m_{\mathbf{k}}\,:\,\mathbf{k}\in\Lambda_{0}\}$ which is a multisequence of generalized cepstral coefficients, called $\nu$-cepstral coefficients, associated to the same random field $y$, see [12]. It is a generalization of the classic logarithmic moments used in [5, 6]. Notice that, the notion of $\nu$-cepstral coefficients is consistent with the objective functional (1) employed in (2). The index set $\Lambda_{0}:=\Lambda\backslash\{\mathbf{0}\}$ is such that $m_{\mathbf{0}}$ is excluded (for technical reasons).

The $\nu$-moment problem (2), which is also referred to as the primal optimization problem, can be used to perform spectral estimation. Assume that a dataset generated by $y$ has been collected. Then, it is possible to compute the sample estimates $\hat{c}_{\mathbf{k}}$ and $\hat{m}_{\mathbf{k}}$ of $c_{\mathbf{k}}$ and $m_{\mathbf{k}}$, respectively. Then, the estimate of $\Phi$ is the solution of (2) in which $c_{\mathbf{k}}$ and $m_{\mathbf{k}}$ are substituted by $\hat{c}_{\mathbf{k}}$ and $\hat{m}_{\mathbf{k}}$.

In [12], we derived the dual optimization problem and showed that, if the dual problem admits an interior-point solution, the optimal spectral density (primal variable) is a rational function of the form

[TABLE]

where $P$ and $Q$ are positive222By positive we mean that $P>0$ for any point on the $d$-torus. In contrast, if $P\geq 0$ and $P=0$ for certain points, we say that $P$ is nonnegative. trigonometric polynomials associated with the Lagrange multipliers. However, it is highly nontrivial to prove that the optimal dual variable lies in the interior of the feasible set. In order to overcome this difficulty, we introduced a regularization term in the dual objective function which depends solely on $P$:

[TABLE]

where $\lambda>0$ is a regularization parameter. Such a regularizer can of course be interpreted as a barrier function (for $P$) since under the regularity condition

[TABLE]

the regularizer takes an infinite value if $P$ has a zero on the $d$-torus. Hence, the optimal $P$ is forced to be an interior point, i.e., a positive polynomial. Obviously, one can always choose $\nu$ such that the condition (5) is met, and indeed in this case, we showed that the optimal $Q$ is also positive so that the optimal form (3) is true for the primal problem (2). The price to pay for the regularization is that the $\nu$-cepstral constraints (2c) are only approximately satisfied with an error that decreases as $\lambda\to 0_{+}$, i.e., the rational function (3) represents an approximate solution to (2). It is worth noting that a process with the spectrum (3) can be generated by the cascade of $\nu$ identical multidimensional filters (see Fig. 1) if $P/Q$ admits a spectral factor corresponding to a state space realization. Clearly, the larger $\nu$ is, the larger the complexity of the realization is. Accordingly, the key point from a practical perspective is to have the possibility to take $\nu$ as small as possible.

Having in mind this important requirement, in this letter we want to face the following problem.

Problem 1**.**

Take the barrier function in the dual problem as

[TABLE]

and show that the corresponding approximate solution of the primal problem exists under the weaker condition $\nu\geq d/2$.

Example 1**.**

Consider an automotive radar system installed in the red car of Fig. 2 that employs coherent linear frequency-modulated pulse trains signals (T) and uses a uniform linear array of receive antennas for the measurement (R). The target (green car in Fig. 2) is identified by the range $r$, the azimuth angle $\alpha$ and the relative velocity $v$. The problem of estimating the target parameters can be formulated as a multidimensional spectral estimation problem with $d=3$ [14, 1]. Under the aforementioned hypothesis on the existence of a “realizable” spectral factor, the estimated model can be realized by a cascade of at least $3$ identical filters using the theory developed in [12], while only $2$ using the theory that we will develop in this letter which makes more efficient the implementation of the model for simulation purposes.

III Dual problem

Let us recall the innocuous condition $\nu\geq 2$ which is assumed throughout this letter. In [12, Section 4], it has been shown that the dual function of the primal problem (2) is

[TABLE]

where:

• •

$\mathbf{q}=\{q_{\mathbf{k}}\,:\,\mathbf{k}\in\Lambda\}$ contains the real coefficients of the nonnegative trigonometric polynomial

[TABLE]

such that $q_{-\mathbf{k}}=q_{\mathbf{k}}$;

• •

similarly, $\mathbf{p}=\{p_{\mathbf{k}}\,:\,\mathbf{k}\in\Lambda_{0}\}$ contains the real coefficients of the nonnegative polynomial

[TABLE]

with $p_{-\mathbf{k}}=p_{\mathbf{k}}$ and $p_{\mathbf{0}}=1$ fixed;

• •

$\langle{\mathbf{q}},\,{\mathbf{c}}\rangle:=\sum_{\mathbf{k}\in\Lambda}q_{\mathbf{k}}c_{\mathbf{k}}$ inner product of two multisequences indexed in $\Lambda$, and $\langle{\mathbf{p}},\,{\mathbf{m}}\rangle$ is understood similarly with the index set replaced by $\Lambda_{0}$.

Since it is rather difficult to prove that the dual problem admits an interior-point solution, see [12], we consider the regularized dual function

[TABLE]

where the bivariate function is defined as

[TABLE]

with $x>0,\ y>0$ and the regularization parameter $\lambda>0$ is fixed. Now let us introduce the feasible sets

[TABLE]

so that $Q\in\mathfrak{P}_{+}$ and $P\in\mathfrak{P}_{+,o}$. The domain of definition of $J_{\nu,\lambda}$ can be extended to the boundary of the feasible set $\mathfrak{P}_{+,o}\times\mathfrak{P}_{+}$ by excluding the zero sets of $P$ and $Q$ from the domain of integration, which does not change the values of the integrals since the zero sets have zero Lebesgue measure. Moreover, $J_{\nu,\lambda}$ may take a value of $\infty$ at some boundary points, and hence it is understood as an extended real-valued function.

Lemma 1**.**

The function $g$ in (11) is strictly convex in the domain $x>0,\ y>0$ $($the first quadrant$)$.

Proof.

We shall prove the claim via the derivative test. After some straightforward computations, we arrive at

[TABLE]

where $\nu(\nu-1)>0$ is a positive integer. It is readily observed that every diagonal element in the Hessian of $g$ is positive, and the first matrix in (12b) is positive semidefinite. After checking the determinant of the Hessian, we conclude that $\nabla^{2}g(x,y)$ is positive definite in the first quadrant and the strict convexity follows. ∎

The next proposition is a direct consequence of Lemma 1.

Proposition 1**.**

The regularized dual function $J_{\nu,\lambda}$ is strictly convex in the closed set $\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$.

Proof.

We only need to show the strict convexity of the integral term in (10b) since the inner products are linear in $(P,Q)$. In the interior of the feasible set, namely $P\in\mathfrak{P}_{+,o}$ and $Q\in\mathfrak{P}_{+}$, the strict convexity of the integral term follows from that of $g$ (see Lemma 1), and this can be seen by a pointwise argument on the integrand (see e.g., Proposition 5.3 in [14]). The same reasoning works if $P$ or $Q$ is on the boundary of the respective feasible set because, once the zero sets of $P$ and $Q$ are excluded from the domain of integration, the function $g(P(e^{i\boldsymbol{\theta}}),Q(e^{i\boldsymbol{\theta}}))$ is well defined and the proof holds verbatim. ∎

We conclude that the regularized dual optimization can be formulated as

[TABLE]

IV A unique interior-point solution under the condition $\nu\geq{d}/{2}$

In this section we shall see how the condition $\nu\geq d/2$ guarantees that the regularized dual problem (13) has an interior-point solution. By Proposition 1, we know that a minimizer of $J_{\nu,\lambda}$ in $\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$ is unique provided that it exists. We shall first establish such existence. Then, under the regularity condition $\nu\geq d/2$ which is weaker than (5), we aim to exclude the possibility that a minimum of $J_{\nu,\lambda}$ may fall on the boundary of the feasible set, which is a weaker version of Lemma 5.8 in [12], using a Byrnes–Gusev–Lindquist-type argument that first appeared in [3] and subsequently in e.g., [6].

IV-A Existence of a minimizer in $\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$

The existence of a solution to (13) can be shown via reasonings similar to the ones in [12]: such existence depends on the following feasibility assumption.

Assumption 1** (Feasibility).**

The given covariances $\{c_{\mathbf{k}}\}_{\mathbf{k}\in\Lambda}$ admit an integral representation

[TABLE]

where $\Phi_{0}$ is a nonnegative function on $\mathbb{T}^{d}$ and is positive on some open ball $B_{1}\subset\mathbb{T}^{d}$.

In order to prove our existence result, we need the following lemmas.

Lemma 2**.**

The unregularized dual function $J_{\nu}$ in (7) and the regularized version $J_{\nu,\lambda}$ in (10) are lower-semicontinuous on $\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$. In particular, they are both continuous on $\mathfrak{P}_{+,o}\times\mathfrak{P}_{+}$.

Proof.

The proof is similar to that of Lemma 5.5 in [12]. The only difference here is that the power of $P$ in the regularizer in (10a) is $\nu$ instead of $\nu-1$ in [12]. ∎

Lemma 3**.**

Suppose that Assumption 1 holds. If a sequence $\{(P_{j},Q_{j})\}_{j\geq 1}\subset\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$ is such that $\|(P_{j},Q_{j})\|\to\infty$ as $j\to\infty$, then $J_{\nu,\lambda}(P_{j},Q_{j})\to\infty$.

Proof.

See the proof of Lemma 5.6 in [12] which uses Assumption 1. ∎

Proposition 2**.**

Under Assumption 1, the regularized dual optimization problem (13) admits a solution.

Proof.

Take a sufficiently large real number $\beta$, so that the sublevel set of the regularized dual function

[TABLE]

is not empty. Then Lemma 2 implies that the sublevel set is closed; Lemma 3 implies that the sublevel set is bounded. Obviously the polynomial pair $(P,Q)$, parametrized by their coefficients, belongs to a finite-dimensional vector space. It follows that the sublevel set is compact. Given the lower-semicontinuity of the objective function $J_{\nu,\lambda}$ (see Lemma 2), a minimizer exists in $J_{\nu,\lambda}^{-1}(-\infty,\beta]$ by the extreme value theorem of Weierstrass. ∎

IV-B Non-optimality of boundary points given $\nu\geq d/2$

In this subsection we prove that the optimal solution of (13) cannot belong to the boundary of $\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$ using arguments which conceptually differ from the ones used for the case $\nu\geq d/2+1$ in [12]. Notice first that if $P\in\partial\mathfrak{P}_{+,o}$ where $\partial$ denotes the boundary of a set, then the regularization term (6) employed in (10a) takes a value of $\infty$ under the condition $\nu\geq d/2$, see Proposition A.4 in [14]. The other term, namely the unregularized function $J_{\nu}(P,Q)$ whose expression is given in (7), is bounded from below under Assumption 1, see the proof of [12, Lemma 5.6]. Therefore, in this case we have $J_{\nu,\lambda}(P,Q)=\infty$ which is certainly not a minimum. Consequently, an optimal $(P,Q)$ must have $P\in\mathfrak{P}_{+,o}$.

Next, we work on the case of $(P,Q)\in\mathfrak{P}_{+,o}\times\partial\mathfrak{P}_{+}$. It is still possible that the integral term in (7) diverges so that $J_{\nu,\lambda}(P,Q)=\infty$, and such a point is obviously not a minimizer. Therefore, we only need to consider points $(P,Q)\in\mathfrak{P}_{+,o}\times\partial\mathfrak{P}_{+}$ such that the function value $J_{\nu,\lambda}(P,Q)$ is finite. Construct the real-valued function

[TABLE]

defined for $t\geq 0$ where $\mathbf{1}$ denotes the constant polynomial taking value one. One can show without difficulty that $f$ “inherits” the strict convexity from $J_{\nu,\lambda}$ (see Proposition 1). It is worth noting that the continuity of $f(t)$ can be extended to $t=0_{+}$ which corresponds to the boundary point $(P,Q)$, in contrast with the general result that $J_{\nu,\lambda}$ is only lower-semicontinuous on the boundary of the feasible set (see Lemma 2). These properties are established in the next statement.

Lemma 4**.**

The function $f(t)$ in (16), where $(P,Q)$ is an arbitrarily fixed point in $\mathfrak{P}_{+,o}\times\partial\mathfrak{P}_{+}$ such that $J_{\nu,\lambda}(P,Q)<\infty$, is strictly convex and continuous in $[0,\infty)$.

Proof.

The strict convexity of $f$ follows directly from Proposition 1, and it remains to show the continuity. For $t>0$, the argument $(P,Q+t\mathbf{1})$ belongs to the interior $\mathfrak{P}_{+,o}\times\mathfrak{P}_{+}$, and the continuity of $f$ follows from that of $J_{\nu,\lambda}$ (see Lemma 2). Hence, we are only concerned with the right continuity of $f$ at $t=0$. By the lower-semicontinuity of $J_{\nu,\lambda}$ at $(P,Q)$ (Lemma 2 again), for any $\varepsilon>0$, there exists $\delta>0$ such that whenever $(P_{1},Q_{1})\in\overline{\mathfrak{P}}_{+,o}\times\overline{\mathfrak{P}}_{+}$ satisfies $\|(P_{1},Q_{1})-(P,Q)\|<\delta$, we have $J_{\nu,\lambda}(P_{1},Q_{1})>J_{\nu,\lambda}(P,Q)-\varepsilon$. We can always choose $t$ sufficiently small to make the argument on the right-hand side of (16) sufficiently close to $(P,Q)$, which leads to the inequality $f(t)>f(0)-\varepsilon$. At the same time, by strict convexity we have $f(t)<tf(1)+(1-t)f(0)=f(0)+t[f(1)-f(0)]$. The last term, namely $t[f(1)-f(0)]$ can be made smaller than $\varepsilon$ (in absolute value) for $t$ sufficiently small. Therefore, we reach the inequality $|f(t)-f(0)|<\varepsilon$ which proves the continuity at $t=0_{+}$ since $\varepsilon>0$ is arbitrarily chosen. ∎

Proposition 3**.**

Fix any point $(P,Q)\in\mathfrak{P}_{+,o}\times\partial\mathfrak{P}_{+}$ such that $J_{\nu,\lambda}(P,Q)<\infty$. Then the function $f(t)$ in (16) has a derivative $f^{\prime}(t)\to-\infty$ as $t\to 0_{+}$. Therefore, $(P,Q)$ cannot be a minimizer of $J_{\nu,\lambda}$.

Proof.

After some computations, we have for $t>0$ that

[TABLE]

Since $P\in\mathfrak{P}_{+,o}$, we know $P_{\min}:=\min_{\boldsymbol{\theta}\in\mathbb{T}^{d}}P(e^{i\boldsymbol{\theta}})>0$. Then the following relation

[TABLE]

holds as $t\to 0_{+}$ by Lebesgue’s monotone convergence theorem [15, p. 21] and Proposition A.4 in [14]. This shows that $f^{\prime}(t)$ in (17) tends to $-\infty$ as $t$ goes to zero from the right. Taking Lemma 4 into account, we conclude that [math] is not a local minimizer of $f(t)$. Consequently, any boundary point $(P,Q)\in\mathfrak{P}_{+,o}\times\partial\mathfrak{P}_{+}$ is not a minimizer of $J_{\nu,\lambda}$, because taking an arbitrarily small step along the direction $(\mathbf{0},\mathbf{1})$, which points towards the interior $\mathfrak{P}_{+,o}\times\mathfrak{P}_{+}$, will result in a decrease of the objective function value. ∎

We summarize what we have got so far in the next theorem.

Theorem 1**.**

Under Assumption 1 and the condition $\nu>d/2$, the optimization problem (13) admits a unique interior-point solution $(\hat{P},\hat{Q})\in\mathfrak{P}_{+,o}\times\mathfrak{P}_{+}$ such that

[TABLE]

In plain words, the spectral density $\hat{\Phi}_{\nu}=(\hat{P}/\hat{Q})^{\nu}$ achieves covariance matching and approximate $\nu$-cepstral matching with errors

[TABLE]

Proof.

The existence of a solution is guaranteed by Proposition 2 and the uniqueness by Proposition 1. Moreover, given the reasoning at the beginning of this subsection and Proposition 3, the optimal $(\hat{P},\hat{Q})$ must be an interior point, i.e., both polynomials are positive. As a consequence, it must satisfy the stationary-point equation $\nabla J_{\nu,\lambda}(P,Q)=\mathbf{0}$, which is equivalent to the conditions in (19). Indeed, this point can be seen by setting the first differential of the regularized dual function

[TABLE]

equal to zero for any direction $(\delta P,\delta Q)$. ∎

Remark 1*.*

It is important to observe that the continuous dependence of the solution $(\hat{P},\hat{Q})$ to (13) on the covariance and $\nu$-cepstral data $(\mathbf{c},\mathbf{m})$ can be established similarly to [12, Sec. 6], so that the optimization problem (13) is in fact well-posed in the sense of Hadamard.

V Numerical simulations

In this section we present some numerical experiments in which the problem (2) is used to reconstruct the spectrum of a $3$-d stationary random field $y(t_{1},t_{2},t_{3})$ starting from a finite set of its covariance lags and $\nu$-cepstral coefficients. We assume that the underlying process $y$ is described as the output of a cascade linear shaping filter $\tilde{W}(z_{1},z_{2},z_{3})=W^{\nu}(z_{1},z_{2},z_{3})$ with a white noise input $e(t_{1},t_{2},t_{3})$, see Fig. 1 with $\nu=2$ and $d=3$. We shall also assume that the transfer function $W$ has the structure consistent with our optimal form (3) for the spectrum. More precisely, we consider the class of rational models such that the numerator and denominator polynomials both have degree one:

[TABLE]

where $\mathbf{z}$ represents $(z_{1},z_{2},z_{3})$ for short. Obviously, the polynomials $a(\mathbf{z})$ and $b(\mathbf{z})$ are described by the respective vectors $\mathbf{a}=[a_{0},\dots,a_{3}]$ and $\mathbf{b}=[b_{0},\dots,b_{3}]$ of coefficients. If the white noise input $e$ has unit variance, then the spectral density of the output process $y$ is $(P/Q)^{\nu}$ where $P(e^{i\boldsymbol{\theta}})=|b(e^{i\boldsymbol{\theta}})|^{2}$ and $Q(e^{i\boldsymbol{\theta}})=|a(e^{i\boldsymbol{\theta}})|^{2}$ are two nonnegative polynomials in $\overline{\mathfrak{P}}_{+}$. Throughout this section, we set the parameter $\nu=2$ to meet the condition $\nu\geq{d}/{2}=3/2$. In such scenario it is worth stressing that the theory developed in [12] does not work. Indeed that theory requires $\nu\geq 3$ to guarantee the existence of an approximate solution to the primal problem (2).

In what follows we consider two models of the form (21). We take two sets of real parameters $(\mathbf{a}_{j},\mathbf{b}_{j})$ with $j=1,2$ such that

[TABLE]

and $\mathbf{b}_{j}=\tilde{\mathbf{b}}_{j}/\|\tilde{\mathbf{b}}_{j}\|$. The last operation of normalization gives $\|\mathbf{b}\|^{2}=1$ which is equivalent to $p_{\mathbf{0}}=1$ for the numerator polynomial $P$. The first model, hereafter called “zeroless model”, corresponds to the process with spectrum $\Phi_{1}=(P_{1}/Q)^{2}$ where the polynomials $Q=|a|^{2}$ and $P_{1}=|b_{1}|^{2}$ are positive on $\mathbb{T}^{3}$; the second model, hereafter called “model with a spectral zero”, corresponds to the process with spectrum $\Phi_{2}=(P_{2}/Q)^{2}$ where $P_{2}=|b_{2}|^{2}$ (and also the spectrum) has a zero at the frequency vector $(\pi,\pi,\pi)$. The index set $\Lambda$ is identified as $\Lambda=\Lambda_{+}-\Lambda_{+}$ with

[TABLE]

Here the set difference is understood as $A-B:=\{x-y:x\in A,\ y\in B\}$. Since $q_{-\mathbf{k}}=q_{\mathbf{k}}$ and $p_{-\mathbf{k}}=p_{\mathbf{k}}$, the total number of variables is $13$.

By means of the discrete formulation described in Section 7 of [12], we consider the following procedure to test the ability to reconstruct the spectra through the solution of (2) for the previous two models:

Fix $N=20$ and discretize $\mathbb{T}^{3}$ into $N^{3}$ regular grid points by gridding the interval $[0,2\pi]$ into $N$ equidistant points in each dimension. 2. 2.

Compute the “true” covariances $\{c_{\mathbf{k}}:\mathbf{k}\in\Lambda\}$ and the $\nu$-cepstral coefficients $\{m_{\mathbf{k}}:\mathbf{k}\in\Lambda_{0}\}$ of $\Phi=|W|^{2}$ via (2b) and (2c) where $\mathrm{d}\mu$ is replaced by a discrete measure with equal mass $1/N^{3}$ on the grid points. 3. 3.

Solve the discrete version of the regularized dual problem (13) using $(\mathbf{c},\mathbf{m})$ computed above and $\lambda>0$ chosen sufficiently small. 4. 4.

Let $(\hat{\mathbf{p}},\hat{\mathbf{q}})$ be the optimal solution to (13) and $(\mathbf{p},\mathbf{q})$ be the polynomial coefficients corresponding to the true spectrum $\Phi=|W|^{2}$. Finally, evaluate the reconstruction error $\|(\hat{\mathbf{p}},\hat{\mathbf{q}})-(\mathbf{p},\mathbf{q})\|$.

In Step 3, the optimization problem is solved using Newton’s method. Some computational details can be found in [12], and suitable modifications on the gradient and Hessian of the objective function are needed for the current problem (13).

The left panel of Fig. 3 shows the reconstruction errors for the two models with different values of the regularization parameter $\lambda=10^{-n}$, $n=0,2,4,6,8,10$. It is readily observed in both cases that the errors decrease monotonically as $\lambda\to 0$. In view of Remark 1, therefore, if we have a dataset generated by a $d$-dimensional process and consider the approximate solution to the problem (2) with $\nu=\lceil{d/2}\rceil$ where $\mathbf{c}$ and $\mathbf{m}$ are replaced by their sample estimators (computed from the data), then the resulting spectral estimator is characterized by a small estimation error provided that $\lambda$ is chosen sufficiently small and the underlying process has a spectrum of the form (3) with $\nu=\lceil{d/2}\rceil$ which agrees with our parameter specification.

In the right panel of Fig. 3, we compare the spectral density of the model with a spectral zero and the reconstructed spectra for values of $\lambda=10^{-10},10^{-8},10^{-6},10^{-4}$ (corresponding to the orange line in the left panel) along a cross section $[\,\cdot\;11\;11\,]$ of the regular grid for $\mathbb{T}^{3}$, where the true spectral zero is located at the index333In Matlab, the array indices start with $1$. $[11,11,11]$. Notice that the other two cross sections $[\,11\;\cdot\;11\,]$ and $[\,11\;11\;\cdot\,]$ of the spectral densities are not shown because they are visually similar to this figure. We conclude that as $\lambda$ goes to zero, the true nonnegative spectrum is well approximated by a positive spectral density with smaller and smaller errors, which is consistent with the left panel.

VI Conclusions

We have considered a multidimensional $\nu$-moment problem which searches a spectral density maximizing the $\nu$-entropy and matching a finite set of covariance lags and $\nu$-cepstral coefficients. A previous work showed that it is possible to guarantee the existence of a rational approximate solution only in the case $\nu\geq d/2+1$ where $d$ is the dimension of the random field described by the multidimensional spectrum. Motivated by the fact that $\nu$ should be chosen as small as possible, we have proposed a different regularization technique which ensures the existence of a rational approximate solution to the $\nu$-moment problem under the weaker condition $\nu\geq d/2$.