Multidimensional Rational Covariance Extension with Approximate Covariance Matching

TL;DR

This paper extends the multidimensional rational covariance extension problem to handle approximate covariance matching, ensuring well-posedness and demonstrating applications in spectral estimation and texture generation.

Contribution

It introduces a framework for approximate covariance matching in RCEP, connecting soft and hard constraints, and proves their well-posedness with practical examples.

Findings

The problems are well-posed under approximate covariance matching.

Soft and hard constraint problems are connected via a homeomorphism.

Applications include spectral estimation and texture synthesis.

Abstract

In our companion paper "Multidimensional rational covariance extension with applications to spectral estimation and image compression" we discussed the multidimensional rational covariance extension problem (RCEP), which has important applications in image processing, and spectral estimation in radar, sonar, and medical imaging. This is an inverse problem where a power spectrum with a rational absolutely continuous part is reconstructed from a finite set of moments. However, in most applications these moments are determined from observed data and are therefore only approximate, and RCEP may not have a solution. In this paper we extend the results to handle approximate covariance matching. We consider two problems, one with a soft constraint and the other one with a hard constraint, and show that they are connected via a homeomorphism. We also demonstrate that the problems are well-posed and illustrate the theory by examples in spectral estimation and texture generation.

Tables

Table 1. Table 1: Norm differences ‖ r ^ − c true ‖ 2 subscript norm ^ 𝑟 subscript 𝑐 true 2 \|\hat{r}-c_{\rm{true}}\|_{2} for different solutions in the first simulation setup.

	Mean	Std.
Biased, exact matching	$3.2374$	$1.7944$
Biased, approximate matching, ME-solution	$3.7886$	$1.3274$
Biased, approximate matching, using true $P$	$3.8152$	$1.6509$
Unbiased, approximate matching, ME-solution	$3.2575$	$1.4721$
Unbiased, approximate matching, using true $P$	$3.2811$	$1.7787$

Table 2. Table 2: Norm differences ‖ r ^ − c true ‖ 2 subscript norm ^ 𝑟 subscript 𝑐 true 2 \|\hat{r}-c_{\rm{true}}\|_{2} for different solutions in the second setup, where all unbiased estimate are outside ℭ ¯ + subscript ¯ ℭ \bar{\mathfrak{C}}_{+} .

	Mean	Std.
Biased, exact matching	$2.9245$	$2.2528$
Biased, approximate matching, ME-solution	$1.9087$	$1.1324$
Biased, approximate matching, using true $P$	$1.8532$	$1.1904$
Unbiased, approximate matching, ME-solution	$1.5018$	$0.6601$
Unbiased, approximate matching, using true $P$	$1.4451$	$0.7296$

Equations

$$c_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i(\mathbf{k},\mathbf{\theta })}d\mu (\mathbf{\theta })\text{for all k∈Λ},$$

$$c:=[c_{\mathbf{k}}\mid \mathbf{k}:=(k_1,\dots ,k_d)\in \Lambda ]$$

$$d\mu (\mathbf{\theta })=\Phi (e^{i\mathbf{\theta }})dm(\mathbf{\theta })+d\nu (\mathbf{\theta }),$$

$$dm(\mathbf{\theta }):=(1/2\pi {)}^d\prod _{j=1}^{d}d{\theta }_j$$

$$\Phi (e^{i\mathbf{\theta }})=\frac{P(e^{i\mathbf{\theta }})}{Q(e^{i\mathbf{\theta }})},p,q\in {\bar{\mathfrak{P}}}_{+}\backslash \{0\},$$

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

$${\mathfrak{C}}_{+}:=\left\{c\mid ⟨c,p⟩>0,\text{for all p∈Pˉ+∖{0}}\right\},$$

$$r_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i(\mathbf{k},\mathbf{\theta })}d\mu (\mathbf{\theta }),\mathbf{k}\in \Lambda ,$$

$$\lambda \mathbb{D}(Pdm,d\mu )+\frac{1}{2}\Vert r-c{\Vert }_2^2$$

$$\mathbb{D}(Pdm,d\mu ):={\int }_{{\mathbb{T}}^d}\left(P\log \frac{P}{\Phi }dm+d\mu -Pdm\right)$$

$$\underset{d\mu \ge 0,r}{\min }$$

$$\mathrm{subjectto}$$

$$\Vert r-c{\Vert }^2\le \lambda .$$

$$\underset{d\mu \ge 0,r}{\min }$$

$$\mathrm{subjectto}$$

$$\Vert r-c{\Vert }_{W^{-1}}^2\le 1,$$

$$⟨q,p⟩={\int }_{{\mathbb{T}}^d}Q(e^{i\mathbf{\theta }})\overline{P(e^{i\mathbf{\theta }})}dm(\mathbf{\theta }),$$

$$d\mu =P(e^{i\mathbf{\theta }})dm(\mathbf{\theta }),p\in {\bar{\mathfrak{P}}}_{+}\backslash \{0\}.$$

$$C(e^{i\mathbf{\theta }}):=\sum _{\mathbf{k}\in \Lambda }{c}_{\mathbf{k}}{e}^{-i(\mathbf{k},\mathbf{\theta })}.$$

$$\underset{d\mu \ge 0}{\min }$$

$$\mathrm{subjectto}$$

$$\underset{q\in {\bar{\mathfrak{P}}}_{+}}{\min }$$

$$\mathrm{supp}(d\hat{\nu })\subseteq \{\mathbf{\theta }\in {\mathbb{T}}^d\mid \hat{Q}(e^{i\mathbf{\theta }})=0\}$$

$$c_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i(\mathbf{k},\mathbf{\theta })}\left(\frac{P}{\hat{Q}}dm+d\hat{\nu }\right),\mathbf{k}\in \Lambda ,$$

$${\hat{c}}_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i(\mathbf{k},\mathbf{\theta })}d\hat{\nu },\mathbf{k}\in \Lambda .$$

$$d\hat{\mu }(\mathbf{\theta })=\frac{P(e^{i\mathbf{\theta }})}{\hat{Q}(e^{i\mathbf{\theta }})}dm(\mathbf{\theta })+d\hat{\nu }(\mathbf{\theta })$$

$$\hat{q}\in {\bar{\mathfrak{P}}}_{+},\hat{c}\in \partial {\mathfrak{C}}_{+},⟨\hat{c},\hat{q}⟩=0$$

$$c_{\mathbf{k}}={\int }_{{\mathbb{T}}^d}{e}^{i(\mathbf{k},\mathbf{\theta })}\frac{P}{\hat{Q}}dm+{\hat{c}}_{\mathbf{k}},\mathbf{k}\in \Lambda ,$$

$$\mathbb{D}(Pdm,d\mu )$$

$$\ge {\int }_{{\mathbb{T}}^d}\left(P(1-\frac{\Phi }{P})dm+\Phi dm-Pdm\right)+{\int }_{{\mathbb{T}}^d}d\nu \ge 0,$$

$$\mathbb{D}(Pdm,d\mu )+\frac{1}{2}\Vert r-c{\Vert }_{W^{-1}}^2$$

$$\mathcal{L}(\Phi ,d\nu ,r,q)=$$

$$=$$

$$\mathbb{D}(Pdm,d\mu )={\int }_{{\mathbb{T}}^d}P(\log P-1)dm-{\int }_{{\mathbb{T}}^d}P\log \Phi dm+r_0,$$

$$\mathcal{L}(\Phi ,d\nu ,r,q)=$$

Full text

\slugger

siconxxxxxxxx–x

1 Introduction

Trigonometric moment problems are ubiquitous in systems and control, such as spectral estimation, signal processing, system identification, image processing and remote sensing [5, 20, 59]. In the (truncated) multidimensional trigonometric moment problem we seek a nonnegative measure $d\mu$ on $\mathbb{T}^{d}$ satisfying the moment equation

[TABLE]

where $\mathbb{T}:=(-\pi,\pi]$, ${\boldsymbol{\theta}}:=(\theta_{1},\ldots,\theta_{d})\in\mathbb{T}^{d}$, and $({\boldsymbol{k}},{\boldsymbol{\theta}}):=\sum_{j=1}^{d}k_{j}\theta_{j}$ is the scalar product in $\mathbb{R}^{d}$. Here $\Lambda\subset\mathbb{Z}^{d}$ is a finite index set satisfying $0\in\Lambda$ and $-\Lambda=\Lambda$. A necessary condition for (1.1) to have a solution is that the sequence

[TABLE]

satisfy the symmetry condition $c_{-{\boldsymbol{k}}}=\bar{c}_{\boldsymbol{k}}$. The space of sequences (1.2) with this symmerty will be denoted $\mathfrak{C}$ and will be represented by vectors $c$, formed by ordering the coefficient in some prescribed manner, e.g., lexiographical. Note that $\mathfrak{C}$ is isomorphic to $\mathbb{R}^{|\Lambda|}$, where $|\Lambda|$ is the cardinality of $\Lambda$. However, as we shall see below, not all $c\in\mathfrak{C}$ are bona fide moments for nonnegative measures $d\mu$.

In many of the applications mentioned above there is a natural complexity constraint prescribed by design specifications. In the context of finite-dimensional systems these constraints often arise in the requirement that transfer functions be rational. This leads to the rational covariance extension problem, which has been studied in various degrees of generality in [25, 26, 36, 53, 54] and can be posed as follows.

Define $e^{i{\boldsymbol{\theta}}}:=(e^{i\theta_{1}},\ldots,e^{i\theta_{d}})$ and let

[TABLE]

where $\bar{\mathfrak{P}}_{+}$ is the closure of the convex cone $\mathfrak{P}_{+}$ of the coefficients $p\in\mathfrak{C}$ corresponding to trigonometric polynomials

[TABLE]

that are positive for all ${\boldsymbol{\theta}}\in\mathbb{T}^{d}$.

The reason for referring to this problem as a rational covariance extension problem is that the numbers (1.2) correspond to covariances $c_{\boldsymbol{k}}:=\operatorname*{\mathbb{E}}\{y({\bf t}+{\boldsymbol{k}})\overline{y({\bf t})}\}$ of a discrete-time, zero-mean, and homogeneous111Homogeneity generalizes stationarity in the case $d=1$. stochastic process $\{y({\bf t});\,{\bf t}\in\mathbb{Z}^{d}\}$. The corresponding power spectrum, representing the energy distribution across frequencies, is defined as the nonnegative measure $d\mu$ on $\mathbb{T}^{d}$ whose Fourier coefficients are the covariances (1.2). A scalar version of this problem ($d=1$) was first posed by Kalman [34] and has been extensively studied and solved in the literature [24, 12, 6, 21, 48, 13, 41, 7, 61, 47]. It has been generalized to more general scalar moment problems [8, 27, 9] and to the multidimensional setting [26, 25, 54, 53, 36]. Also worth mentioning here is work by Lang and McClellan [39, 40, 45, 46, 38, 37] considering the multidimensional maximum entropy problem, which hence has certain overlap with the above literature.

The multidimensional rational covariance extension problem posed above has a solution if and only if $c\in\mathfrak{C}_{+}$, where $\mathfrak{C}_{+}$ is the open convex cone

[TABLE]

where $\langle c,p\rangle:=\sum_{{\boldsymbol{k}}\in\Lambda}c_{\boldsymbol{k}}\bar{p}_{\boldsymbol{k}}$ is the inner product in $\mathfrak{C}$ (Theorem 2.5). However, the covariances $[c_{\boldsymbol{k}}\mid{\boldsymbol{k}}:=(k_{1},\ldots,k_{d})\in\Lambda]$ are generally determined from statistical data. Therefore the condition $c\in\mathfrak{C}_{+}$ may not be satisfied, and testing this condition is difficult in the multidimensional case. Therefore we may want to find a positive measure $d\mu$ and a corresponding $r\in\mathfrak{C}_{+}$, namely

[TABLE]

so that $r$ is close to $c$ in some norm, e.g., the Euclidian norm $\|\centerdot\|_{2}$. This is an ill-posed inverse problem which in general has an infinite number of solutions $d\mu$. As we already mentioned, we are interested in rational solutions (1.3), and to obtain such solutions we use regularization as in [54]. Hence, we seek a $d\mu$ that minimizes

[TABLE]

subject to (1.5), where $\lambda\in\mathbb{R}$ is a regularization parameter and

[TABLE]

is the nomalized Kullback-Leibler divergence [33, ch. 4] [15, 61]. As will be explained in Section 2, $\mathbb{D}(Pdm,d\mu)$ is always nonnegative and has the property $\mathbb{D}(Pdm,Pdm)=0$.

In this paper we shall consider a more general problem in the spirit of [22]. To this end, for any Hermitian, positive definite matrix $M$, we define the weighted vector norm $\|x\|_{M}:=(x^{*}Mx)^{1/2}$ and consider the problem

[TABLE]

which is the same as the problem above with $W=\lambda I$. We shall refer to $W$ as the weight matrix.

Using the same principle as in [57], we shall also consider the problem to minimize $\mathbb{D}(Pdm,d\mu)$ subject to (1.5) and the hard constraint

[TABLE]

Since (1.5) are bona fide moments and hence $r\in\mathfrak{C}_{+}$, while $c\not\in\mathfrak{C}_{+}$ in general, this problem will not have a solution if the distance from $c$ to $\mathfrak{C}_{+}$ is greater than $\sqrt{\lambda}$. Hence the choice of $\lambda$ must be made with some care. Analogously with the rational covariance extension with soft constraints in (1.7), we shall consider the more general problem

[TABLE]

to which we shall refer as the rational covariance extension problem with hard constraints. Again this problem reduces to the simpler problem by setting $W=\lambda I$.

As we shall see, the soft-constrained problem (1.7) always has a solution, while the hard-constrained problem (1.9) may fail to have a solution for some weight matrices $W$. However, in Section 7 we show that the two problems are in fact equivalent in the sense that whenever (1.9) has a solution there is a corresponding $W$ in (1.7) that gives the same solution, and any solution of (1.7) can also be obtained from (1.9) by a suitable choice of $W$. The reason for considering both formulations is that one formulation might be more suitable than the other for the particular application at hand. For example, an absolute error estimate for the covariances is more naturally incorporated in the formulation with hard constraints. A possible choice of the weight matrix $W$ in either formulation would be the covariance matrix of the estimated moments, as suggested in [22]. This corresponds to the Mahalanobis distance and could be a natural way to incorporate uncertainty of the covariance estimates in the spectral estimation procedure.

Previous work in this direction can be found in [58, 22, 10, 57, 35], where [58, 35, 10] consider the problem of selecting an appropriate covariances sequence to match in a given confidence region. The two approximation problems considered here are similar to the ones considered in [57] and [22]. (For more details, also see [3, Ch. B].)

We begin in Section 2 by reviewing the regular multidimensional rational covariance extension problem for exact covariance matching in a broader perspective. In Section 3 we present our main results on approximate rational covariance extension with soft constraints, and in Section 4 we show that the dual solution is well-posed. In Section 5 we investigate conditions under which there are solutions without a singular part. The approximate rational covariance extension with hard constraints is considered in Section 6, and in Section 7 we establish a homeomorhism between the weight matrices in the two problems, showing that the problems are actually equivalent when solutions exist. We also show that under certain conditions the homeomorphism can be extended to hold between all sets of parameters, allowing us to carry over results from the soft-constrained setting to the hard-constrained one. In Section 8 we discuss the properties of various covariance estimators, in Section 9 we give a 2D example from spectral estimation, and in Section 10 we apply our theory to system identification and texture reconstruction. Some of the results of this paper were announced in [55] without proofs.

2 Rational covariance extension with exact matching

The trigonometric moment problem to determine a positive measure $d\mu$ satisfying (1.1) is an inverse problem that has a solution if and only if $c\in\bar{\mathfrak{C}}_{+}$ [36, Theorem 2.3], where $\bar{\mathfrak{C}}_{+}$ is the closure of $\mathfrak{C}_{+}$, and then in general it has infinitely many solutions. However, the nature of possible rational solutions (1.3) will depend on the location of $c$ in $\bar{\mathfrak{C}}_{+}$. To clarify this point we need the following lemma.

Lemma 2.1**.**

$\bar{\mathfrak{P}}_{+}\setminus\{0\}\subset\mathfrak{C}_{+}$*. *

Proof 2.2**.**

Obviously the inner product $\langle q,p\rangle:=\sum_{{\boldsymbol{k}}\in\Lambda}q_{\boldsymbol{k}}\bar{p}_{\boldsymbol{k}}$ can be expressed in the integral form

[TABLE]

*and therefore $\langle q,p\rangle>0$ for all $q,p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$, as $P$ and $Q$ can have zeros only on sets of measure zero. Hence the statement of the lemma follows. *

Therefore, under certain particular conditions, the multidimensional rational covariance extension problem has a very simple solution with a polynomial spectral density, namely

[TABLE]

Proposition 2.3**.**

The multidimensional rational covariance extension problem has a unique polynomial solution (2.2) if and only if $c\in\bar{\mathfrak{P}}_{+}\backslash\{0\}$, namely $P=C$, where

[TABLE]

The proof of Proposition 2.3 is immediate by noting that any such $C$ is a bona fide spectral density and noting that $c_{\boldsymbol{k}}=\int_{\mathbb{T}^{d}}e^{i({\boldsymbol{k}},{\boldsymbol{\theta}})}C(e^{i{\boldsymbol{\theta}}})dm({\boldsymbol{\theta}})$.

As seen from the following result presented in [36, Section 6], the other extreme occurs for $c\in\partial\mathfrak{C}_{+}:=\bar{\mathfrak{C}}_{+}\setminus\mathfrak{C}_{+}$, when only singular solutions exist.

Proposition 2.4**.**

*For any $c\in\partial\mathfrak{C}_{+}$ there is a solution $d\mu$ of (1.1) with support in at most $|\Lambda|-1$ points. There is no solution with a absolutely continuous part $\Phi dm$. *

However, for any $c\in\mathfrak{C}_{+}$, there is a rational solution (1.3) parametrized by $p\in\bar{\mathfrak{P}}_{+}\backslash\{0\}$, as demonstrated in [54] by considering a primal-dual pair of convex optimization problems. In that paper the primal problem is a weighted maximum entropy problem, but as also noted in [54, Sec. 3.2], it is equivalent to

[TABLE]

where $\Phi dm$ is the absolutely continuous part of $d\mu$. This amounts to minimizing the (regular) Kullback-Leibler divergence between $Pdm$ and $d\mu$, subject to $d\mu$ matching the given data [27, 54]. In the present case of exact covariance matching, this problem is equivalent to minimizing (1.6) subject to (1.1), since $P$ is fixed and the total mass of $d\mu$ is determined by the [math]:th moment $c_{\boldsymbol{0}}=\int_{\mathbb{T}^{d}}d\mu$. Hence both $\int_{\mathbb{T}^{d}}d\mu$ and $\int_{\mathbb{T}^{d}}Pdm$ are constants in this case. Hence problem (1.7) is the natural extension of (2.3) for the case where the covariance sequence is not known exactly.

The primal problem (2.3) is a problem in infinite dimensions, but with a finite number of constraints. The dual to this problem will then have a finite number of variables but an infinite number of constraints and is given by

[TABLE]

In particular, Theorem 2.1 in [54], based on corresponding analysis in [36], reads as follows.

Theorem 2.5**.**

Problem (2.3) has a solution if and only if $c\in\mathfrak{C}_{+}$. For every $c\in\mathfrak{C}_{+}$ and $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$ the functional in (2.4) is strictly convex and has a unique minimizer $\hat{q}\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$. Moreover, there exists a unique $\hat{c}\in\partial\mathfrak{C}_{+}$ and a (not necessarily unique) nonnegative singular measure $d\hat{\nu}$ with support

[TABLE]

such that

[TABLE]

For any such $d\hat{\nu}$, the measure

[TABLE]

*is an optimal solution to the problem (2.3). Moreover, $d\hat{\nu}$ can be chosen with support in at most $|\Lambda|-1$ points, where $|\Lambda|$ is the cardinality of the index set $\Lambda$. *

If $c\in\partial\mathfrak{C}_{+}$, only a singular measure with finite support would match the moment condition (Proposition 2.4). In this case, the problem (2.3) makes no sense, since any feasible solution has infinite objective value.

In [36] we also derived the KKT conditions

[TABLE]

which are necessary and sufficient for optimality of the primal and dual problems.

Since (2.3) is an inverse problem, we are interested in how the solution depends on the parameters of the problem. From Propositions 7.3 and 7.4 in [36] we have the following result.

Proposition 2.6**.**

*Let $c$, $p$ and $\hat{q}$ be as in Theorem 2.5. Then the map $(c,p)\mapsto\hat{q}$ is continuous. *

To get a full description of well-posedness of the solution we would like to extend this continuity result to the map $(c,p)\mapsto(\hat{q},\hat{c})$. However, such a generalization is only possible under certain conditions. The following result is a consequence of Proposition 2.6 and [54, Corollary 2.3].

Proposition 2.7**.**

*Let $c$, $p$, $\hat{q}$ and $\hat{c}$ be as in Theorem 2.5. Then, for $d\leq 2$ and all $(c,p)\in\mathfrak{C}_{+}\times\mathfrak{P}_{+}$, the mapping $(c,p)\to(\hat{q},\hat{c})$ is continuous. *

Corollary 2.3 in [54] actually ensures that $\hat{c}=0$ for $d\leq 2$ and $p\in\mathfrak{P}_{+}$. However, in Section 4.2 we present a generalization of Proposition 2.7 to cases with $d\geq 3$, where then $\hat{c}$ may be nonzero. (The proof of this generalization can be found in [52].) Here we shall also consider an example where continuity fails when $p$ belongs to the boundary $\partial\mathfrak{P}_{+}:=\bar{\mathfrak{P}}_{+}\setminus\mathfrak{P}_{+}$, i.e., the corresponding nonnegative trigonometric polynomial $P(e^{i{\boldsymbol{\theta}}})$ is zero in at least one point.

3 Approximate covariance extension with soft constraints

To handle the case with noisy covariance data, when $c$ may not even belong to $\mathfrak{C}_{+}$, we relax the exact covariance matching constraint (1.1) in the primal problem (2.3) to obtain the problem (1.7). In this case it is natural to reformulate the objective function in (2.3) to include a term that also accounts for changes in the total mass of $d\mu$. Consequently, we have exchanged the objective function in (2.3) by the normalized Kullback-Leibler divergence (1.6) plus a term that ensures approximate data matching.

Using the normalized Kullback-Leibler divergence, as proposed in [33, ch. 4] [15, 61], is an advantage in the approximate covariance matching problem since this divergence is always nonnegative, precisely as is the case for probability densities. To see this, observe that, in view of the basic inequality $x-1\geq\log x$,

[TABLE]

since $d\nu$ is a nonnegative measure. Moreover, $\mathbb{D}(Pdm,Pdm)=0$, as can be seen by taking $d\mu=Pdm$ in (1.6).

The problem under consideration is to find a nonnegative measure $d\mu=\Phi dm+d\nu$ minimizing

[TABLE]

subject to (1.5). To derive the dual of this problem we consider the corresponding maximization problem and form the Lagrangian

[TABLE]

where $q:=[q_{\boldsymbol{k}}\mid{\boldsymbol{k}}:=(k_{1},\ldots,k_{d})\in\Lambda]$ are Lagrange multitipliers and $Q$ is the corresponding trigonometric polynomial (1.4). However,

[TABLE]

and therefore

[TABLE]

where $e:=[e_{\boldsymbol{k}}]_{{\boldsymbol{k}}\in\Lambda}$, $e_{\boldsymbol{0}}=1$ and $e_{{\boldsymbol{k}}}=0$ for ${\boldsymbol{k}}\in\Lambda\setminus\{\boldsymbol{0}\}$, and hence $r_{\boldsymbol{0}}=\langle r,e\rangle$.

In deriving the dual functional

[TABLE]

to be minimized, we only need to consider $q\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$, as $\varphi$ will take infinite values for $q\not\in\bar{\mathfrak{P}}_{+}$. In fact, following along the lines of [54, p. 1957], we note that, if $Q(e^{i{\boldsymbol{\theta}}_{0}})<0$, (3) will tend to infinity when $\nu({\boldsymbol{\theta}}_{0})\to\infty$. Moreover, since $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$, there is a neighborhood where $P(e^{i{\boldsymbol{\theta}}})>0$, letting $\Phi$ tend to infinity in this neighborhood, (3) will tend to infinity if $Q\equiv 0$. We also note that the nonnegative function $\Phi$ can only be zero on a set of measure zero; otherwise the first term in (3) will be $-\infty$.

The directional derivative222Formally, the Gateaux differential [44]. of the Lagrangian (3) in any feasible direction $\delta\Phi$, i.e., any direction $\delta\Phi$ such that $\Phi+\varepsilon\delta\Phi\geq 0$ for sufficiencly small $\varepsilon>0$, is easily seen to be

[TABLE]

In particular, the direction $\delta\Phi:=\Phi\,\text{sign}(P-Q\Phi)$ is feasible since $(1\pm\varepsilon)\Phi\geq 0$ for $0<\varepsilon<1$. Therefore, any maximizing $\Phi$ must satisfy $\int_{\mathbb{T}^{d}}|P-Q\Phi|dm\leq 0$ and hence (1.3b). Moreover, a maximizing choice of $d\nu$ will require that

[TABLE]

as this nonnegative term can be made zero by the simple choice $d\nu\equiv 0$, and consequently (2.5) must hold. Finally, the directional derivative

[TABLE]

is zero for all $\delta r\in\mathfrak{C}$ if

[TABLE]

Inserting this together with (1.3b) and (3.3) into (3) then yields the dual functional

[TABLE]

Consequently the dual of the (primal) optimization problem (1.7) is equivalent to

[TABLE]

Theorem 3.1**.**

For every $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$ the functional in (3.5) is strictly convex and has a unique minimizer $\hat{q}\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$. Moreover, there exists a unique $\hat{r}\in\mathfrak{C}_{+}$, a unique $\hat{c}\in\partial\mathfrak{C}_{+}$ and a (not necessarily unique) nonnegative singular measure $d\hat{\nu}$ with support

[TABLE]

such that

[TABLE]

and the measure

[TABLE]

*is an optimal solution to the primal problem (1.7). Moreover, $d\hat{\nu}$ can be chosen with support in at most $|\Lambda|-1$ points. *

Proof 3.2**.**

The objective functional $\mathbb{J}$ of the dual problem (3.5) can be written as the sum of two terms, namely

[TABLE]

where $\tilde{c}\in\mathfrak{C}_{+}$. The functional $\mathbb{J}_{1}$ is strictly convex (Theorem 2.5), and trivially the same holds for $\mathbb{J}_{2}$ since it is a positive definite quadratic form. Consequently, $\mathbb{J}=\mathbb{J}_{1}+\mathbb{J}_{2}$ is strictly convex, as claimed. Moreover, $\mathbb{J}_{1}$ is lower semicontinuous [54, Lemma 3.1] with compact sublevel sets $\mathbb{J}_{1}^{-1}(-\infty,\rho]$ [54, Lemma 3.2]. Likewise, $\mathbb{J}_{2}$ is continuous with compact sublevel sets. Therefore $\mathbb{J}$ is lower semicontinuous with compact sublevel sets and therefore has a minimum $\hat{q}$, which must be unique by strict convexity.

In view of (3.4), the optimal value of $r$ is given by

[TABLE]

and is hence unique. Since therefore the linear term $c+W(q-e)$ in the gradient of $\mathbb{J}$ takes the value $\hat{r}$ at the optimal point, the analysis in [54, sect. 3.1.5] applies with obvious modifications, showing that there is a $\hat{c}\in\bar{\mathfrak{C}}_{+}$, which then must be unique, such that

[TABLE]

Moreover, there is a discrete measure $d\hat{\nu}$ with support in at most $|\Lambda|-1$ points such that (3.7b) holds; see, e.g., [36, Proposition 2.4]. Then (3.7a) holds as well. In view of (3.3),

[TABLE]

and consequently $\hat{c}\in\partial\mathfrak{C}_{+}$, and the support of $d\hat{\nu}$ must satisfy (3.6).

Finally, let $r$ be given in terms of $d\mu$ by (1.5), and let $\mathbb{I}(d\mu)$ be the corresponding primal functional in (1.7). Then, for any such $d\mu$,

[TABLE]

*and hence $d\hat{\mu}$ is an optimal solution to the primal problem (1.7), as claimed. *

We collect the KKT conditions in the following corollary.

Corollary 3.3**.**

The conditions

[TABLE]

*are necessary and sufficient conditions for optimality of the dual pair (1.7) and (3.5) of optimization problems. *

4 On the well-posedness of the soft-constrained problem

In the previous sections we have shown that the primal and dual optimization problems are well-defined. Next we investigate the well-posedness of the primal problem as an inverse problem. Thus, we first establish continuity of the solutions $\hat{q}$ in terms of the parameters $W$, $c$ and $p$.

4.1 Continuity of $\hat{q}$ with respect to $c$, $p$ and $W$

We start considering the continuity of the optimal solution with respect to the parameters. The parameter set of interest is

[TABLE]

Theorem 4.1**.**

Let

[TABLE]

*Then the map $({c,p,W})\mapsto\hat{q}:=\operatorname*{arg\,min}_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{{c,p,W}}(q)$ is continuous on $\mathcal{P}$. *

Proof 4.2**.**

Following the procedure in [36, Proposition 7.3] we use the continuity of the optimal value (Lemma A.1) to show continuity of the optimal solution. To this end, let $({c^{(k)},p^{(k)},W^{(k)}})$ be a sequence of parameters in $\mathcal{P}$ converging to $({c,p,W})\in\mathcal{P}$ as $k\to\infty$. Moreover, defining $\mathbb{J}_{k}(q):=\mathbb{J}_{c^{(k)},p^{(k)},W^{(k)}}(q)$ and $\mathbb{J}(q):=\mathbb{J}_{c,p,W}(q)$ for simplicity of notation, let $\hat{q}_{k}=\operatorname*{arg\,min}_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{k}(q)$ and $\hat{q}=\operatorname*{arg\,min}_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}(q)$. By Lemma A.1, $(\hat{q}_{k})$ is bounded, and hence there is a subsequence, which for simplicity we also call $(\hat{q}_{k})$, converging to a limit $q_{\infty}$. If we can show that $q_{\infty}=\hat{q}$, then the theorem follows. To this end, choosing a $q_{0}\in\mathfrak{P}_{+}$, we have

[TABLE]

Consequently, by Lemma A.1,

[TABLE]

*However $\hat{q}_{k}+\varepsilon q_{0}\in\mathfrak{P}_{+}$, and, since $({c,p,W},q)\mapsto\mathbb{J}_{{c,p,W}}(q)$ is continuous in $\mathcal{P}\times\mathfrak{P}_{+}$, we obtain *

[TABLE]

*Letting $\varepsilon\to 0$ in (4.3), we obtain the inequality $\mathbb{J}(\hat{q})\geq\mathbb{J}(q_{\infty})$. By strict convexity of $\mathbb{J}$ the optimal solution is unique, and hence $\hat{q}=q_{\infty}$. *

4.2 Continuity of $\hat{c}$ with respect to $\hat{q}$

We have now established continuity from $({c,p,W})$ to $\hat{q}$. In the same way as in Proposition 2.7 we are also interested in continuity of the map $({c,p,W})\mapsto(\hat{q},\hat{c})$. This would follow if we could show that the map from $\hat{q}$ to $\hat{c}$ is continuous. From the KKT condition (3.11c), it is seen that $\hat{r}$ is continuous in $c$, $W$ and $\hat{q}$. In view of (3.11b), i.e.,

[TABLE]

continuity of $\hat{c}$ would follow if $\int_{\mathbb{T}^{d}}P\hat{Q}^{-1}dm$ is continuous in $(p,\hat{q})$ whenever it is finite. If $p\in\mathfrak{P}_{+}$, this follows from the continuity the map $\hat{q}\mapsto\hat{Q}^{-1}$ in $L_{1}(\mathbb{T}^{d})$. For the case $d\leq 2$, this is trivial since if $\int_{\mathbb{T}^{d}}\hat{Q}^{-1}dm$ is finite, then $\hat{q}\in\mathfrak{P}_{+}$ and $\hat{Q}$ is bounded away from zero (cf., Proposition 2.7). However, for the case $d>2$ the optimal $\hat{q}$ may belong to the boundary $\partial\mathfrak{P}_{+}$, i.e., $\hat{Q}$ is zero in some point. The following proposition shows $L_{1}$ continuity of $\hat{q}\mapsto\hat{Q}^{-1}$ for certain cases.

Proposition 4.3**.**

*For $d\geq 3$, let $\hat{q}\in\bar{\mathfrak{P}}_{+}$ and suppose that the Hessian $\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\,\hat{Q}$ is positive definite in each point where $\hat{Q}$ is zero. Then $\hat{Q}^{-1}\in L_{1}(\mathbb{T}^{d})$ and the mapping from the coefficient vector $q\in\bar{\mathfrak{P}}_{+}$ to $Q^{-1}$ is $L_{1}$ continuous in the point $\hat{q}$. *

The proof of this proposition is given in [52]. From Propositions 4.3 and 2.7 the following continuity result follows directly.

Corollary 4.4**.**

*For all $c\in\mathfrak{C},p\in\mathfrak{P}_{+},W>0$, the mapping $(c,p,W)\to(\hat{q},\hat{c})$ is continuous in any point $({c,p,W})$ for which the Hessian $\nabla_{{\boldsymbol{\theta}}{\boldsymbol{\theta}}}\,\hat{Q}$ is positive definite in each point where $\hat{Q}$ is zero. *

The condition $p\in\mathfrak{P}_{+}$ is needed, since we may have pole-zero cancelations in $P/\hat{Q}$ when $p\in\partial\mathfrak{P}_{+}$, and then $\int_{\mathbb{T}^{d}}P/\hat{Q}dm$ may be finite even if $\hat{Q}^{-1}\not\in L_{1}$. The following example shows that this may lead to discontinuities in the map $p\mapsto\hat{c}$ (cf. Example 3.8 in [36]).

Example 4.5**.**

Let

[TABLE]

where $dm=d\theta/2\pi$ and $d\nu_{0}$ is the singular measure $\delta_{0}(\theta)d\theta$ with support in $\theta=0$. Since $d\mu:=2dm+d\nu_{0}$ is positive, $c\in\bar{\mathfrak{C}}_{+}$. Moreover, since

[TABLE]

we have that $c\in\mathfrak{C}_{+}$ (see, e.g., [41, p. 2853]). Thus we know [54, Corollary 2.3] that for each $p\in\mathfrak{P}_{+}$ we have a unique $\hat{q}\in\mathfrak{P}_{+}$ such that $P/\hat{Q}$ matches $c$, and hence $\hat{c}=0$. However, for $p=2(-1,2,-1)^{\prime}$ we have that $\hat{q}=(-1,2,-1)^{\prime}$ and $\hat{c}=(1,1,1)^{\prime}$ (Theorem 2.5). Then, for the sequence $(p_{k})$, where $p_{k}=2(-1,2+1/k,-1)\in\mathfrak{P}_{+}$, we have $\hat{c}_{k}=0$, so

[TABLE]

*which shows that the mapping $p\to\hat{c}$ is not continuous. *

5 Tuning to avoid a singular part

In many situations we prefer solutions where there is no singular measure $d\nu$ in the optimal solution. An interesting question is therefore for what prior $P$ and weight $W$ we obtain $d\hat{\nu}=0$. The following result provides a sufficient condition.

Proposition 5.1**.**

Let $c\in\mathfrak{C}$ and let $p$ be the Fourier coefficients of the prior $P$. If the weight satisfies333Here $\|A\|_{2,1}=\max_{c\neq 0}\|Ac\|_{1}/\|c\|_{2}$ denotes the subordinate (induced) matrix norm.

[TABLE]

then the optimal solution of (1.7) is on the form

[TABLE]

*i.e., the singular part $d\hat{\nu}$ vanishes. *

Remark 5.2**.**

Note that for a scalar weight, $W=\lambda I$ the bound (5.1) simplifies to

[TABLE]

*where $|\Lambda|$ is the cardinality of index set $\Lambda$. *

For the proof of Proposition 5.1 we need the following lemma.

Lemma 5.3**.**

Condition (5.1) implies

[TABLE]

*where $\hat{r}$ is the optimal value of $r$ in problem (1.7). *

Proof 5.4**.**

Let

[TABLE]

be the cost function of problem (1.7), and let $(d\hat{\mu},\hat{r})$ be the optimal solution. Clearly, $\mathbb{I}(Pdm,p)\geq\mathbb{I}(d\hat{\mu},\hat{r})$, and consequently

[TABLE]

since $\mathbb{D}(Pdm,d\hat{\mu})\geq 0$ and $\mathbb{D}(Pdm,Pdm)=0$. Therefore,

[TABLE]

*which is less than one by (5.1). Hence (5.1) implies (5.3). *

Proof 5.5** (Proof of Proposition 5.1).**

Suppose the optimal solution has a nonzero singular part $d\hat{\nu}$, and form the directional derivative of (5.4) at $(d\hat{\mu},\hat{r})$ in the direction $-d\hat{\nu}$. Then $\Phi$ in (1.3a) does not vary, and

[TABLE]

where

[TABLE]

Then $|\delta r_{\boldsymbol{k}}|\leq\int d\hat{\nu}$ for all ${\boldsymbol{k}}\in\Lambda$, and hence

[TABLE]

by (5.3) (Lemma 5.3). Consequently,

[TABLE]

*whenever $d\hat{\nu}\neq 0$, which contradicts optimality. Hence $d\hat{\nu}$ must be zero. *

The condition of Proposition 5.1 is just sufficient and is in general conservative. To illustrate this, we consider a simple one-dimensional example ($d=1$).

Example 5.6**.**

Consider a covariance sequence $(1,c_{1})$, where $c_{1}\neq 0$, and a prior $P(e^{i\theta})=1-\cos\theta$, and set $W=\lambda I$. Then, since

[TABLE]

the sufficient condition (5.2) for an absolutely continuous solution is

[TABLE]

We want to investigate how restrictive this condition is.

Clearly we will have a singular part if and only if $\hat{Q}=q_{0}P$, in which case we have

[TABLE]

for some $\beta>0$. In fact, it follows from $\langle\hat{c},\hat{q}\rangle=0$ in (3.11a) that $\hat{c}_{1}=\hat{c}_{0}$. Moreover, (3.11b) and (3.11c) yield

[TABLE]

By eliminating $\beta$, we get

[TABLE]

and solving for $q_{0}$ yields

[TABLE]

(note that $\lambda>0$ and $q_{0}>0$). Again, using (3.11c) we have

[TABLE]

We are interested in $\lambda$ for which $\beta>0$, i.e.,

[TABLE]

which is equivalent to the two conditions

[TABLE]

which could be seen by noting that the left member of (5.6) must be positive and then squaring both sides. To find out whether this has a solution we consider three cases, namely $c_{1}<-1/2$, $-1/2<c_{1}<0$, and $c_{1}>0$. For $c_{1}<-1/2$, condition (5.7) becomes $2c_{1}<\lambda<1+5c_{1}$, which is impossible since $1+5c_{1}<2c_{1}$. Condition (5.7) cannot be satisfied when $-1/2<c_{1}<0$, because then $\lambda$ would be negative which contradicts $\lambda>0$. When $c_{1}>0$, Condition (5.7) is satisfied if and only if $\lambda<2c_{1}$.

Consequently, there is no singular part if either $c_{1}$ is negative or

[TABLE]

*This shows that the condition (5.5) is not tight. *

6 Covariance extension with hard constraints

The alternative optimization problem (1.9) amounts to minimizing $\mathbb{D}(Pdm,d\mu)$ subject to the hard constraint $\|r-c\|_{W^{-1}}^{2}\leq 1$, where $r_{\boldsymbol{k}}=\int_{\mathbb{T}^{d}}e^{i({\boldsymbol{k}},{\boldsymbol{\theta}})}d\mu$. Hard constraints of this type were used in [57] in the context of entropy maximization. In general the data $c\not\in\bar{\mathfrak{C}}_{+}$, whereas, by definition, $r\in\bar{\mathfrak{C}}_{+}$. Consequently, a necessary condition for the existence of a solution is that $\bar{\mathfrak{C}}_{+}$ and the strictly convex set

[TABLE]

have a nonempty intersection. In the case that $\mathfrak{S}_{W}\cap\bar{\mathfrak{C}}_{+}\subset\partial\mathfrak{C}_{+}$, this intersection only contains one point [44, Section 3.12]. In this case, any solution to the moment problem contains only a singular part (Proposition 2.4), and then the primal problem (1.9) has a unique feasible point $r$, but the objective function is infinite. Moreover, $\mathbb{D}(Pdm,d\mu)\geq 0$ is strictly convex with $\mathbb{D}(Pdm,Pdm)=0$, so if $p\in\mathfrak{S}_{W}$ then (1.9) has the trivial unique optimal solution $d\hat{\mu}=Pdm$, and $\hat{r}=p$. The remaining case, $p\not\in\mathfrak{S}_{W}\cap\mathfrak{C}_{+}\neq\emptyset$ needs further analysis.

To this end, setting $d\mu=\Phi dm+d\nu$, we consider the Lagrangian

[TABLE]

where $\gamma\geq 0$. Therefore, in view of (3.1),

[TABLE]

where, as before, $e:=[e_{\boldsymbol{k}}]_{{\boldsymbol{k}}\in\Lambda}$, $e_{\boldsymbol{0}}=1$ and $e_{{\boldsymbol{k}}}=0$ for ${\boldsymbol{k}}\in\Lambda\setminus\{\boldsymbol{0}\}$, and hence $r_{\boldsymbol{0}}=\langle r,e\rangle$. This Lagrangian differs from that in (3) only in the last term that does not depend on $\Phi$. Therefore, in deriving the dual functional

[TABLE]

we only need to consider $q\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$, and a first variation in $\Phi$ yields (1.3b) and (3.3). The directional derivative

[TABLE]

is zero for

[TABLE]

Thus inserting (1.3b) and (3.3) and (6.3) into (6) yields the dual functional

[TABLE]

to be minimized over all $q\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$ and $\gamma\geq 0$. Since $\frac{d\varphi}{d\gamma}=-\frac{1}{4\gamma^{2}}\|q-e\|_{W}^{2}+1$, there is a stationary point

[TABLE]

that is nonnegative as required.

For $\gamma=0$ we must have $q=e$, and consequently $\varphi(q,\gamma)$ tends to zero as $\gamma\to 0$. By weak duality zero is therefore a lower bound for the minimization problem (1.9), and $\mathbb{D}(Pdm,d\hat{\mu})=0$, which corresponds to the trivial unique solution $d\hat{\mu}=Pdm$ and $\hat{r}=p$ mentioned above. This solution is only feasible if $p\in\mathfrak{S}_{W}$. Therefore we can restrict our attention to the case $\gamma>0$. Inserting (6.5) into (6.4) and removing the constant term $c_{\boldsymbol{0}}$, we obtain the modified dual functional

[TABLE]

Moreover, combining (6.3) and (6.5), we obtain

[TABLE]

which also follows from complementary slackness since $\gamma>0$ and restricts $r$ to the boundary of $\mathfrak{S}_{W}$.

Theorem 6.1**.**

Suppose that $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$, $p\not\in\mathfrak{S}_{W}$ and $\mathfrak{S}_{W}\cap\mathfrak{C}_{+}\neq\emptyset$. Then the modified dual problem

[TABLE]

has a unique solution $\hat{q}\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$. Moreover, there exists a unique $\hat{r}\in\mathfrak{C}_{+}$, a unique $\hat{c}\in\partial\mathfrak{C}_{+}$ and a (not necessarily unique) nonnegative singular measure $d\hat{\nu}$ with support

[TABLE]

such that

[TABLE]

and the measure

[TABLE]

is an optimal solution to the primal problem (1.9). Moreover,

[TABLE]

and $d\hat{\nu}$ can be chosen with support in at most $|\Lambda|-1$ points.

*If $p\in\mathfrak{S}_{W}$, the unique optimal solution is $d\hat{\mu}=Pdm$, and then $\hat{r}=p$. If $\mathfrak{S}_{W}\cap\bar{\mathfrak{C}}_{+}\subset\partial\mathfrak{C}_{+}$, any solution to the moment problem will have only a singular part. Finally, if $\mathfrak{S}_{W}\cap\bar{\mathfrak{C}}_{+}=\emptyset$, then the problem (1.9) will have no solution. *

Proof 6.2**.**

We begin by showing that the functional $\mathbb{J}$ has a minimum under the stated conditions. To this end, we first establish that the functional $\mathbb{J}$ has compact sublevel sets $\mathbb{J}^{-1}(-\infty,\rho]$, i.e., $\|q\|_{\infty}$ is bounded for all $q$ such that $\mathbb{J}(q)\leq\rho$, where $\rho$ is sufficiently large for the sublevel set to be nonempty. The functional (6.6) can be decomposed in a linear and a logarithmic term as

[TABLE]

where $h(q):=\langle c,q-e\rangle+\|q-e\|_{W}$. The integral term will tend to $-\infty$ as $\|q\|_{\infty}\to\infty$. Therefore we need to have the linear term to tend to $+\infty$ as $\|q\|_{\infty}\to\infty$, in which case we can appeal to the fact that linear growth is faster than logarithmic growth. However, if $c\not\in\bar{\mathfrak{C}}_{+}$ as is generally assumed, there is a $q\in\bar{\mathfrak{P}}_{+}$ such that $\langle c,q\rangle<0$, so we need to ensure that the positive term $\|q-e\|_{W}$ dominates.

Let $\tilde{r}\in\mathfrak{S}_{W}\cap\mathfrak{C}_{+}\neq\emptyset$. Then, by Theorem 2.5, there is a positive measure $d\tilde{\mu}=\tilde{\Phi}dm+d\tilde{\nu}$ with a nonzero $\tilde{\Phi}$ such that

[TABLE]

and $\tilde{r}$ satisfies the constraints in the primal problem (1.9). Consequently,

[TABLE]

for all $q\in\bar{\mathfrak{P}}_{+}$ and $\gamma\geq 0$, which in particular implies that

[TABLE]

Now, if there is a $q\in\bar{\mathfrak{P}}_{+}$ such that $h(q)\leq 0$, then $\mathbb{J}(\lambda q)\to-\infty$ as $\lambda\to\infty$, which contradicts (6.13). Therefore, $h(q)>0$ for all $q\in\bar{\mathfrak{P}}_{+}$. Then, since $h$ is continuous, it has a minimum $\varepsilon$ on the compact set $K:=\{q\in\bar{\mathfrak{P}}_{+}\setminus\{0\}\mid\|q-e\|_{\infty}=1\}$. As $e\not\in K$, $\epsilon>0$. Therefore,

[TABLE]

since $\|Q\|_{\infty}\leq|\Lambda|\|q\|_{\infty}$ [54, Lem. A.1]. Likewise,

[TABLE]

since $Q/\|Q\|_{\infty}\leq 1$. Hence

[TABLE]

Comparing linear and logarithmic growth we see that the sublevel set is bounded from above and below. Moreover, a trivial modification of [54, Lemma 3.1] shows that $\mathbb{J}$ is lower semi-continuous, and hence $\mathbb{J}^{-1}(-\infty,\rho]$ is compact. Consequently, the problem (6.8) has an optimal solution $\hat{q}$.

Next we show that $\hat{q}$ is unique. For this we return to the original dual problem to find a minimum of (6.4). The solution $\hat{q}$ is a minimizer of $\varphi(q,\hat{\gamma})$, where

[TABLE]

and $\mathbb{J}(\hat{q})=\varphi(\hat{q},\hat{\gamma})+c_{\boldsymbol{0}}$. To show that $\varphi$ is strictly convex, we form the Hessian

[TABLE]

and the quadratic form

[TABLE]

which is positive for all nonzero $(x,\xi)$, since $(q-e)\neq 0$ and $\gamma>0$. Consequently, $\varphi$ has a unique minimizer $(\hat{q},\hat{\gamma})$, where $\hat{q}$ is the unique minimizer of $\mathbb{J}$.

It follows from (6.3) and (6.5) that

[TABLE]

which consequently is unique. Moreover, $h(\hat{q})=\langle\hat{r},\hat{q}\rangle-\hat{r}_{\boldsymbol{0}}$, and hence we can follow the same line of proof as in Theorem 3.1 to show that there is a unique $\hat{c}\in\partial\mathfrak{C}_{+}$ such that $\langle\hat{c},\hat{q}\rangle=0$ and a positive discrete measure $d\hat{\nu}$ with support in $|\Lambda|-1$ points so that (6.9) and (6.10) hold. Next, let $\mathbb{I}(d\mu)=-\mathbb{D}(Pdm,d\mu)$ be the primal functional in (1.9), where $d\mu$ is restricted to the set of positive measures $d\mu:=\Phi dm+d\nu$ such that $r$, given by (1.5), satisfies the constraint $\|r-c\|_{W}\leq 1$. In view of (6.12),

[TABLE]

*for any such $d\mu$, and hence $d\hat{\mu}$ is an optimal solution to the primal problem (1.9). Finally, the cases $p\in\mathfrak{S}_{W}$, $\mathfrak{S}_{W}\cap\bar{\mathfrak{C}}_{+}\subset\partial\mathfrak{C}_{+}$, and $\mathfrak{S}_{W}\cap\bar{\mathfrak{C}}_{+}=\emptyset$ have already been discussed above. *

Corollary 6.3**.**

Suppose that $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$ and $\mathfrak{S}_{W}\cap\mathfrak{C}_{+}\neq\emptyset$. The KKT conditions

[TABLE]

*are necessary and sufficient conditions for optimality of the dual pair (1.9) and (6.8) of optimization problems. *

The corollary follows by noting that, if $p\in\mathfrak{S}_{W}$, then we obtain the trivial solution $\hat{q}=e$, which corresponds to the primal optimal solution $d\hat{\mu}=Pdm$.

Proposition 6.4**.**

The condition

[TABLE]

*is sufficient for the pair (1.9) and (6.8) of dual problems to have optimal solutions. *

Proof 6.5**.**

*If $W>cc^{*}$, then $(q-e)^{*}W(q-e)\geq\langle c,q-e\rangle^{2}$ with equality only for $q=e$. Hence, if $q\neq e$, $\|q-e\|_{W}>|\langle c,q-e\rangle|$, i.e., $h(q)>0$ for all $q\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$ except $q=e$. Then we proceed as in the proof of Theorem 6.1. *

Remark 6.6**.**

Condition (6.17) guarantees that $0\in{\rm int}(\mathfrak{S}_{W})$ and hence in particular that $\mathfrak{S}_{W}\cap\mathfrak{C}_{+}\neq\emptyset$ as required in Theorem 6.1. To see this, note that $0\in\bar{\mathfrak{C}}_{+}$ and that $r=0$ satisfies the hard constraint in (1.9) if $c^{*}W^{-1}c\leq 0$. However, since $W>cc^{*}$, there is a $W_{0}>0$ such that $W=W_{0}+cc^{*}$. Then the well-known Matrix Inversion Lemma (see, e.g., [42, p. 746]) yields

[TABLE]

and therefore

[TABLE]

*which establishes that $0\in{\rm int}(\mathfrak{S}_{W})$. However, for $\mathfrak{S}_{W}\cap\mathfrak{C}_{+}$ to be nonempty, $r=0$ need not be contained in this set. Hence, condition (6.17) is not necessary, although it is easily testable. In fact, this provides an alternative proof of Proposition 6.4. *

7 On the equivalence between the two problems

Clearly $\mathfrak{S}_{W}\cap\mathfrak{C}_{+}$ is always nonempty if $c\in\mathfrak{C}_{+}$. Then both the problem (1.7) with soft constraints and the problem (1.9) with hard constraints have a solution for any choice of W. On the other hand, if $c\not\in\mathfrak{C}_{+}$, the problem with soft constraints will always have a solution, while the problem with hard constraints may fail to have one for certain choices of $W$. However, if the weight matrix in the hard-constrained problem – let us denote it $W_{\rm hard}$ – is chosen in the set $\mathcal{W}:=\{W>0\mid\mathfrak{S}_{W}\cap\mathfrak{C}_{+}\neq\emptyset,p\not\in\mathfrak{S}_{W}\}$, then it can be seen from Corollaries 3.3 and 6.3 that we obtain exactly the same solution $\hat{q}$ in the soft-constrained problem by choosing

[TABLE]

We note that (7.1) can be written $W_{\rm hard}=\alpha W_{\rm soft}$, where $\alpha:=\|\hat{q}-e\|_{W_{\rm hard}}$. Therefore, substituting $W_{\rm hard}$ in (7.1), we obtain

[TABLE]

which yields $\alpha=\|\hat{q}-e\|_{W_{\rm soft}}^{2}$. Hence the inverse of (7.1) is given by

[TABLE]

By Theorem 4.1 $\hat{q}$ is continuous in $W_{\rm soft}$, and hence, by (7.2), the corresponding $W_{\rm hard}$ varies continuously with $W_{\rm soft}$. In fact, this can be strengthened to a homeomorphism between the two weight matrices.

Theorem 7.1**.**

*The map (7.1) is a homeomorphism between $\mathcal{W}$ and the space of all (Hermitian positive definite) weight matrices, and the inverse is given by (7.2). *

Proof 7.2**.**

By [11, Lemma 2.3], a continuous map between two spaces of the same dimension is a homeomorphism if and only if it is injective and proper, i.e., the preimage of any compact set is compact. To see that $\mathcal{W}$ is open, we observe that $\mathfrak{S}_{W}$ is continuous in $W$ and that $\mathfrak{C}_{+}$ is an open set. As noted above, the map (7.2) – let us call it $f$ – is continuous and also injective, as it can be inverted. Hence it only remains to show that $f$ is proper. To this end, we take a compact set $K\subset\mathcal{W}$ and show that $f^{-1}(K)$ is also compact. There are two ways this could fail. First, the preimage could contain a singular semidefinite matrix. However this is impossible by (7.2), since $\|\hat{q}\|_{\infty}$ is bounded for $W_{\text{hard}}\in K$ (Lemma A.3) and a nonzero scaling of a singular matrix cannot be nonsingular. Secondly, $\|W_{\rm soft}\|_{F}$ could tend to infinity. However, this is also impossible. To see this, we first show that there is a $\kappa>0$ such that $\|p-r\|_{W_{\rm hard}^{-1}}\geq\kappa$ for all $r\in\mathfrak{S}_{W_{\rm hard}}$ and all $W_{\rm hard}\in K$. To this end, we observe that the minimum of $\|p-r\|_{W^{-1}}$ over all $W\in K$ and $r$ satisfying the constraint $\|r-c\|_{W^{-1}}\leq 1$ is bounded by

[TABLE]

*by the triangle inequality $\|p-r\|_{W^{-1}}\geq\|p-c\|_{W^{-1}}-\|c-r\|_{W^{-1}}\geq\|p-c\|_{W^{-1}}-1$. The minimum is attained, since $K$ is compact, and positive, since $p\not\in\bigcup_{W\in K}\mathfrak{S}_{W}$. Now, from Corollary 6.3 we see that $\hat{q}=e$ if and only if $\hat{r}=p$. The map from $\hat{q}\mapsto\hat{r}$ is continuous in $q=e$. In fact, $\hat{Q}$ is uniformly positive in a neighborhood of $e$ and hence the corresponding $\hat{c}=0$. Due to this continuity, if $\hat{q}\to e$, then $\hat{r}\to p$, which cannot happen since $\|p-r\|_{W^{-1}}\geq\kappa$ for all $W\in K$. Thus, since $\|\hat{q}-e\|_{W}$ is bounded away from zero, the preimage $f^{-1}(K)$ of $K$ is bounded. Finally, consider a convergent sequence $(W_{k})$ in $f^{-1}(K)$ converging to a limit $W_{\infty}$. Since the sequence is bounded and cannot converge to a singular matrix, we must have $W_{\infty}>0$, i.e., $W_{\infty}\in f^{-1}(\mathcal{W})$. By continuity, $f(W_{k})$ tends to the limit $f(W_{\infty})$, which must belong to $K$ since it is compact. Hence the preimage $W_{\infty}$ must belong to $f^{-1}(K)$. Therefore, $f^{-1}(K)$ is compact as claimed. *

It is illustrative to consider the simple case when $W=\lambda I$. Then the two maps (7.1) and (7.2) become

[TABLE]

Whereas the range of $\lambda_{\rm soft}$ is the semi-infinite interval $(0,\infty)$, for the homeomorphism to hold $\lambda_{\rm hard}$ is confined to

[TABLE]

where $\lambda_{\rm min}$ is the distance from $c$ to the cone $\bar{\mathfrak{C}}_{+}$ and $\lambda_{\rm max}=\|c-p\|$. When $\lambda_{\rm soft}\to\infty$, $\lambda_{\rm hard}\to\lambda_{\rm max}$ and $\hat{q}\to e$. If $\lambda_{\rm hard}\geq\lambda_{\rm max}$, then the coresponding problem has the trivial unique solution $\hat{q}=e$, corresponding to the primal solution $d\hat{\mu}=Pdm$.

Note that Theorem 7.1 implies that some continuity results in one of the problems can be automatically transferred to the other problem. In particular, we have the following result.

Theorem 7.3**.**

Let

[TABLE]

*Then the map $W\mapsto\hat{q}:=\operatorname*{arg\,min}_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{W}(q)$ is continuous. *

Proof 7.4**.**

*The theorem follows by noting that $W\mapsto\hat{q}:=\operatorname*{arg\,min}_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{W}(q)$ can be seen as a composition of two continuous maps, namely the one in Theorem 4.1 and the one in Theorem 7.1. *

Next we shall vary also $c$ and $p$, and to this end we introduce a more explicit notation for $\mathfrak{S}_{W}$ and $\mathcal{W}$, namely $\mathfrak{S}_{c,W}=\mathfrak{S}_{W}$ in (6.1) and

[TABLE]

Then the corresponding set of parameters (4.1) for the problem with hard constraints is given by

[TABLE]

the interior of which is

[TABLE]

Theorem 7.1 can now be modified accordingly to yield the following theorem, the proof of which is deferred to the appendix.

Theorem 7.5**.**

Let the map $(c,p,W_{\rm hard})\mapsto W_{\rm soft}$ be given by (7.1) and the map $(c,p,W_{\rm soft})\mapsto W_{\rm hard}$ by (7.2). Then the map that sends $(c,p,W_{\rm hard})\in{\rm int}(\mathcal{P}_{\rm hard})$ to $(c,p,W_{\rm soft})\in{\rm int}(\mathcal{P})$ is a homeomorphism.

Note that this theorem is not a strict amplification of Theorem 7.1 as we have given up the possibility for $p$ to be on the boundary $\partial\mathfrak{P}_{+}$. The same is true for the following modification of Theorem 7.3.

Theorem 7.6**.**

*Let $\mathbb{J}_{{c,p,W}}(q)$ be as in (7.4). Then the map $({c,p,W})\mapsto\hat{q}:=\operatorname*{arg\,min}_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{{c,p,W}}(q)$ is continuous on ${\rm int}(\mathcal{P}_{\rm hard})$. *

Proof 7.7**.**

*The theorem follows immediately by noting that $({c,p,W}_{\text{hard}})\mapsto\hat{q}$ can be seen as a composition of two continuous maps, namely $({c,p,W}_{\text{hard}})\mapsto({c,p,W}_{\text{soft}})$ of Theorem 7.5 and $({c,p,W}_{\text{soft}})\mapsto\hat{q}$ of Theorem 4.1. *

Theorem 7.6 is a counterpart of Theorem 4.1 for the problem with hard constraints, except that $p$ is restricted to the interior $\mathfrak{P}_{+}$. It should be possible to extend the result to hold for all $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$ via a direct proof along the lines of the proof of Theorem 4.1.

8 Estimating covariances from data

For a scalar stationary stochastic process $\{y(t);\,t\in\mathbb{Z}\}$, it is well-known that the biased covariance estimate

[TABLE]

based on an observation record $\{y_{t}\}_{t=0}^{N-1}$, yields a positive definite Toeplitz matrix, which is equivalent to $c\in\mathfrak{C}_{+}$ [2, pp. 13-14] In fact, these estimates correspond to the ones obtained from the periodogram estimate of the spectrum (see, e.g., [59, Sec. 2.2]). On the other hand, the Toeplitz matrix of the unbiased estimate

[TABLE]

is in general not positive definite.

The same holds in higher dimensions ($d>1$) where the observation record is $\{y_{\boldsymbol{t}}\}_{{\boldsymbol{t}}\in\mathbb{Z}^{d}_{N}}$ with

[TABLE]

The unbiased estimate is then given by

[TABLE]

and the biased estimate by

[TABLE]

where we define $y_{\boldsymbol{t}}=0$ for ${\boldsymbol{t}}\notin\mathbb{Z}^{d}_{N}$. The sequence of unbiased covariance estimates does not in general belong to $\mathfrak{C}_{+}$, but the biased covariance estimates yields $c\in\mathfrak{C}_{+}$ also in the multidimensional setting. In fact, this can be seen by noting that the biased estimate corresponds to the Fourier coefficients of the periodogram [18, Sec. 6.5.1], i.e., if the estimates $c_{\boldsymbol{k}}$ are given by (8.2), then

[TABLE]

where $\mathbb{Z}^{d}_{N}-\mathbb{Z}^{d}_{N}$ denotes the Minkowski set difference. This leads to the following lemma.

Lemma 8.1**.**

*Given the observed data $\{y_{{\boldsymbol{t}}}\}_{{\boldsymbol{t}}\in\mathbb{Z}^{d}_{N}}$, let $\{c_{\boldsymbol{k}}\}_{{\boldsymbol{k}}\in\Lambda}$ be given by (8.2). Then $c\in\mathfrak{C}_{+}$. *

Proof 8.2**.**

Given $\{y_{{\boldsymbol{t}}}\}_{{\boldsymbol{t}}\in\mathbb{Z}^{d}_{N}}$, let $c=\{c_{\boldsymbol{k}}\}_{{\boldsymbol{k}}\in\mathbb{Z}^{d}_{N}}$, where $c_{\boldsymbol{k}}$ be given by (8.2). In view of (2.1) and (8.3) we have

[TABLE]

*which is positive for all $p\in\bar{\mathfrak{P}}_{+}\setminus\{0\}$. Consequently $c\in\mathfrak{C}_{+}$. *

An advantage of the approximate procedures to the rational covariance extension problem is that they can also be used for cases where the biased estimate is not available, e.g., where the covariance is estimated from snapshots.

9 Application to spectral estimation

As long as we use the biased estimate (8.2), we may apply exact covariance matching as outlined in Section 2, whereas in general approximate covariance matching will be required for biased covariance estimates. However, as will be seen in the following example, approximate covariance matching may sometimes be better even if $c\in\mathfrak{C}_{+}$.

In this application it is easy to determine a bound on the acceptable error in the covariance matching, so we use the procedure with hard constraints. Given data generated from a two-dimensional stochastic system, we test three different procedures, namely (i) using the biased estimate and exact matching, (ii) using the biased estimate and the approximate matching (1.9), and (iii) using the unbiased estimate and the approximate matching (1.9). The procedures are then evaluated by checking the size of the error between the matched covariances and the true ones from the dynamical system.

9.1 An example

Let $y_{(t_{1},t_{2})}$ be the steady-state output of a two-dimensional recursive filter driven by a white noise input $u_{(t_{1},t_{2})}$. Let the transfer function of the recursive filter be

[TABLE]

where $\Lambda_{+}=\{(k_{1},k_{2})\in\mathbb{Z}^{2}\mid 0\leq k_{1}\leq 2,0\leq k_{2}\leq 2\}$ and the coefficients are given by $b_{(k_{1},k_{2})}=B_{k_{1}+1,k_{2}+1}$ and $a_{(k_{1},k_{2})}=A_{k_{1}+1,k_{2}+1}$, where

[TABLE]

The spectral density $\Phi$ of $y_{(t_{1},t_{2})}$, which is shown in Fig. 1 and is similar to the one considered in [54], is given by

[TABLE]

and hence the index set $\Lambda$ of the coefficients of the trigonometric polynomials $P$ and $Q$ is given by $\Lambda=\Lambda_{+}-\Lambda_{+}=\{(k_{1},k_{2})\in\mathbb{Z}^{2}\,|\;|k_{1}|\leq 2,|k_{2}|\leq 2\}$. Using this example, we perform two different simulation studies.

9.2 First simulation study

The system was simulated for $500$ time steps along each dimension, starting from $y_{(t_{1},t_{2})}=u_{(t_{1},t_{2})}=0$ whenever either $t_{1}<0$ or $t_{2}<0$. Then covariances were estimated from the $9\times 9$ last samples, using both the biased and the unbiased estimator. With this covariance data we investigate the three procedures (i), (ii) and (iii) described above. In each case, both the maximum entropy (ME) solutions and solutions with the true numerator are computed.444Maximum entropy: $P\equiv 1$. True numerator: $P=P_{\text{true}}$. The weighting matrix is taken to be $W=\lambda I$, where $\lambda$ is $\lambda_{\text{biased}}:=\|c_{\text{true}}-c_{\text{biased}}\|_{2}^{2}$ in procedure (ii) and $\lambda_{\text{unbiased}}:=\|c_{\text{true}}-c_{\text{unbiased}}\|_{2}^{2}$ in procedure (iii).555Note that this is the smallest $\lambda$ for which the true covariance sequence belongs to the uncertainty set $\{r\,|\,\|r-c\|_{2}^{2}\leq\lambda\}$. The norm of the error666Here we use the norm of the covariance estimation error as measure of fit. However, note that this is not the only way to compare accuracy of the different methods. The reason for this choice is that comparing the accuracy of the spectral estimates is not straightforward since it depends on the selected metric or distortion measure. between the matched covariances and the true ones, $\|\hat{r}-c_{\text{true}}\|_{2}$, is shown in Table. 1. The means and standard deviations are computed over the $100$ runs.

The biased covariance estimates belong to the cone $\mathfrak{C}_{+}$ (Lemma 8.1), and therefore procedure (i) can be used. The corresponding error in Table 1 is the statistical error in estimating the covariance. This error is quite large because of a short data record. Using approximate covariance matching in this case seems to give a worse match. However, approximate matching of the unbiased covariances gives as good a fit as exact matching of the biased ones.

9.3 Second simulation study

In this simulation the setup is the same as the previous one, except that the simulation data has been discarded if the unbiased estimate belongs to $\bar{\mathfrak{C}}_{+}$. To obtain $100$ such data sets, $414$ simulations of the system were needed. (As a comparison, in the previous experiment $23$ out of the $100$ runs resulted in an unbiased estimate outside $\bar{\mathfrak{C}}_{+}$.) Again, the norm of the error between matched covariances and the true ones is shown in Table 2, and the means and standard deviations are computed over the $100$ runs.

As before, the biased covariance estimates belong to the cone $\mathfrak{C}_{+}$, and therefore procedure (i) can be used. Comparing this with the results from procedure (ii) suggests that there may be an advantage not to enforce exact matching, although we know that the data belongs to the cone. Regarding procedure (iii), we know that the unbiased covariance estimates do not belong to the cone $\bar{\mathfrak{C}}_{+}$, hence we need to use approximate covariance matching. In this example, this procedure turns out to give the smallest estimation error.

10 Application to system identification and texture reconstruction

Next we apply the theory of this paper to texture generation via Wiener system identification. Wiener systems form a class of nonlinear dynamical systems consisting of a linear dynamic part composed with a static nonlinearity as illustrated in Figure 2. This is a subclass of so called block-oriented systems [4], and Wiener system identification is a well-researched area (see, e.g., [32] and references therein) that is still very active [43, 60, 1]. Here, we use Wiener systems to model and generate textures.

Using dynamical systems for modeling of images and textures is not new and has been considered in, e.g., [14, 50]. The setup presented here is motivated by [23], where thresholded Gaussian random fields are used to model porous materials for design of surface structures in pharmaceutical film coatings. Hence we let the static nonlinearity, call it $f$, be a thresholding with unknown thresholding parameter $\tau$. In our previous work [52] we applied exact covariance matching to such a problem. However, in general there is no guarantee that the estimated covariance sequence $c$ belongs to the cone $\mathfrak{C}_{+}$. Consequently, here we shall use approximate covariance matching instead.

The Wiener system identification can be separated into two parts. We start by identifying the nonlinear part. Using the notations of Figure 2, let $\{u_{\boldsymbol{t}};\,{\boldsymbol{t}}\in\mathbb{Z}^{d}\}$ be a zero-mean Gaussian white noise input, and let $\{x_{\boldsymbol{t}};\,{\boldsymbol{t}}\in\mathbb{Z}^{d}\}$ be the stationary output of the linear system, which we assume to be normalized so that $c_{\mathbf{0}}:=\operatorname*{\mathbb{E}}[x_{\boldsymbol{t}}^{2}]=1$. Moreover, let $y_{\boldsymbol{t}}=f(x_{\boldsymbol{t}})$ where $f$ is the static nonlinearity

[TABLE]

with unknown thresholding parameter $\tau$. Since $\operatorname*{\mathbb{E}}[y_{\boldsymbol{t}}]=1-\phi(\tau)$, where $\phi(\tau)$ is the Gaussian cumulative distribution function, an estimate of $\tau$ is given by $\tau_{\rm est}=\phi^{-1}(1-\operatorname*{\mathbb{E}}[y_{\boldsymbol{t}}])$.

Now, let $c_{\boldsymbol{k}}^{x}:=\operatorname*{\mathbb{E}}[x_{{\boldsymbol{t}}+{\boldsymbol{k}}}x_{\boldsymbol{t}}]$ be the covariances of $x_{\boldsymbol{t}}$, and let $c_{\boldsymbol{k}}^{y}:=\operatorname*{\mathbb{E}}[y_{{\boldsymbol{t}}+{\boldsymbol{k}}}y_{\boldsymbol{t}}]-\operatorname*{\mathbb{E}}[y_{{\boldsymbol{t}}+{\boldsymbol{k}}}]\operatorname*{\mathbb{E}}[y_{\boldsymbol{t}}]$ be the covariances of $y_{\boldsymbol{t}}$. As was explained in [52], by using results from [51] one can obtain a relation between $c_{\boldsymbol{k}}^{y}$ and $c_{\boldsymbol{k}}^{x}$, given by

[TABLE]

This is an invertible map, which we compute numerically, and given $\tau_{\rm est}$ we can thus get estimates of the covariances $c_{\boldsymbol{k}}^{x}$ from estimates of the covariances $c_{\boldsymbol{k}}^{y}$. However, even if $c^{y}$ is is a biased estimate so that $c^{y}\in\mathfrak{C}_{+}$, $c^{x}$ may not be a bona fide covariance sequence.

10.1 Identifying the linear system

Solving (1.7) or (1.9) for a given sequence of covariance estimates $c$, we obtain an estimate of the absolutely continuous part of the power spectrum $\Phi$ of that process. In the case $d=1$, $\Phi=P/Q$ can be factorized as

[TABLE]

which provides a transfer function of a corresponding linear system, which fed by a white noise input will produce an autoregressive-moving-average (ARMA) process with an output signal with precisely the power distribution $\Phi$ in steady state. For $d\geq 2$, a spectral factorization of this kind is not possible in general [19], but instead there is always a factorization as a sum-of-several-squares [17, 28],

[TABLE]

the interpretation of which in terms of a dynamical system is unclear when $m>1$. Therefore we resort to a heuristic and apply the factorization procedure in [29, Theorem 1.1.1] although some of the conditions required to ensure the existence of a spectral factor may not be met. (See [54, Section 7] for a more detailed discussion.)

10.2 Simulation results

The method, which is summarized in Algorithm 1, is tested on some textures from the Outex database [49] (available online at http://www.outex.oulu.fi/). These textures are color images and have thus been converted to binary textures by first converting them to black-and-white and then thresholding them.777The algorithm has been implemented and tested in Matlab, version R2015b. The textures have been normalized to account for light inhomogenities using a reference image available in the database. The conversion from color images to black-and-white images was done with the built-in function rgb2gray, and the threshold level was set to the mean value of the maximum and minimum pixel value in the black-and-white image. Three such textures are shown in Figure 3a through 3c.

In this example there is no natural bound on the error, so we use the problem with soft constraints, for which we choose the weight $W=\lambda I$ with $\lambda=0.01$ for all data sets. Moreover, we do maximum-entropy reconstructions, i.e., we set the prior to $P\equiv 1$. The optimization problems are then solved by first discretizing the grid $\mathbb{T}^{2}$, in this case in $50\times 50$ points (cf. [54, Theorem 2.6]), and solving the corresponding problems using the CVX toolbox [31, 30]. The reconstructions are shown in Figures 3d - 3f. Each reconstruction seems to provide a reasonable visual representation of the structure of the corresponding original. This is especially the case for the second texture.

11 Conclusions

In this work we extend the results of our previous paper [54] on the multidimensional rational covariance extension problem to allow for approximate covariance matching. We have provided two formulations to this problem, and we have shown that they are connected via a homeomorphism. In both formulations we have used weighted 2-norms to quantify the missmatch of the estimated covariances. However, we expect that by suitable modifications of the proofs similar results can be derived for other norms, since all norms have directional derivatives in each point [16, p. 49].

These results provide a procedure for multidimensional spectral estimation, but in order to obtain a complete theory for multidimensional system identification and realization theory there are still some open problems, such as spectral factorization and interpretations in terms of multidimensional stochastic systems, as briefly discussed in Section 10.1.

Appendix A Deferred proofs

Let $B_{\rho}(x^{(0)})$ denote the closed ball $\{x\in X\mid\|x-x^{(0)}\|_{X}\leq\rho\}$, where $X$ is either a set of vectors or a set of matrices depending on then context. The norm $\|\cdot\|_{X}$ is the Euclidean norm for vectors and Frobenius norm for matrices.

Lemma A.1**.**

*Let $\mathcal{P}$ be given by (4.1) and $\mathbb{J}_{{c,p,W}}$ by (4.2). Furthermore, let $\hat{q}:=\min_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{{c,p,W}}(q)$. Then the map $({c,p,W})\mapsto\mathbb{J}_{{c,p,W}}(\hat{q})$ is continuous for $({c,p,W})\in\mathcal{P}$. Moreover, for any compact $K\subset\mathcal{P}$, the corresponding set of optimal solutions $\hat{q}$ is bounded. *

Proof A.2**.**

The proof follows along the lines of Lemma 7.2 and Proposition 7.4 in [36]. Let $({c^{(0)},p^{(0)},W^{(0)}})\in\mathcal{P}$ be arbitrary and let

[TABLE]

where $\rho>0$ is chosen so that $\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})\subset\mathcal{P}$, i.e., $\rho<\|p^{(0)}\|_{2}$ and $W>0$ for all $\|W-W^{(0)}\|_{F}\leq\rho$. First we will show that the minimizer $\hat{q}_{{c,p,W}}$ of $\mathbb{J}_{{c,p,W}}$ is bounded for all $({c,p,W})\in\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$. To this end, note that by optimality

[TABLE]

and hence $\mathbb{J}_{c,p,W}(p)$ is bounded from above on the compact set $\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$. Consequently, by using the same inequality as in the proof of [36, Lemma 7.1], we see that

[TABLE]

Due to norm equivalence between $\|Q\|_{\infty}$ and $\|q\|_{W}$, and since the quadratic term is dominating, the norm of $\hat{q}_{{c,p,W}}$ is bounded in the set $\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$.

Now, let $K\subset\mathcal{P}$ be compact. We want to show that $\hat{q}$ is bounded on $K$. Assume it is not. Then let $({c^{(k)},p^{(k)},W^{(k)}})\in K$ be a sequence with $\|\hat{q}_{k}\|\to\infty$. Since $K$ is compact there is a converging subsequence $({c^{(k)},p^{(k)},W^{(k)}})\to({c,p,W})\in K$ with $\|\hat{q}_{k}\|\to\infty$. Since $({c,p,W})\in K$ there is a $\rho>0$ such that $\tilde{B}_{\rho}({c,p,W})\subset\mathcal{P}$. However, all but finitely many points $({c^{(k)},p^{(k)},W^{(k)}})$ belong to $\tilde{B}_{\rho}({c,p,W})$, and since $\hat{q}_{k}$ is bounded for all $({c^{(k)},p^{(k)},W^{(k)}})\in\tilde{B}_{\rho}({c,p,W})$, we cannot have $\|\hat{q}_{k}\|\to\infty$.

Next, let $({c^{(1)},p^{(1)},W^{(1)}}),({c^{(2)},p^{(2)},W^{(2)}})\in\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$ and let $\hat{q}_{1},\hat{q}_{2}\in\bar{\mathfrak{P}}_{+}$ be the unique minimizers of $\mathbb{J}_{{c^{(1)},p^{(1)},W^{(1)}}}$ and $\mathbb{J}_{{c^{(2)},p^{(2)},W^{(2)}}}$, respectively. Choose a $q_{0}\in\mathfrak{P}_{+}$ and note that $Q_{0}$ is strictly positive and bounded. By optimality,

[TABLE]

for all $\varepsilon>0$. Hence, if we can show that, for any $\delta>0$, there is an $\varepsilon>0$ and a $\tilde{\rho}>0$ such that

[TABLE]

hold whenever $\|c^{(1)}-c^{(2)}\|_{2}\leq\tilde{\rho}$, $\|p^{(1)}-p^{(2)}\|_{2}\leq\tilde{\rho}$ and $\|W^{(1)}-W^{(2)}\|_{F}\leq\tilde{\rho}$, then this would imply that

[TABLE]

showing that the optimal value is continuous in ${c^{(1)},p^{(1)},W^{(1)}}$. The lower bound is obtained by using (A.2a) and (A.1b), and the upper bound is obtained from (A.1a) and (A.2b). To prove (A.2a), we note that

[TABLE]

Next we observe that

[TABLE]

*by the KKT conditions (3.11) and the fact that $q_{0}\in\mathfrak{P}_{+}$, $\hat{c}_{1}\in\bar{\mathfrak{C}}_{+}$. Hence $\varepsilon$ can be selected small enough for the second and fourth term in (A.3) each to be bounded by $\delta/5$ for any $({c^{(1)},p^{(1)},W^{(1)}}),({c^{(2)},p^{(2)},W^{(2)}})\in B_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$. Each of the remaining terms can now be bounded by $\delta/5$ by selecting $\tilde{\rho}$ sufficiently small. Hence (A.2a) follows. This also proves (A.2b). *

Lemma A.3**.**

*Let $\mathcal{P}_{\rm hard}$ be given by (7.5) and $\mathbb{J}_{{c,p,W}}$ by (6.6). Furthermore, let $\hat{q}:=\min_{q\in\bar{\mathfrak{P}}_{+}}\mathbb{J}_{{c,p,W}}(q)$. Then for any compact $K\subset\mathcal{P}_{\rm hard}$, the corresponding set of optimal solutions $\hat{q}$ is bounded. *

Proof A.4**.**

The proof follows closely the proof of the corresponding part of Lemma A.1. Let $({c^{(0)},p^{(0)},W^{(0)}})\in\mathcal{P}_{\rm hard}$ be arbitrary and let

[TABLE]

where $\rho>0$ is chosen so that $\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})\subset\mathcal{P}_{\rm hard}$. To see that the minimizer $\hat{q}_{{c,p,W}}$ of $\mathbb{J}_{{c,p,W}}$ is bounded for all $({c,p,W})\in\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$, first note that by optimality

[TABLE]

and hence $\mathbb{J}_{c,p,W}(p)$ is bounded from above on the compact set $\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$.

Now let $h(q):=\langle c,q-e\rangle+\|q-e\|_{W}$, as in the proof of Theorem  6.1. Following the same line of argument as in that proof, we see that $h(q)>0$ for all $q\in\bar{\mathfrak{P}}_{+}$ and $({c,p,W})\in\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$. Since $h$ is continuous in the arguments $(q,{c,p,W})$, it has a minimum $\varepsilon>0$ on the compact set of tuples $(q,{c,p,W})$ such that $q\in\bar{\mathfrak{P}}_{+}\setminus\{0\},\|q-e\|_{\infty}=1,$ and $({c,p,W})\in\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$ hold. Thus the second half of inequality (6.14) still holds, i.e.,

[TABLE]

for all $q$. This is true in particular for $\hat{q}_{{c,p,W}}$, thus

[TABLE]

*Since the linear growth dominates the logarithmic growth, the norm of $\hat{q}_{{c,p,W}}$ is bounded on the set $\tilde{B}_{\rho}({c^{(0)},p^{(0)},W^{(0)}})$. The proof now follows verbatim from the argument in the second paragraph in the proof of Lemma A.1. *

Proof A.5** (Proof of Theorem 7.5).**

This is a modification of the proof of Theorem 7.1, again utilizing [11, Lemma 2.3], where we replace the map $f$ defined by $\mathcal{W}\ni W_{\rm hard}\mapsto W_{\rm soft}\in\{W\mid W>0\}$ and redefine it with the map ${\rm int}(\mathcal{P}_{\rm hard})\ni(c,p,W_{\rm hard})\mapsto(c,p,W_{\rm soft})\in{\rm int}(\mathcal{P})$. To show that $f$ is a homeomorphism we need to show that the map is proper. To this end, we take a compact set $K\subset{\rm int}(\mathcal{P}_{\rm hard})$ and show that $f^{-1}(K)$ is also compact. Again, there are two ways this could fail. First, the preimage could contain a singular semidefinite matrix. However this is impossible by (7.2), since $\|\hat{q}\|_{\infty}$ is bounded for $(c,p,W_{\text{hard}})\in K$ (Lemma A.3) and a nonzero scaling of a singular matrix cannot be nonsingular. Secondly, $\|W_{\rm soft}\|_{F}$ could tend to infinity. However, this is also impossible. To see this, we first show that there is a $\kappa>0$ such that $\|p-r\|_{W_{\rm hard}^{-1}}\geq\kappa$ for all $r\in\mathfrak{S}_{c,W_{\rm hard}}$ and all $(c,p,W_{\rm hard})\in K$. Again, using the triangle inequality $\|p-r\|_{W_{\rm hard}^{-1}}\geq\|p-c\|_{W_{\rm hard}^{-1}}-\|c-r\|_{W_{\rm hard}^{-1}}$, we observe that the minimum of $\|p-r\|_{W_{\rm hard}^{-1}}$ over all $(c,p,W_{\rm hard})\in K$ and $r$ satisfying the constraint $\|r-c\|_{W_{\rm hard}^{-1}}\leq 1$ is bounded by

[TABLE]

*The minimum is attained, as $K$ is compact, and positive, since $p\not\in\bigcup_{(c,W_{\rm hard})\in K}\mathfrak{S}_{c,W}$. The remaining part of the proof now follows with minor modifications from the proof of Theorem 7.1 by noting that $\hat{q}$ is bounded away from $e$, and hence the preimage $f^{-1}(K)$ is bounded. Therefore the limit of a sequence in the preimage must belong to $f^{-1}(K)$, and hence $f^{-1}(K)$ is compact as claimed. *