Synchronous deautoconvolution algorithm for discrete-time positive signals via I-divergence approximation

TL;DR

This paper introduces a new iterative algorithm for approximating nonnegative signals using autoconvolution, employing I-divergence as the criterion, with proven convergence properties and demonstrated effectiveness through numerical experiments.

Contribution

It develops a novel synchronous deautoconvolution algorithm based on I-divergence minimization and relaxation techniques, with proven convergence to Kuhn-Tucker points.

Findings

Algorithm converges to Kuhn-Tucker points.

Numerical experiments validate the algorithm's effectiveness.

The method effectively handles nonnegativity constraints.

Abstract

We pose the problem of the optimal approximation of a given nonnegative signal $y_t$ with the scalar autoconvolution $(x*x)_t$ of a nonnegative signal $x_t$, where $x_t$ and $y_t$ are signals of equal length. The $\mathcal{I}$-divergence has been adopted as optimality criterion, being well suited to incorporate nonnegativity constraints. To find a minimizer we derive an iterative descent algorithm of the alternating minimization type. The algorithm is based on the lifting of the original problem to a larger space, a relaxation technique developed by Csisz\'ar and Tusn\'ady in [Statistics \& Decisions (S1) (1984), 205--237] which, in the present context, requires the solution of a hard partial minimization problem. We study the asymptotic behavior of the algorithm exploiting the optimality properties of the partial minimization problems and prove, among other results, that its limit points are Kuhn-Tucker points of the original minimization problem. Numerical experiments illustrate the results.

Equations

$$(x\ast x{)}_i=\sum _{j=-\infty }^{\infty }{x}_{i-j}{x}_j=\sum _{j=0}^i{x}_{i-j}{x}_j,i\ge 0.$$

$$y=(y_0,\dots ,y_n),$$

$$\mathcal{I}(u,v)=\sum_{i}\Big{(}u_{i}\log\frac{u_{i}}{v_{i}}-u_{i}+v_{i}\Big{)}\,,$$

$$\mathcal{I}=\mathcal{I}(x):=\mathcal{I}(y||x*x)=\sum_{i=0}^{n}\Big{(}y_{i}\log\frac{y_{i}}{(x*x)_{i}}-y_{i}+(x*x)_{i}\Big{)}\,.$$

$$\nabla\!_{k}\mathcal{I}(x):=\frac{\partial\mathcal{I}(x)}{\partial x_{k}}=-2\sum_{j=0}^{n-k}\,\Big{(}\frac{y_{j+k}}{(x*x)_{j+k}}-1\Big{)}\,x_{j}.$$

$${\hat{y}}_i$$

$${\rho }_i$$

$$\nabla {}_k\mathcal{I}(x)$$

$$\sum _{j=k}^n{\rho }_j{x}_{j-k}=-\frac{1}{2}\nabla {}_k\mathcal{I}(x)+\sum _{j=k}^n{x}_{j-k}.$$

$$x_0$$

$$x_1$$

$$x_2$$

$$x_0$$

$$x_1$$

$$x_2$$

$$\frac{\partial\mathcal{I}}{\partial x_{2}}=\frac{\big{(}y_{0}+\frac{y_{1}}{2}\big{)}\big{(}\frac{y_{1}^{2}}{4}-y_{0}y_{2}\big{)}}{\big{(}y_{2}+\frac{y_{1}}{2}\big{)}^{2}\sqrt{y_{0}+y_{1}+y_{2}}}\geq 0\,.$$

$$\mbox{{\boldmath$\mathcal{Y}$}}:=\Big{\{}\,\mathbf{Y}:\,\,\,\textstyle{\sum_{j}}\mathbf{Y}_{ij}=y_{i}\,,\quad\text{and}\quad\mathbf{Y}_{ij}=0,\,\,\text{for}\,\,i<j\,\Big{\}}\,,$$

$$\mbox{{\boldmath$\mathcal{W}$}}:=\Big{\{}\,\mathbf{W}:\,\,\,\mathbf{W}_{ij}=x_{i-j}x_{j}\,,\,\,\text{if}\,\,j\leq i\,,\quad\text{and}\quad\mathbf{W}_{ij}=0\,,\,\,\text{for}\,\,i<j\,\Big{\}}\,,$$

$$\mathbf{Y}=\left(\begin{array}{ccccccc}{\mathbf{Y}}_{00}& 0& 0& 0& 0& 0& 0\\ {\mathbf{Y}}_{10}& {\mathbf{Y}}_{11}& 0& 0& 0& 0& 0\\ {\mathbf{Y}}_{20}& {\mathbf{Y}}_{21}& {\mathbf{Y}}_{22}& 0& 0& 0& 0\\ {\mathbf{Y}}_{30}& {\mathbf{Y}}_{31}& {\mathbf{Y}}_{32}& {\mathbf{Y}}_{33}& 0& 0& 0\\ {\mathbf{Y}}_{40}& {\mathbf{Y}}_{41}& {\mathbf{Y}}_{42}& {\mathbf{Y}}_{43}& {\mathbf{Y}}_{44}& 0& 0\\ {\mathbf{Y}}_{50}& {\mathbf{Y}}_{51}& {\mathbf{Y}}_{52}& {\mathbf{Y}}_{53}& {\mathbf{Y}}_{54}& {\mathbf{Y}}_{55}& 0\\ {\mathbf{Y}}_{60}& {\mathbf{Y}}_{61}& {\mathbf{Y}}_{62}& {\mathbf{Y}}_{63}& {\mathbf{Y}}_{64}& {\mathbf{Y}}_{65}& {\mathbf{Y}}_{66}\end{array}\right),$$

$$\mathbf{W}=\left(\begin{array}{ccccccc}{x}_0{x}_0& 0& 0& 0& 0& 0& 0\\ x_1{x}_0& x_0{x}_1& 0& 0& 0& 0& 0\\ x_2{x}_0& x_1{x}_1& x_0{x}_2& 0& 0& 0& 0\\ x_3{x}_0& x_2{x}_1& x_1{x}_2& x_0{x}_3& 0& 0& 0\\ x_4{x}_0& x_3{x}_1& x_2{x}_2& x_1{x}_3& x_0{x}_4& 0& 0\\ x_5{x}_0& x_4{x}_1& x_3{x}_2& x_2{x}_3& x_1{x}_4& x_0{x}_5& 0\\ x_6{x}_0& x_5{x}_1& x_4{x}_2& x_3{x}_3& x_2{x}_4& x_1{x}_5& x_0{x}_6\end{array}\right).$$

$$\mathcal{I}(M||N):=\sum_{i,j}\Big{(}M_{ij}\log\frac{M_{ij}}{N_{ij}}-M_{ij}+N_{ij}\Big{)}\,.$$

$$\min_{x\in\mathbb{R}^{n+1}_{+}}\mathcal{I}(x)=\min_{\mathbf{Y}\in\mbox{{\boldmath$\mathcal{Y}$}}}\min_{\mathbf{W}\in\mbox{{\boldmath$\mathcal{W}$}}}\mathcal{I}(\mathbf{Y}||\mathbf{W}),$$

$${\mathbf{Y}}_{ij}^{\star }$$

$$=\frac{x_j{x}_{i-j}}{{\sum }_j{x}_j{x}_{i-j}}{y}_i=\frac{x_j{x}_{i-j}}{(x\ast x{)}_i}{y}_i,$$

$$\mathcal{I}(\mathbf{Y}||\mathbf{W})=\mathcal{I}(\mathbf{Y}||{\mathbf{Y}}^{\star })+\mathcal{I}({\mathbf{Y}}^{\star }||\mathbf{W}),$$

$${\tilde{\mathbf{Y}}}_j=\sum _{i=j}^n{\mathbf{Y}}_{ij}+\sum _{i=j}^n{\mathbf{Y}}_{i,i-j}.$$

$$2x_j\sum _{i=0}^{n-j}{x}_i={\tilde{\mathbf{Y}}}_j,j=0,\dots ,n.$$

$$\mathcal{I}(\mathbf{Y}||\mathbf{W})=\mathcal{I}(\mathbf{Y}||{\mathbf{W}}^{\star })+\mathcal{I}({\mathbf{W}}^{\star }||\mathbf{W}).$$

$$\sum _{i=0}^n(x^{\star }\ast x^{\star }{)}_i=\sum _{ij}{\mathbf{Y}}_{ij}.$$

$$\mathcal{I}(\mathbf{Y}||\mathbf{W})=\sum _{ij}{\mathbf{Y}}_{ij}\log \frac{{\mathbf{Y}}_{ij}}{{\mathbf{W}}_{ij}}-\sum _{ij}{\mathbf{Y}}_{ij}+\sum _{ij}{\mathbf{W}}_{ij}.$$

$$F(\mathbf{W})$$

$$=-\sum _{i=0}^n\sum _{j=0}^i{\mathbf{Y}}_{ij}\log (x_j{x}_{i-j})+\sum _{i=0}^n\sum _{j=0}^i{x}_j{x}_{i-j}$$

$$=-\sum _{i=0}^n\sum _{j=0}^i{\mathbf{Y}}_{ij}\log x_j-\sum _{i=0}^n\sum _{j=0}^i{\mathbf{Y}}_{i,i-j}\log x_j+\sum _{i=0}^n\sum _{j=0}^i{x}_j{x}_{i-j}$$

$$\displaystyle=-\sum_{j=0}^{n}\Big{(}\sum_{i=j}^{n}\mathbf{Y}_{ij}\Big{)}\log x_{j}-\sum_{j=0}^{n}\Big{(}\sum_{i=j}^{n}\mathbf{Y}_{i,i-j}\Big{)}\log x_{j}+\sum_{i=0}^{n}\sum_{j=0}^{i}x_{j}x_{i-j}$$

$$=-\sum _{j=0}^n{\tilde{\mathbf{Y}}}_j\log x_j+\sum _{i=0}^n\sum _{j=0}^i{x}_j{x}_{i-j},$$

Full text

Keywords: autoconvolution, inverse problem, positive system, I-divergence, alternating minimization

AMS subject classification: 93B30, 94A17

1 Introduction

The focus of the paper is on the inverse problem of positive autoconvolution. Given the data sequence of nonnegative observations $y_{t}$, $t$ in an appropriate time domain, find a nonnegative signal $x_{t}$ such that the autoconvolutions $(x*x)_{t}:=\sum_{i=0}^{t}x_{i}x_{t-i}$ form an optimal approximation of the $y_{t}$. We call such $x_{t}$ the positive deautoconvolution of $y_{t}$. The problem of finding the deautoconvolution can be posed in two different setups. The first is to approximate the data sequence $y_{t}$, $t=0,\ldots,n$, with the autoconvolutions $(x*x)_{t}$ of a finite support signal $x_{t}$, $t=0,\ldots,m$. Since the support of $(x*x)_{t}$ is $[0,2m]$, in this setup $n=2m$, i.e. the support of the deautoconvolution is half of the support of $y_{t}$, hence we call this setup the half support case. An algorithm for finding the deautoconvolution $x_{t}$ in the half support case has been proposed, and its properties have been discussed, in the companion paper [10]. The second approach is to approximate the data sequence $y_{t}$, $t=0,\ldots n$ with the initial segment of the autoconvolution $(x*x)_{t}$, $t=0,\dots,n$ of a signal $x_{t}$ of long, unspecified length. In this case the problem is to find the initial segment of the deautoconvolution, i.e. the signal $x_{t}$, $t=0,\ldots,n$. In this approach the data sequence $y_{t}$ and the deautoconvolution $x_{t}$ have common support, $t=0,\ldots,n$, hence this is called the full support or synchronous case. The present paper deals with the synchronous setup. Despite obvious similarities between the two setups, it turns out that the algorithm to construct the deautoconvolution is considerably more involved in the synchronous setup, and the results strongly differ between the two cases. A similar situation, with two alternative setups, arises also for the simpler problem of nonnegative deconvolution, in which one approximates optimally the given data sequence $y_{t}$ with the convolution $(h*x)_{t}$, where $h_{t}$ is known and $x_{t}$ is to be determined ($y_{t}$, $h_{t}$ and $x_{t}$ being all nonnegative). Compare with [8], [9].

Extensive work has been dedicated to the inverse problem of autoconvolution for functions on the real line, mostly emphasizing the functional analytic aspects and motivating its interest in a variety of applications in physics and engineering. We briefly mention in this context the references [5], [12], [11], [6]. The main differences between the cited literature for deconvolution of functions and this paper are that we consider approximation problems, rather than looking for exact solutions which exist only exceptionally, and that our (time) domain is discrete rather than the real line. Moreover, the nonnegativity constraint, that we impose on all signals, is a crucial feature of the present work. Some earlier contributions share, at least in part, our point of view, e.g. the papers [1], [2], contain an algorithm of the same type as ours, but valid in the half support case. For more detailed comments on  [1], [2], see [10].

The purpose of this paper is threefold. First we pose the problem of time domain approximation of a given discrete time nonnegative data sequence with finite autoconvolutions. Following the choice made in other optimization problems, we opt for the I-divergence as a criterion, which as argued in [4] (see also [13]), is the natural choice under nonnegativity constraints. We provide a result on the existence of the minimizer of the approximation criterion. Then we propose an iterative algorithm to find the best approximation, and present a number of results on its behavior. Our approach – alternating minimization – is based on earlier works like [7] on nonnegative matrix factorization, [8], [9], on linear convolutional problems, in contrast to the inherently nonlinear autoconvolution.

The inherent nonconvexity, and nonlinearity of the problem make the analysis of the asymptotic behavior difficult, even more so because of the existence of several local minima of the objective function and possibly identification issues as well. The main result in this respect is Theorem 4 which states that limit points of the algorithm satisfy the Kuhn-Tucker optimality conditions, which is, in formulation, similar to the corresponding result in [10]. In the contest of nonlinear non convex optimization problems, like the present one, to the best of our knowledge there are no asymptotic results that improve on our Theorem 4.

A brief summary of the paper follows. In Section 2 we state the full support minimization problem, show the existence of a solution and give some of its properties. In Section 3 the original problem is lifted into a higher dimensional setting, thus making it amenable to alternating minimization. The optimality properties (Pythagoras rules) of the ensuing partial minimization problems are discussed here. It turns out that one of these partial minimization problems has a very complex solution, whose precise representation will be derived in the appendix. In Section 4 we derive the iterative minimization algorithm combining the solutions of the partial minimizations, and analyse its properties. In particular we show that limit points of the algorithm are Kuhn-Tucker points of the original optimization problem. In the concluding Section 5 we present numerical experiments that shed light on the theoretical results concerning the behaviour of the algorithm and also highlight some unexpected phenomena.

2 Problem statement and initial results

In the paper we consider real valued signals $x:\mathbb{Z}\rightarrow\mathbb{R}$, maps $i\mapsto x_{i}$, that vanish for $i<0$, i.e., $x_{i}=0$ for $i<0$. The support of $x$ is the discrete time interval, $[0,n]$, where n=\inf\{\,k:\,\,x_{i}=0,\,\,\,\text{for i>k}\,\}, if the infimum is finite (and then a minimum), or $[0,\infty)$ otherwise. The autoconvolution of $x$ is the signal $x*x$, vanishing for $i<0$, and satisfying,

[TABLE]

Notice that if the support of $x$ is finite $[0,n]$, the support of $x*x$ is $[0,2n]$. In this case, when computing $(x*x)_{i}$ for $i>n$, the summation in Equation (1) has non zero addends only in the range $i-n\leq j\leq n$, as $x_{i-j}=0$ and $x_{j}=0$ for $i-j>n$ and $j>n$ respectively. If the signal $x$ is nonnegative, i.e. $x_{i}\geq 0$ for all $i\in\mathbb{Z}$, the autoconvolution (1) is too.

The paper addresses the inverse problem of autoconvolution. Given a finite nonnegative data sequence

[TABLE]

find a nonnegative signal $x$ whose autoconvolution $x*x$ best approximates $y$. Since the signals involved are nonnegative, it is natural to adopt the I-divergence as the approximation criterion, see [3, 4]. Recall that the I-divergence between two nonnegative vectors $u$ and $v$ of equal length is

[TABLE]

if $u_{i}=0$ whenever $v_{i}=0$, and $\mathcal{I}(u,v)=\infty$ if there exists an index $i$ with $u_{i}>0$ and $v_{i}=0$. It is known that $\mathcal{I}(u,v)\geq 0$, with equality iff $u=v$.

The approximation problem depends on the constraints imposed on the support of $x$. In Problem 2 below, the signal $x$ is assumed to have the same length of the data sequence, $x=(x_{0},\dots,x_{n})$ and produces the approximation problem specified below, where we write $x*x\in\mathbb{R}^{n+1}$ for the restriction to $[0,n]$ of the convolution defined in (1). This version of the problem arises naturally when the observed data sequence is the initial segment, $y=(y_{0},\ldots,y_{n})$, of a signal $y$ of long unspecified length, modeled with the initial segment of the autoconvolution $x*x$ of a signal $x$ of long unspecified length.

• Problem 2.1

Given $y\in\mathbb{R}^{n+1}_{+}$ minimize, over $x\in\mathbb{R}^{n+1}_{+}=[0,\infty)^{n+1}$,

[TABLE]

The objective function (2) is non convex and nonlinear in $x$, the existence of a minimizer is therefore not immediately clear. Our first result settles in the affirmative the question of the existence. The issue of uniqueness remains open, but we have evidence of the existence of multiple local minima of $\mathcal{I}(x)$. See Section 5 for numerical examples.

• Proposition 2.2

Let $y_{0}>0$, then Problem 2 admits a solution. If all $y_{k}$ are strictly positive, then the minimum is uniquely attained.

• Proof

Deferred to Section 4 as we need properties of the algorithm of that section. $\Box$

For later reference we compute the components of the gradient $\nabla\mathcal{I}(x)$ of the objective function (2).

[TABLE]

We will also often use the notations

[TABLE]

and as a result we have the following alternative ways of writing Equation (3),

[TABLE]

and then also

[TABLE]

The expressions in (3) are highly nonlinear in $x$ and solving the first order optimality conditions $\nabla\mathcal{I}(x)=0$ to find the stationary points of (2), will not result in closed form solutions except in trivial cases. This observation calls for a numerical approach to the optimization, which we will present in Section 4. First a useful property of minimizers of $\mathcal{I}$.

• Proposition 2.3

Minimizers $x^{\star}$ are such that $\sum_{j=0}^{n}(x^{\star}*x^{\star})_{j}=\sum_{j=0}^{n}y_{j}$.

• Proof

Suppose that $x^{\star}$ is a minimizer of $\mathcal{I}(x)$ and let $f(\alpha)=\mathcal{I}(\alpha x^{\star})$, $\alpha>0$. Then $f^{\prime}(1)=0$. A direct computation of $f(\alpha)$ gives after differentiation $f^{\prime}(\alpha)=-\frac{2}{\alpha}\sum_{i=0}^{n}y_{i}+2\alpha\sum_{i=0}(x^{\star}*x^{\star})_{i}$. Hence $f^{\prime}(1)=0$ yields the result. $\Box$

Let us call $x$ an interior solution of the minimization of $\mathcal{I}(x)$ if $\nabla\mathcal{I}(x)=0$, where in a point $x\in\mathbb{R}^{n+1}_{+}$ with some coordinate $x_{j}=0$ partial derivatives are understood as right derivatives. A boundary solution occurs at a point $x$ with some $x_{j}=0$ and with a partial derivative $\nabla\!_{j}\mathcal{I}(x)$ that is strictly positive.

We now give some examples, where the minimizer can be explicitly calculated. These are exceptional cases.

• Example 2.4

Let $n=1$. We will look at four cases.

(i) $y_{0}=y_{1}=0$. Here $\mathcal{I}(x)=x_{0}^{2}+2x_{0}x_{1}$ and we observe that $\mathcal{I}$ is minimized for $x_{0}=0$ and arbitrary $x_{1}$. The minimum is not uniquely attained.

(ii) $y_{0}=0<y_{1}$. Then $\mathcal{I}(x)=y_{1}\log\frac{y_{1}}{2x_{0}x_{1}}-y_{1}+x_{0}^{2}+2x_{0}x_{1}$. Let $x_{0},x_{1}$ be such that $2x_{0}x_{1}=y_{1}$, then $\mathcal{I}(x)=x_{0}^{2}$ and we see that the infimum of $\mathcal{I}(x)$ is not attained.

(iii) $y_{1}=0<y_{0}$. Then $\mathcal{I}(x)=y_{0}\log\frac{y_{0}}{x_{0}^{2}}-y_{0}+x_{0}^{2}+2x_{0}x_{1}$. Here one shows that the minimium is uniquely attained for a boundary point $x$ with $x_{0}=\sqrt{y_{0}}$ and $x_{1}=0$, which gives $\mathcal{I}(x)=0$. At the minimum we have that the gradient $\nabla\mathcal{I}(x)=(0,2x_{0})$.

(iv) $y_{0},y_{1}>0$. The minimum is uniquely attained at an interior point $x$ with $x_{0}=\sqrt{y_{0}}$ and $x_{1}=\frac{y_{1}}{2\sqrt{y_{0}}}$. At the minimum we have $\nabla\mathcal{I}(x)=(0,0)$.

The message of this example is that $y_{0}=0$ (cases (i) and (ii)) displays undesired properties when searching for a minimum, and that minimizers at the boundary of the domain occur.

• Example 2.5

Let $n=2$ and $y_{0}>0$. Then it depends on the specific values of the $y_{k}$ whether the problem has an interior solution. This happens if $y_{1}^{2}<4y_{0}y_{2}$, which gives the minimizing

[TABLE]

and then the divergence equals zero. If $y_{1}^{2}\geq 4y_{0}y_{2}$, then the solution is quite different. In this case the minimizers become

[TABLE]

Of course the two formulas coincide in the boundary case $y_{1}^{2}=4y_{0}y_{2}$, which implies $y_{0}+y_{1}+y_{2}=(\sqrt{y_{0}}+\sqrt{y_{2}})^{2}$, as they should. Note that in the second case the partial derivatives $\frac{\partial\mathcal{I}}{\partial x_{j}}$ are zero at the minimizer for $j=0,1$, whereas

[TABLE]

The next proposition shows that internal solutions to Problem 2 are exceptions rather than the rule.

• Proposition 2.6

Assume $y_{0}>0$. A point $x\in\mathbb{R}^{n+1}_{+}$ is an internal solution iff it gives a perfect match with the $y_{k}$, i.e. $(x*x)_{k}=y_{k}$ for all $k=0,\ldots,n$.

• Proof

Consider Equation (3). From this equation one immediately sees that an exact match produces an internal solution. To show the converse, note first that $x_{0}>0$ and let $\nabla\mathcal{I}(x)=0$. Then, by virtue of (3), for all $k=0,\ldots,n$ one has \sum_{j=0}^{n-k}\big{(}\frac{y_{j+k}}{(x*x)_{j+k}}-1\big{)}x_{j}=0. For $k=n$ this equation reads $(\frac{y_{n}}{(x*x)_{n}}-1)x_{0}=0$ and hence $y_{n}=(x*x)_{n}$. For $k=n-1$, the equation becomes $(\frac{y_{n-1}}{(x*x)_{n-1}}-1)x_{0}+(\frac{y_{n}}{(x*x)_{n}}-1)x_{1}=0$, from which we can now deduce $y_{n-1}=(x*x)_{n-1}$. Counting $k$ down to zero yields the assertion. $\Box$

3 Lifting

Following the approach of the companion paper [10] we recast Problem 2 as a double minimization problem by lifting it into a larger space, a necessary step to derive the minimization algorithm. Related work with similar techniques is [7] to analyse a nonnegative matrix factorization problem, and closer to the present paper [8], [9], noting that the latter references treat linear convolutional problems, whereas the autoconvolution is inherently nonlinear.

The ambient spaces for the lifted problem are the subsets $\mathcal{Y}$ and $\mathcal{W}$, defined below, of the set of matrices $\mathbb{R}^{(n+1)\times(n+1)}_{+}$,

[TABLE]

where $y=(y_{0},\dots,y_{n})\in\mathbb{R}^{n+1}_{+}$ is the given data vector, and

[TABLE]

where $x=(x_{0},\ldots x_{n})$ are nonnegative parameters. The structure of the matrices in $\mathcal{Y}$ and $\mathcal{W}$ is shown below for $n=6$,

[TABLE]

[TABLE]

The interpretation is as follows. The matrices \mathbf{Y}\in\mbox{{\boldmath\mathcal{Y}}} and \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}} have common support on (sub)diagonals. The row marginal (i.e. the column vector of row sums) of any \mathbf{Y}\in\mbox{{\boldmath\mathcal{Y}}} coincides with the given data vector $y$. The elements of the $\mathbf{W}$ matrices factorize and their row marginals are the autoconvolutions $(x*x)_{i}$, $i=0,\ldots,n$. Note also that the elements \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}} exhibit the symmetry $\mathbf{W}_{ij}=\mathbf{W}_{i,i-j}$, $i\geq j$.

We introduce two partial minimization problems over the subsets $\mathcal{Y}$ and $\mathcal{W}$. Recall that the I-divergence between two nonnegative matrices of the same sizes $M,N\in\mathbb{R}^{p\times q}_{+}$ is defined as

[TABLE]

• Problem 3.1

Minimize $\mathcal{I}(\mathbf{Y}||\mathbf{W})$ over \mathbf{Y}\in\mbox{{\boldmath\mathcal{Y}}} for given \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}}.

• Problem 3.2

Minimize $\mathcal{I}(\mathbf{Y}||\mathbf{W})$ over \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}} for given \mathbf{Y}\in\mbox{{\boldmath\mathcal{Y}}}.

The relation between the original Problem 2 and the partial minimization Problems 3 and 3 is as follows.

• Proposition 3.3

It holds that

[TABLE]

where in all three cases the minima are attained.

• Proof

The proof is almost the same as that of Proposition 3.8 in [10]. $\Box$

The partial minimization problems are instrumental in deriving the algorithm in Section 4 for which one needs that these two problems have solutions in closed form. For Problem 3 it can easily be given, Lemma 3. On the other hand, the solution to Problem 3, see Lemma 3 that relies on results in Appendix A, turns out to be substantially more complicated and therefore radically differs from the solution to the companion Problem 3.2 in [10]. The remainder of this section is devoted to properties of the solutions to the partial minimization problems.

• Lemma 3.4

Problem 3 has the explicit minimizer \mathbf{Y}^{\star}\in\mbox{{\boldmath\mathcal{Y}}} given by

[TABLE]

for \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}} and $i\geq j$. Moreover, $\mathcal{I}(\mathbf{Y}^{\star}||\mathbf{W})=\mathcal{I}(y||x*x)$ and the Pythagorean identity

[TABLE]

holds for any \mathbf{Y}\in\mbox{{\boldmath\mathcal{Y}}}.

• Proof

By direct computation and almost verbatim the same as the proof of Lemma 3.3 in [10]. $\Box$

• Remark 3.5

The optimal $\mathbf{Y}^{\star}$ in Lemma 3 exhibits the symmetry $\mathbf{Y}^{\star}_{ij}=\mathbf{Y}^{\star}_{i,i-j}$ for all $i\geq j$, which is the same symmetry as enjoyed by the \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}}.

• Lemma 3.6

Consider Problem 3 for given \mathbf{Y}\in\mbox{{\boldmath\mathcal{Y}}} and let for $j=0,\ldots,n$,

[TABLE]

The minimizing $x_{i}^{\star}$ uniquely exist and are the solutions to the nonlinear system of equations

[TABLE]

The Pythagorean identity holds for the solutions, i.e. for every \mathbf{W}\in\mbox{{\boldmath\mathcal{W}}} it holds that

[TABLE]

Furthermore, it holds that

[TABLE]

• Proof

We compute

[TABLE]

It follows that we have to minimize

[TABLE]

We compute the gradient of $F$,

[TABLE]

Setting $\frac{\partial F}{\partial x_{j}}=0$ for all $j$ gives equations (9).

Next we show that this system of equations has a unique solution. Let

[TABLE]

and $k=\lfloor\frac{n}{2}\rfloor$. In this new notation we have to solve

[TABLE]

We have to verify Condition (22) of Lemma A for the $r_{j}$ given by (11). The result then follows by application of Propositions A.1 or A.2. We write

[TABLE]

Consider the first term in $S^{2}$,

[TABLE]

and of it the first sum only. We split

[TABLE]

Here the first term on the right is our goal. We distinguish two cases, $n$ odd, $n=2k+1$ and $n$ is even, $n=2k$. Consider first $n$ is odd. The second term on the right in (14) becomes

[TABLE]

Next we work on the second sum in (13). We develop with $n=2k+1$

[TABLE]

which coincides with (15). This finishes the case $n=2k+1$ and we proceed with $n=2k$. The second term on the right in (14) now becomes

[TABLE]

Next we work on the second sum in (13). We develop with $n=2k$

[TABLE]

which seems different from (16), the difference being the summand when $j=k$. If we take $j=k$ in (16), we get for the sum over $i$ the empty sum $\sum_{i=2k+1}^{2k}\mathbf{Y}_{ij}$, which is zero. Hence the final term of the last display indeed coincides with (16), as desired. This finishes the case $n=2k$ for the first term in parentheses for $S^{2}$. The proof for the second term in parentheses in $S^{2}$ is similar.

To prove the Pythagorean identity (10), we compute both the quantities $\mathcal{I}(\mathbf{Y}||\mathbf{W})-\mathcal{I}(\mathbf{Y}||\mathbf{W}^{\star})$ and $\mathcal{I}(\mathbf{W}^{\star}||\mathbf{W})$ and show that they are the same. We put $\widetilde{\mathcal{I}}(\mathbf{Y}||\mathbf{W})=\sum_{i}\sum_{j\leq i}\mathbf{Y}_{ij}\log\frac{\mathbf{Y}_{ij}}{\mathbf{W}_{ij}}$ and a similar convention for the other modified divergences. It is then sufficient to show that $\widetilde{\mathcal{I}}(\mathbf{Y}||\mathbf{W})-\widetilde{\mathcal{I}}(\mathbf{Y}||\mathbf{W}^{\star})$ and $\widetilde{\mathcal{I}}(\mathbf{W}^{\star}||\mathbf{W})$ are the same. Recall $\mathbf{W}_{ij}=x_{j}x_{i-j}$ for $i\geq j$ and zero otherwise. For the elements $\mathbf{W}^{\star}_{ij}$ one has similar expressions. A direct computation now gives

[TABLE]

Similar computations lead to

[TABLE]

Since the $x_{j}^{\star}$ solve the system of equations (9), the term in parentheses equals $\widetilde{\mathbf{Y}}_{j}$ and we are done with proving the Pythagorean identity. As the optimal $x^{\star}$ satisfies (9), we can sum the equalities over all $j$ to get by swapping the summation order $2\sum_{i=0}^{n}(x^{\star}*x^{\star})=\sum_{ij}\tilde{\mathbf{Y}}_{ij}$. But $\sum_{ij}\tilde{\mathbf{Y}}_{ij}$ just equals $2\sum_{ij}\mathbf{Y}_{ij}$ by the special structure of the matrix $\mathbf{Y}$. The last claim follows. $\Box$

• Remark 3.7

Under the symmetry condition on $\mathbf{Y}$ as in Remark 3, Equation (9) reduces to

[TABLE]

This means that, for every $j$, the sum of the $j$-th column of the matrix $\mathbf{W}$ equals the sum of the $j$-th column of the matrix $\mathbf{Y}$ and, by symmetry, that the sum of the corresponding subdiagonal of the matrix $\mathbf{W}$ equals that of the sum of its counterpart of the matrix $\mathbf{Y}$.

4 The algorithm

This section contains the key results of the paper. We derive an optimal algorithm producing iterates $x^{t}\in\mathbb{R}^{n+1}_{+}$, $t\geq 0$, aiming at approximating a solution to Problem 2. Results of Appendix A are used throughout.

Start with a vector $x^{0}=(x^{0}_{0},\ldots,x^{0}_{n})$ with positive elements. An alternating concatenation of the solutions to Problems 3 and 3 yields the following. All relevant quantities of the algorithm have an index $t$ at the $t$-th step. Suppose the $t$-iterate $x^{t}$ is given, and then also $\mathbf{W}^{t}$. According to Lemma 3, with $\mathbf{W}^{t}$ as input, and recalling (4), (5), the optimal $\mathbf{Y}^{t}$ is given by

[TABLE]

In the next step, one uses the $\mathbf{Y}^{t}$ as input to compute according to Lemma 3 the optimal $x^{t+1}$ and $\mathbf{W}^{t+1}$. So, one has to solve the analog of (9), i.e. finding the $x^{t+1}_{j}$ as the solutions to

[TABLE]

or, in view of Remark 3, the solutions to

[TABLE]

The $x^{t+1}_{j}$ can be taken as the $x^{\star}_{j}$ in Proposition A.1 for even $n$, or Proposition A.2 for odd $n$. Although explicit expressions are complicated, the above procedure is amenable to a practically implementable algorithm. The following algorithm is the result of this procedure.

Algorithm 4.1**.**

The data are given as $(y_{0},\ldots,y_{n})$. Initiate the algorithm at a point $x^{0}\in(0,\infty)^{n+1}$ and let $x^{t}=(x^{t}_{0},\ldots,x^{t}_{n})$ be the values at the $t$-th iteration of the algorithm. To compute the values $x^{t+1}_{i}$ at the $t+1$-th iteration:

Compute the $\widehat{y}^{t}_{j}$ and $\rho^{t}_{j}$ as (see (4), (5))

[TABLE] 2. 2.

Compute $r^{t}_{j}$ from the $x^{t}_{j}$ according to

[TABLE] 3. 3.

Compute the $B^{t}_{j}$ and $E^{t}_{j}$ from the $r^{t}_{j}$ according to (23) and (24). 4. 4.

Use Proposition A.1 ($n$ even) or Proposition A.2 ($n$ odd) to compute the $x^{t+1}_{m}$ as the $x^{\star}_{m}$ in those propositions from the $B^{t}_{j}$ and $E^{t}_{j}$.

Here we present an alternative description of the algorithm, only for the case where $n$ is even, $n=2k$ to get more insight in the way Algorithm 4.1 is structured. This also forms the basis of a practical implementation. Suppose all $x^{t}_{j}$ and related quantities are given at the $t$-th iteration, e.g. $S^{t}=\sum_{j=0}^{k}x^{t}_{j}$, as well as $x^{t+1}_{k-l},\ldots,x^{t+1}_{k+l}$ for some $l\in\{0,\ldots,k-1\}$. Then one computes $x^{t+1}_{k+l+1}$ recursively (in $l$) as the solution to

[TABLE]

and using this value, one computes $x^{t+1}_{k-l-1}$ from

[TABLE]

Hence starting from the previous iterate $x^{t}$ and the ‘middle’ value $x^{t+1}_{k}$ (the computation of which follows from Proposition A.1), the above two equations produce the other $x^{t+1}_{j}$.

We proceed with some properties of Algorithm 4.1. The portmanteau Proposition 4 below summarizes some useful properties of the algorithm, in particular it quantifies the update gain at each iteration step.

• Proposition 4.2

The iterates $x^{t},\,t\geq 0$, of Algorithm 4.1 satisfy the following properties.

1. (i)

If $x^{0}>0$ componentwise, then $x^{t}>0$ componentwise, for all $t>0$. 2. (ii)

The iterates $x^{t}$ belong to the set $\mathcal{S}=\{x\in\mathbb{R}^{n+1}_{+}:\sum_{i=0}^{n}(x*x)_{i}=\sum_{i=0}^{n}y_{i}\}$ for all $t>0$, in agreement with Proposition 2. Additionally, the iterates $x^{t}$ also satisfy $\sum_{j=0}^{n}x^{t}_{j}\nabla_{j}\mathcal{I}(x^{t})=0$ for $t\geq 1$. 3. (iii)

$\mathcal{I}(y||x^{t}*x^{t})$ decreases at each iteration, in fact one has

[TABLE]

where one can use the expression

[TABLE]

As a corollary, $\mathcal{I}(\mathbf{W}^{t+1}||\mathbf{W}^{t})$ vanishes asymptotically. 4. (iv)

If $y=x^{t}*x^{t}$ then $x^{t+1}=x^{t}$, i.e. perfect matches are fixed points of the algorithm. 5. (v)

The update relation can be written in a formula –albeit an implicit expression– in terms of the gradient,

[TABLE] 6. (vi)

If $\nabla\mathcal{I}(x^{t})=0$ then $x^{t+1}=x^{t}$, i.e. stationary points of $\mathcal{I}(x)$ are fixed points of the algorithm. Additionally, if $x^{t}_{j}=0$, then also $x^{t+1}_{j}=0$.

• Proof

(i) is obvious.

(ii) Recall that $x^{t+1}$ is for $t\geq 0$ the solution to (18). Summing over $j$ one obtains $\sum_{j=0}^{n}x^{t+1}_{j}\sum_{i=0}^{n-j}x^{t+1}_{i}=\sum_{i=0}^{n}y_{i}$ as the row sums of $\mathbf{Y}^{t}$ are the $y_{i}$. By swapping the summation order, one sees $\sum_{j=0}^{n}x^{t+1}_{j}\sum_{i=0}^{n-j}x^{t+1}_{i}=\sum_{i=0}^{n}(x^{t+1}*x^{t+1})_{i}$. The first result follows. Using this in (20) one obtains $\sum_{j=0}^{n}x^{t}_{j}\nabla_{j}\mathcal{I}(x^{t})=0$ for $t\geq 1$, as desired.

(iii) Recall that $x^{t+1}_{j}\sum_{i=j}^{n}x^{t+1}_{i-j}=r^{t}_{j}$ for all $j$, which follows from Equation (12) with $r=r^{t}$ and $x=x^{t+1}$. Property (ii), the iterates $x^{t}$ belong to the set $\mathcal{S}$, allows us to write, also using $r^{t}_{j}$ as in (18),

[TABLE]

It follows from (19) that the sum $\sum_{t=0}\mathcal{I}(y||x^{t}*x^{t})-\mathcal{I}(y||x^{t+1}*x^{t+1})$ is finite and so is then $\sum_{t=0}^{\infty}\mathcal{I}(\mathbf{W}^{t+1}||\mathbf{W}^{t})$, which implies $\mathcal{I}(\mathbf{W}^{t+1}||\mathbf{W}^{t})\to 0$ for $t\to\infty$.

(iv) In this case one has $\rho_{j}=1$ for all $j$. Hence, (18) reads

[TABLE]

As this system of equations has a unique solution in terms of the variables $x^{t+1}_{j}$ by virtue of Propositions A.1 and A.2, the solution is necessarily $x^{t+1}_{j}=x^{t}_{j}$ for all $j$.

(v) In view of (6) the sum $\sum_{i=j}^{n}\mathbf{Y}^{t}_{ij}=x^{t}_{j}\sum_{i=j}^{n}x^{t}_{i-j}\rho_{i}^{t}$ in (18) can be written as x^{t}_{j}\big{(}-\frac{1}{2}\nabla_{j}\mathcal{I}(x^{t})+\sum_{i=j}^{n}x^{t}_{i-j}\big{)}, hence the $x^{t+1}_{j}$ also satisfy (20).

(vi) If some $x^{t}_{m}=0$, then by the multiplicative nature of the algorithm (see (20)), also $x^{t+1}_{m}=0$, so in the further analysis we assume that all $x^{t}_{m}>0$. By Propositions A.1 and A.2 the $x^{t+1}_{j}$ are unique solutions to the system of equations. If $\nabla\mathcal{I}(x^{t})=0$, then the equations to solve become

[TABLE]

and therefore $x^{t+1}=x^{t}$. $\Box$

Knowing properties of Algorithm 4.1 we now present a proof of Proposition 2.

• Proof

of Proposition 2. As $y_{0}>0$, one can only have finite divergence if $x_{0}>0$, because of the term $y_{0}\log\frac{y_{0}}{(x*x)_{0}}$ in the divergence. Furthermore, $\sum_{k=0}^{n}(x*x)_{k}\geq x_{0}\sum_{k=0}^{n}x_{k}$. Hence if some $x_{k}\to\infty$, the divergence tends to infinity as well. Therefore, the minimum has to be sought in a sufficiently large compact $K_{0}$ set of $\mathbb{R}^{n+1}_{+}$. We use Proposition 2 and various results from Proposition 4. Let $x=x^{0}$ be an arbitrary point in $\mathbb{R}^{n+1}_{+}$. Then one step of Algorithm 4.1 yields $x^{1}\in\mathbb{R}^{n+1}_{+}$ satisfying $\mathcal{I}(x^{1})\leq\mathcal{I}(x)$. Hence we can limit $K_{0}$ to those $x$ which also satisfy $\sum_{j=0}^{n}(x*x)_{j}=\sum_{j=0}^{n}y_{j}$, which gives a compact set as well, call it $K_{1}$. Note that $\mathcal{I}(x)=\sum_{k:y_{k}>0}(y_{k}\log\frac{y_{k}}{(x*x)_{k}}-y_{k})+\sum_{k}(x*x)_{k}$. To find the minimum, we can restrict ourselves even further to those $x$ for which $(x*x)_{k}\geq\varepsilon$ for all $k$ such that $y_{k}>0$, by choosing $\varepsilon$ sufficiently small and positive. This implies that we restrict the finding of the minimizers to an even smaller compact set $K_{2}$ on which $\mathcal{I}$ is continuous. This proves the existence of a minimizer. $\Box$

The Algorithm 4.1 shares some properties with its counterpart in [10]. We mention only a few related to its large running time behavior. We follow ideas taken from [8] and [10].

• Lemma 4.3

Let all data $y_{i}$ be positive. Limit points of the sequence $(x^{t})$ are fixed points of Algorithm 4.1. Moreover, if along a subsequence $(t_{j})$ one has the convergence $x^{t_{j}}\to x^{\infty}$, then also $x^{t_{j}+1}\to x^{\infty}$.

• Proof

Let $\mathcal{I}^{t}=\mathcal{I}(y||x^{t}*x^{t})$. Recall from Proposition 4(iii) that the divergence $\mathcal{I}(\mathbf{W}^{t}||\mathbf{W}^{t+1})\to 0$. Let $x^{\infty}$ be a limit point of the algorithm and $x^{\infty+1}$ its next iterate, with corresponding $\mathbf{W}^{\infty}$ and $\mathbf{W}^{\infty+1}$, then by continuity of the divergence in both arguments $\mathcal{I}(\mathbf{W}^{\infty}||\mathbf{W}^{\infty+1})=0$, so $\mathbf{W}^{\infty}=\mathbf{W}^{\infty+1}$. We have for the $00$-elements $\mathbf{W}^{\infty}_{00}=\mathbf{W}^{\infty+1}_{00}$, so $x^{\infty}_{0}=x^{\infty+1}_{0}\geq 0$. As we have assumed $y_{0}>0$, we must have $x^{\infty}_{0}=x^{\infty+1}_{0}>0$, otherwise $\mathcal{I}(y_{0}||(x^{\infty}*x^{\infty})_{0})=\mathcal{I}(y_{0}||(x^{\infty}_{0})^{2})$ would be infinite. Because $\mathbf{W}^{\infty}_{0j}=x^{\infty}_{0}x^{\infty}_{j}$ and $\mathbf{W}^{\infty+1}_{0j}=x^{\infty+1}_{0}x^{\infty+1}_{j}$ are equal we can now conclude $x^{\infty}_{j}=x^{\infty+1}_{j}$ for $j>0$.

The second assertion follows from the just given proof. $\Box$

Numerical examples, see Section 5, indicate that depending on the initial values of the algorithm convergence of the iterates takes place to a point that is a local minimizer of the objective function. Comparing the full support Problem 2 to the half support problem discussed in [10], one observes that in the present case there are about twice as many parameters involved when the number of observations is kept the same. Therefore it is not surprising that multiple local minima exist and that convergence of the algorithm to a unique limit, that should also be a global minimizer is at least questionable. The next theorem is therefore the best possible result on the convergence properties of the algorithm.

• Theorem 4.4

Limit points of the sequence $(x^{t})$ are Kuhn-Tucker points of the minimization Problem 2.

• Proof

The main idea here is to reason as in the proof of Proposition 4.8 in [10]. First we invoke the second assertion of Lemma 4 to deduce by taking limits in (20) that

[TABLE]

Hence one has

[TABLE]

If $x^{\infty}_{j}>0$, then necessarily $\nabla_{j}\mathcal{I}(x^{\infty})=0$. If $x^{\infty}_{j}=0$, the analysis is more delicate to deduce that $\nabla_{j}\mathcal{I}(x^{\infty})\geq 0$ and analogous to the proof of Proposition 4.8 in [10]. Further details are omitted. $\Box$

5 Numerical examples

In this section we illustrate the behaviour of Algorithm 4.1 reviewing the results of numerical experiments. All figures are collected at the end of the paper. In the first experiment the data sequence $y$, with $m=20$, was obtained as the autoconvolution $y=x*x$ of a randomly generated signal $x$, with components $x_{j}$ uniformly distributed in the (real) interval $[1,10]$. We investigated whether in this ”exact” case the algorithm is capable of retrieving the true sequence $x$ generating the data $y=x*x$. Figures 1, 2, and 3 share the same data sequence $y$, obtained as specified above. The three figures show the results of three runs of Algorithm 4.1, with $T=2000$, started from three different (randomly generated) initial conditions $x^{0}$. In each figure the top graph shows, in distinct colors, the trajectories of the iterates of the components $x^{t}_{j}$, plotted against the iteration number $t\in[1,T]$. The diamonds at the right end of the graph are the true $x_{j}$ values that generated the data $y$. The second graph shows the superimposed plots of the data generating signal $x$, and of the reconstructed signal $x^{T}$, at the last iteration, both plotted against their component number (in the figure labelled $j=1,2,\dots,m+1$ instead of $j=0,1,\dots,m$). The third graph shows the decreasing sequence $\mathcal{I}(y||x^{t}*x^{t})$. The fourth and last graph shows the superimposed plots of the data vector $y$ and of the reconstructed convolution $x^{T}*x^{T}$, at the last iteration, both against the component number (labelled $j=1,\dots,2m+1$). Figure 1 shows a typical behavior, with the Algorithm producing limit points $x^{\infty}$ that coincide with the generating vector $x$. Another common behavior is shown in Figure 2 which shows convergence (the iterates in the top graph stabilize for $t\rightarrow\infty$) but the limit point $x^{\infty}$ does not coincide with the generating $x$. An atypical but interesting behavior is shown in Figure 3 where, after approximately 1600 iterations, the iterates $x^{t}$ suddenly switch basin of attraction as indicated by the sudden dip in the divergence $\mathcal{I}(y||x^{t}*x^{t})$.

For the second experiment the data $y$ are not produced by autoconvolution but randomly generated. The elements of $y$ are defined as $y_{k}=(k+1)u_{k}$, for $k=0,\ldots,n$, where the $u_{k}$ are generated as IID random numbers from a uniform distribution on $[1,2K^{2}]$, whose expectation is approximately $K^{2}$. Here $K$ is an integer, chosen with a convenient value. This procedure is motivated by comparison with the first experiment, the exact case in which $y=x*x$ with the $x_{k}$’s randomly generated as uniform on $[1,2K]$. In that case the expectation of $x_{i}x_{j}$ (for $i\neq j$) is roughly equal to $K^{2}$ and the autoconvolution $(x*x)_{k}$ will then have expectation approximately equal to $(k+1)K^{2}$. The choice to generate the $y_{k}$ as described above has been made to generate components $y_{k}$ which have an upward trend, as is the case when $y_{k}=(x*x)_{k}$. It would not make a lot of sense to model observations $y_{k}$ as the output of an autoconvolution system, if e.g. the components $y_{k}$ show a decreasing behavior. Figures 4 and 5 share the same data sequence $y$, generated as specified above with $m=12$ and $K=5$. The two figures show the results of two runs of Algorithm 4.1, with $T=2000$, started from two different (randomly generated) initial conditions. In each figure the top graph shows, in distinct colors, the trajectories of the iterates of the components $x^{t}_{j}$, plotted against the iteration number $t\in[1,T]$. The second graph shows the decreasing sequence $\mathcal{I}(y||x^{t}*x^{t})$. The third and last graph shows the superimposed plots of the data vector $y$ and of the fitted convolution $x^{T}*x^{T}$, at the last iteration, both against the component number (labelled $j=1,\dots,2m+1$). Figure 4 shows a typical behavior of the algorithm. The iterates $x_{k}^{t}$ converge quickly to a limit point $x^{\infty}$, as shown by their stabilizing trajectories in the top graph and quick decrease in the divergence $\mathcal{I}(y||x^{t}*x^{t})$. A different run on the same random data $y$, starting from different initial conditions $x^{0}$, produces Figure 5. Also in this case the iterates $x_{k}^{t}$ (top graph) stabilize for $t\rightarrow\infty$, but at a much slower rate, as indicated by the divergence which only after about 1000 iterations dips to reach a limit value much higher than in the first run. The numerical behaviors, both in the exact and in the random data cases confirm the strong dependence of the limit points of the algorithm from the initial conditions, as is expected from this class of algorithms.

Appendix A The solution to Equation (12)

Recall from Lemma 3 that for Problem 3 one has to solve, for a special choice of $r_{j}$ as in (11), the equations given by (12). We repeat these equations here, but in this section we allow for arbitrary positive $r_{j}$. So, we consider

[TABLE]

and want to find the solutions denoted $x_{j}^{\star}$ with positive $r_{j}$ given. As the analytic expressions for the solution turn out to be different for even $n$ and odd $n$, we treat them separately in the next two sections. First some preliminaries that are generally valid.

• Lemma A.1

Consider the system of equations (21). Assume

[TABLE]

and $S=\sqrt{S^{2}}$. Let $x_{j}^{\star}$ be the solutions. Then it holds that $H_{n}:=\sum_{j\leq n/2}x^{\star}_{j}=S$.

• Proof

For convenience we write $k=\lfloor\frac{n}{2}\rfloor$, so $S^{2}=\sum_{j\leq k}r_{j}-\sum_{j\geq k+1}r_{j}$. Add now the first $k+1$ equations of (21) and subtract the last $n-k$. The resulting right hand side is then exactly $S^{2}$ and the left hand side becomes

[TABLE]

Consider the double sum in $H_{n}^{2}$, which is equal to

[TABLE]

which yields for the solutions $x_{j}^{\star}$ the identity H_{n}^{2}=\big{(}\sum_{i=0}^{k}x^{\star}_{i}\big{)}^{2}=S^{2}. $\Box$

Below we need for $j\in\{0,\ldots,n\}$ the quantities

[TABLE]

A.1 The case of $n$ even

Assume $n$ is even, $n=2k$. Consider the equation in (21) for $j=k$,

[TABLE]

Lemma A gives $x_{k}^{\star}=r_{k}/S$. The complete solution to that equation is the content of Proposition A.1. Before formulating the results of that proposition in full detail, we outline the underlying structure. Starting from the ‘middle’ term $x^{\star}_{k}=\frac{r_{k}}{S}$ one sets up two recursions, for $x^{\star}_{k+l}$ and for $x^{\star}_{k-l}$, with $0\leq l\leq k$ that are solved in an alternating way. The explicit solution is presented in the proposition.

• Proposition A.2

Assume $n=2k$ and condition (22). Then the system of equations (21) has the unique solutions $x^{\star}_{j}$ given by the following expressions. For $j=k$ one has $x^{\star}_{k}=\frac{r_{k}}{S}$. For $m>k$ one has

[TABLE]

and for $m\leq k$

[TABLE]

• Proof

Recall $n=2k$. It follows from Lemma A that for the solutions one necessarily has $\sum_{j\leq k}x^{\star}_{j}=S$. The following procedure shows uniqueness. Knowing $x_{k}^{\star}$, we look at the equation for $x_{k+1}$, rewritten as $x_{k+1}(S-x^{\star}_{k})=r_{k+1}$, which gives $x^{\star}_{k+1}$. Then we consider $x_{k-1}(S+x^{\star}_{k+1})=r_{k-1}$, from which we can now find $x^{\star}_{k-1}$. Iterating this procedure gives all solutions $x_{j}^{\star}$.

Having shown the uniqueness, we now verify that the $x_{m}^{\star}$ given by Equations (25) and (26) satisfy the equations. First we claim for $p\leq k$

[TABLE]

One verifies this equality easily for $p=0$ by application of Lemma A, and we proceed by induction. Suppose this formula holds true for some $p-1$, where $p\geq 1$. Then we compute $S-\sum_{l=0}^{p}x^{\star}_{k-l}=S-\sum_{l=0}^{p-1}x^{\star}_{k-l}-x^{\star}_{k-p}$, which becomes

[TABLE]

and this equals

[TABLE]

which in turn is equal to

[TABLE]

The last expression equals

[TABLE]

since $r_{k-p}=B_{k-p}-B_{k-p-1}$. But this is just (27) and proves the claim.

Another computation gives

[TABLE]

We can now compute the product $(S-\sum_{l=0}^{p-1}x^{\star}_{k-l})x^{\star}_{k+p}$. From (27) and (28) we see that the product becomes $r_{k+p}$, as desired.

In the next part of the proof we need $S+\sum_{l=1}^{p}x^{\star}_{k+l}$, $p\geq 0$, which we claim to be equal to

[TABLE]

As above, this can be proved by induction. It is trivial to check this for $p=1$. Assuming the identity to hold for $p-1$ ($p\geq 1$), we compute

[TABLE]

which proves the claim. The next fact to verify is the value of the product x^{\star}_{k-p}\big{(}S+\sum_{l=1}^{p}x^{\star}_{k+l}\big{)}. For that we need a result in the opposite direction (compared to (28)). For $l\geq 1$ one has

[TABLE]

A combination of (29) and (30) shows that this product equals $r_{k-p}$. $\Box$

A.2 The case of $n$ odd

Here we consider the case where $n$ is odd, $n=2k+1$. Recall the result of Lemma A, $\sum_{j=0}^{k}x_{j}^{\star}=S$, with now $S^{2}=\sum_{j=0}^{k}r_{j}-\sum_{j=k+1}^{n}r_{j}=B_{k}-E_{n-k}$. The equation in (12) for $j=k+1$ is $x_{k+1}\sum_{i=0}^{k}x_{i}=r_{k+1}$. Hence $x_{k+1}^{\star}=\frac{r_{k+1}}{S}$. The procedure to get the other solutions $x_{j}^{\star}$ is as in Section A.1, again a double system of recursions is set up. For instance one gets $x_{k}^{\star}=\frac{r_{k}}{S+x_{k+1}^{\star}}=r_{k}\frac{S}{S^{2}+r_{k+1}}=r_{k}\frac{S}{B_{k}-E_{n-k+1}}$, and $x_{k+2}^{\star}=\frac{r_{k+2}}{S-x_{k}^{\star}}=\frac{r_{k+2}}{S}\frac{B_{k}-E_{n-k+1}}{B_{k-1}-E_{n-k+1}}$. The general case is detailed in the following proposition

• Proposition A.3

Let $n=2k+1$. Then $x_{k+1}^{\star}=\frac{r_{k+1}}{S}$. For $m\geq k+1$ one has

[TABLE]

and for $m\leq k$ one has (with $E_{n+1}=0$)

[TABLE]

• Proof

Let the $x_{m}^{\star}$ be given by Equations (31) and (32). Instead of Equation (27) we now have for $p\leq k$

[TABLE]

For $p=0$ this is verified by application of Lemma A. For other values of $p$ one proves this by induction.

Instead of Equation (29) we now have for $p\leq k+1$

[TABLE]

which can again be proved by induction with respect to $p$.

One then verifies that the postulated expressions for the $x_{m}^{\star}$ satisfy the equations by arguments similar to those in the proof of Proposition A.1. $\Box$