# Signal estimation from correlation vectors with structured low-rank approximation

Antonio Fazzi, Julien Flamant, Konstantin Usevich

► To cite this version:

Antonio Fazzi, Julien Flamant, Konstantin Usevich. Signal estimation from correlation vectors with structured low-rank approximation. EUSIPCO 2025 - 33rd European Signal Processing Conference, Sep 2025, Isola delle Femmine – Palermo, Italy. hal-05122986

HAL Id: hal-05122986

https://hal.science/hal-05122986v1

Submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Signal estimation from correlation vectors with structured low-rank approximation

Antonio Fazzi

Université de Lorraine, CNRS, CRAN
F-54000 Nancy, France
antonio.fazzi@univ-lorraine.fr

Julien Flamant

Université de Lorraine, CNRS, CRAN
F-54000 Nancy, France
julien.flamant@cnrs.fr

Konstantin Usevich

Université de Lorraine, CNRS, CRAN
F-54000 Nancy, France
konstantin.usevich@cnrs.fr

**Abstract**—The paper focuses on the reconstruction of a multivariate signal from its noisy correlation vectors, a problem appearing in phase retrieval and blind channel identification. In the noiseless case, the signals can be retrieved from greatest common divisors of the polynomials associated to observed vectors. For the noisy case, by exploiting the properties of Sylvester matrices, we propose a new reconstruction approach that exploits all the available data and preserves matrix structure in the associated low-rank approximation problem. By doing so, we achieve an improvement of the reconstruction performance with respect to existing methods.

**Index Terms**—Signals estimation, Structured Low-Rank Approximation, Approximate GCD, Sylvester matrices.

## I. INTRODUCTION

This paper deals with the estimation of signals from their auto-correlation and cross-correlation functions. This problem arises notably in blind channel identification [1], [2] and phase retrieval [3], [4]. The problem can be mathematically formulated as follows. Consider $R$ discrete signals $\mathbf{x}_i \in \mathbb{C}^{N_i}$, $i = 1, \dots, R$ and let $\gamma_{ij} = \mathbf{x}_i \star \mathbf{x}_j \in \mathbb{C}^{N_i+N_j-1}$ denote the vector encoding the auto- or cross-correlation between the $(i, j)$-th signal pair. The estimation problem thus consist in

$$ \text{find } \{\mathbf{x}_i\}_{i=1}^R \text{ given } \{\gamma_{ij} = \mathbf{x}_i \star \mathbf{x}_j\}_{i,j=1}^R. \quad (1) $$

For $R = 1$, the problem amounts at the classical spectral factorization problem [5] of univariate polynomials. The case $R = 2$ has been studied in several publications [6]–[8] while the general case $R \ge 2$ can be interpreted as a rank-one factorization problem of a matrix polynomial built from correlation functions [9]. The study of theoretical properties of the problem (1) relies on a classical polynomial reformulation thanks to z-transform-like expansions, which is particularly useful for studying uniqueness of solutions. It is well known (see e.g., [10]) that the case $R = 1$ does not admit a unique solution; uniqueness can be enforced through additional constraints (such as minimum phase) – see [4] for a review of such strategies. On the other hand, the case $R \ge 2$ enjoys strong uniqueness properties for generic signals [9]. Moreover, it can be shown that the problem admits a unique algebraic solution as the greatest common divisor (GCD) of correlation polynomials: this feature has been exploited recently [7] to

This work was supported in part by the ANR through grants RICOCHET ANR-21-CE48-0013 and ATEMPORAL ANR-23-CE48-0007.

devise an efficient algebraic reconstruction strategy for the polarimetric phase retrieval problem.

Existing GCD-based reconstruction strategies [6], [7] for (1) suffer from poor performance in noisy scenarios. This work addresses this limitation by studying novel structured low-rank approximation (SLRA) approaches for the problem (1). We propose several Sylvester-like matrix structures to take into account the specificities of the correlation measurements in (1). In addition, the proposed approach leverages recent advances in approximate GCD computations [11], [12], by efficiently reformulating the SLRA optimization problem as an ODE system. Numerical experiments demonstrate that taking into account the structure of (1) together with dedicated SLRA algorithms yields significant improved reconstruction performance in noisy scenarios, while maintaining reasonable computational burden.

**Paper organization:** Section II reviews the classic polynomial representation of (1) and the properties of Sylvester matrices. Section III introduces the SLRA problem, relevant Sylvester-like matrix structures, and the ODE-based gradient system algorithm. Numerical experiments are presented in Section IV and Section V gathers concluding remarks.

## II. PRELIMINARIES

To simplify the presentation, we assume from now on $R = 2$ and $N_1 = N_2 = N$, which is sufficient to demonstrate the rationale of the approach. The case $N_1 \ne N_2$ can be recovered by careful adaption of dimensions in subsequent definitions. In addition, the proposed tools can be extended to $R > 2$ signals by using techniques presented in [13].

### A. Polynomial representation

Consider two signals $\mathbf{x}_1, \mathbf{x}_2 \in \mathbb{C}^N$ and define their polynomial representation $X_1, X_2 \in \mathbb{C}_{\le N-1}[z]$ such that $X_i(z) = \sum_{k=0}^{N-1} x_i[k]z^k$ for $i = 1, 2$. In addition, let us define the polynomials with reversed and conjugated coefficients $\tilde{X}_i(z) = \sum_{k=0}^{N-1} x_i^*[N - 1 - k]z^k$ where $\cdot^*$ denotes the complex conjugate operation. It can be observed that $\tilde{X}_i(z) = z^{N-1}X_i^*((z^*)^{-1})$. Note that these polynomial representations are much alike to the standard z-transform of discrete signals, except for the use of positive powers of $z$ – a standard convention in the Fourier phase retrieval literature [4], [6], [10]. Similarly, for any $(i, j)$ pair, one can define the correlation
polynomial $\Gamma_{ij}(z) = \sum_{k=0}^{2N-2} \gamma_{ij}[k-N+1]z^k$. Standard calculations (see e.g., [4]) show that $\Gamma_{ij}(z) = X_i(z)\tilde{X}_j(z)$. Therefore, the problem (1) is equivalent to the following polynomial factorization problem

$$ \text{find } X_1(z), X_2(z) \text{ given } \{\Gamma_{ij}(z) = X_i(z)\tilde{X}_j(z)\}_{i,j=1}^2. \quad (2) $$

The formulation (2) allows for a complete characterization of solutions to the original problem (1) [9]. In particular, when polynomials $X_1(z)$, $X_2(z)$ are co-prime (equivalently, when all $\{\Gamma_{ij}(z)\}$ are co-prime) then the problem admits a unique solution, which can be expressed using the greatest common divisor of correlation polynomials, as stated below.

**Proposition 1** (GCD-based recovery [7]). Let $X_1, X_2 \in \mathbb{C}_{\le N-1}[z]$ such that $\gcd(X_1, X_2) = 1$. Then $X_1(z)$ and $X_2(z)$ can be uniquely recovered from $\{\Gamma_{ij}(z)\}_{i,j=1}^2$ as

$$ X_1(z) = c_1 \gcd(\Gamma_{11}, \Gamma_{12}) \quad \text{and} \quad X_2 = c_2 \gcd(\Gamma_{21}, \Gamma_{22}) \quad (3) $$

where $c_1, c_2 \in \mathbb{C}$ can be determined explicitly (up to global phase, i.e., a common scalar complex factor of modulus 1) from correlation polynomials.

Proposition 1 has two key consequences. First, since for generic polynomials $\gcd(X_1, X_2) = 1$ almost surely, this means that the recovery problem (2) (and thus (1) as well) admits a unique solution of the form (3). Second, it suggests an algebraic reconstruction strategy based on numerical methods for computing the (approximate) greatest common divisor (GCD) of polynomials [14]. These rely on Sylvester-like matrices and their kernel properties, which are reviewed next.

## B. Sylvester matrices and their (right) kernel

Low-rank approximation of Sylvester (or Sylvester-like matrices) matrices [15]–[17] is a common approach to find approximate greatest common divisors in numerical linear algebra. First, for a polynomial $P(z) \in \mathbb{C}_{\le L}[z]$ with coefficients $p \in \mathbb{C}^{L+1}$, we define the multiplication matrix as

$$ \mathbf{M}_K(P) = \begin{bmatrix} p_0 & & & \\ & \ddots & & \\ & & p_L & \\ & & & \\ & & & p_L \end{bmatrix} \in \mathbb{C}^{(L+K+1) \times (K+1)}, $$

which is a matrix representation of multiplication of $P(z)$ by any polynomial of degree $\le K$. For two polynomials $P(z), Q(z) \in \mathbb{C}_{\le L}[z]$ and a number $D$, let

$$ \mathbf{S}_D(P, Q) = [\mathbf{M}_{L-D}(Q) - \mathbf{M}_{L-D}(P)] \quad (4) $$

denote the modified<fnref n="1" /> Sylvester (subresultant) matrix. Then the following proposition is known.

**Proposition 2** (Right kernel property [18], simplified). Consider two polynomials $P, Q \in \mathbb{C}_{\le L}[z]$ such that $\gcd(P, Q) =$

## Algorithm 1 Right kernel Sylvester

**Input:** Sylvester subresultant matrix $\mathbf{S}$ with $2N$ columns, an estimate $\hat{C}$ of $\|[\hat{\mathbf{x}}_1, \hat{\mathbf{x}}_2]\|_2$

1. Take $\mathbf{v} = \mathbf{v}_{2N} \in \mathbb{C}^{2N}$ to be the $2N$-th right singular vector of $\mathbf{S}$ (corresponding to the $2N$-th singular value)

2. Partition $\mathbf{v}$ as $\mathbf{v} = [\mathbf{v}_1, \mathbf{v}_2]$.

3. Set $\hat{\mathbf{x}}_1 = \hat{C}\mathbf{v}_1$ and $\hat{\mathbf{x}}_2 = \hat{C}\mathbf{v}_2$

**Output:** Estimates $\hat{\mathbf{x}}_1$ and $\hat{\mathbf{x}}_2$

$H \in \mathbb{C}_{\le D}[z]$ and let $P(z) = H(z)U(z)$, $Q(z) = H(z)V(z)$, where $\bar{U}, V \in \mathbb{C}_{\le L-D}[z]$ are quotient polynomials. Then the right kernel of $\mathbf{S}_D(P, Q)$ has dimension one and is spanned by the concatenation of the coefficients of quotient polynomials:

$$ \mathbf{S}_D(P, Q) \begin{bmatrix} \mathbf{u} \\ \mathbf{v} \end{bmatrix} = 0. $$

Propositions 1 and 2 lead to a practical reconstruction strategy for the problem (2), called right kernel Sylvester approach, summarized in Algorithm 1. It enables the recovery of polynomials $X_1(z)$ and $X_2(z)$ up to a global phase factor from an input Sylvester matrix with $2N$ columns constructed from correlation polynomials $\{\Gamma_{ij}(z)\}$. A simple choice, used e.g., in [7], is to take $P(z) = \Gamma_{11}(z)$ and $Q(z) = \Gamma_{21}(z)$ in (4) leading to

$$ \mathbf{S}_{\text{partial}} := \mathbf{S}_{N-1}(\Gamma_{11}, \Gamma_{21}) \in \mathbb{C}^{(3N-2) \times 2N}, \quad (5) $$

where the subscript refers to the partial use of correlation polynomials in (2). Note that a typical estimate of the norm in Algorithm 1 is $\hat{C} = \sqrt{\gamma_{11}^2[0] + \gamma_{22}^2[0]}$. While being very computationally efficient, i.e., it only requires the computation of one SVD, this approach exhibits poor reconstruction performance in noisy scenarios [7], as the Sylvester matrix structure is not preserved by the low-rank SVD approximation. The next section precisely address this issue through several structured low-rank approximation algorithms.

## III. PROPOSED APPROACH

In practical settings, one has only access to noisy measurements of correlations polynomials $\{\Gamma_{ij}(z)\}$ in (2). In this paper, we assume that the $(i, j)$-th correlation polynomial is given by

$$ \Gamma_{ij}(z) = \Gamma_{ij}^0(z) + \alpha N_{ij}(z), \quad (6) $$

where $\Gamma_{ij}^0(z) = X_i(z)\tilde{X}_j(z)$ is the noiseless correlation polynomial and where $N_{ij} \in \mathbb{C}_{\le 2N-2}[z]$ is a noise polynomial whose coefficients are i.i.d. complex circular Gaussian variables. Finally, the parameter $\alpha \ge 0$ controls the relative noise contribution in (6) and can be interpreted as the standard deviation of the noise polynomial coefficients.

### A. Improved Sylvester-like matrix structure

As explained above, any method exploiting the right-kernel properties of Sylvester matrices (such as Algorithm 1) requires matrix structures such that Proposition 2 (or variations) holds. A standard example is $\mathbf{S}_{\text{partial}}$ (see Eq. (5)); another straightforward example would be $\mathbf{S}'_{\text{partial}} = \mathbf{S}_{N-1}(\Gamma_{12}, \Gamma_{22})$, which

<footnote n="1">¹The definition used here differs from the standard one [14] by the ordering and changing sign of the blocks. While it does not affect rank properties, it simplifies the reconstruction procedure in Algorithm 1.</footnote>
uses the complementary pair of correlation polynomials and encodes the same quotient polynomials as $S_{partial}$.

In this paper, we propose to fully exploit all the polynomials in $\{\Gamma_{ij}(z)\}$ by considering the block Sylvester matrix

$$S_{\text{full}} = \begin{bmatrix} M_{N-1}(\Gamma_{21}) & -M_{N-1}(\Gamma_{11}) \\ M_{N-1}(\Gamma_{22}) & -M_{N-1}(\Gamma_{12}) \end{bmatrix} \in \mathbb{C}^{(6N-4) \times 2N}. \quad (7)$$

The motivation of this choice is that, despite the two pairs of polynomials having different GCDs (see Proposition 1), the two block rows of $S_{full}$ share the same one dimensional right kernel (that is, the polynomials have the same cofactors $X_1, X_2$), so does the overall structured matrix in (7). By encoding all available information in a single matrix $S_{full}$, one can expect improved robustness to noise compared to partial structures as in (5).

## B. Structured low-rank approximation

The issue with using just the plain SVD in Algorithm 1 is that it does not preserve the (Sylvester) matrix structure in the low-rank approximation. This is known to be suboptimal in signal reconstruction problems and approximate polynomials GCD computations [11], [13]. To circumvent this issue, we propose to use the powerful structured low-rank approximation (SLRA) [19] approach, which we detail below.

Let $S \subset \mathbb{C}^{m \times n}$ be a set (affine subspace) of matrices of some fixed structure (e.g., Sylvester, Hankel, etc). We will refer to $S$ as the matrix structure. Then the SLRA problem (with rank reduction by 1 [13], [19]) is to find a closest rank-deficient structured matrix to a given one.

**Problem (SLRA with rank reduction by 1).** Given a structured matrix $S \in S$, compute

$$\min_{\hat{S} \in S} \|S - \hat{S}\|_F, \quad \text{s.t.} \quad \sigma_{\min}(\hat{S}) = 0, \quad (8)$$

where $\sigma_{\min}(\cdot)$ denotes the smallest singular value.

Rank reduction by one suits well our needs since both the matrices $S_{partial}$ and $S_{full}$ are tall, and have one dimensional right kernel under assumptions of Proposition 2.

As we discussed before, there are several matrix structures that encode the polynomial factorization property. We consider the following three matrix structures:

1) $S_{partial}$: the ordinary Sylvester structure (5);

2) $S_{full}$: the block Sylvester structure in (7), that exploits all elements of $\{\Gamma_{ij}(z)\}$;

3) $S_{full}^{corr}$: the block Sylvester structure from (7), with a symmetry imposed on the elements of the matrix, which reflects the symmetry possessed by the coefficients of the correlation polynomials $\{\Gamma_{ij}(z)\}$ (see (2)):

$$S_{\text{full}}^{\text{corr}} = \{S_{\text{full}} : \tilde{\Gamma}_{ij}(z) = \Gamma_{ji}(z), \forall i, j\}.$$

The matrix structure $S_{full}^{corr}$ is particularly interesting as it is the smallest matrix structure of the three that contains autocorrelation matrix polynomials. Note that $S_{full}^{corr}$ is an $\mathbb{R}$-linear subspace, unlike $S_{partial}$ and $S_{full}$ which are $\mathbb{C}$-linear subspaces; this fact has some repercussions for the numerical solution of SLRA (see the next subsection).

## C. An ODE-based method for block Sylvester SLRA

Many methods exist for SLRA with affine matrix structure (see [11, Section 1.2] for a summary). In this paper, we employ the ODE-based methodology proposed recently [11], which consists in reformulating the SLRA problem as matrix ODE. The advantages of such a methodology are twofold: (a) it already proved to be efficient for Sylvester-like structures [20], [21], complex polynomials [11, Section 5.2] and block structures [22]; and (b) it is easy to implement since we only need to change the projection $P_S(\cdot)$ on the matrix structure. We briefly recall here the main ingredients of the ODE-based algorithm for the SLRA problem (8) and we refer the interested reader to [11], [12] for proofs or other details.

Let $S$ be a full rank structured matrix. This is the case, e.g., when building $S_{partial}$ or $S_{full}$ from noisy correlation polynomials (6). Considering the SLRA problem (8), the basic idea of the ODE-based method [11] is to search for a solution $\hat{S} = S + \epsilon E$ which is a linear perturbation of the input matrix $S$. The perturbation $\epsilon E$ is such that $E$ is a normalized (with respect to the Frobenius norm) structured matrix that minimizes the smallest singular value $\sigma_{\min}(\hat{S})$ over the ball of perturbation matrices $\epsilon E$ of norm $\epsilon$.

Starting from an initial guess $\hat{S}_0 = S + \epsilon_0 E_0$, this splitting allows for updating $\epsilon$ and $E$ iteratively on two levels:

1) at an inner level, the norm $\epsilon$ is fixed, and we look for the structured matrix $E(\epsilon)$ that minimizes $\sigma_{\min}(S + \epsilon E(\epsilon))$;

2) at an outer level, if the computed $\hat{S}$ does not satisfy (at least approximately) $\sigma_{\min}(\hat{S}) \approx 0$, we update $\epsilon$.

At the inner level, we look for a zero of the derivative of $\sigma_{\min}$. By using results about derivatives of singular values, the computation of $E = E(\epsilon)$ can be reformulated [11] as finding stationary points of the ODE:

$$\dot{\mathbf{E}} = -\mathcal{P}_S(\mathbf{u}\mathbf{v}^H) + \langle \mathbf{E}, \mathcal{P}_S(\mathbf{u}\mathbf{v}^H) \rangle \mathbf{E}. \quad (9)$$

where $\mathbf{u}$ and $\mathbf{v}$ are the singular vectors associated with $\sigma_{\min}, \cdot^H$ is the standard conjugate transpose operator and $\langle \cdot, \cdot \rangle$ denotes the usual Frobenius inner product for matrices. At the outer level, we look for $\bar{\epsilon}$ that annihilates $\sigma_{\min}(S + \epsilon E)$. This is computed by the Newton method coupled with a bisection step, to look for possible better solutions once we find an admissible one. We summarize the procedure in Algorithm 2, where computational steps are explained in [11], [12].

The orthogonal projections on $S_{partial}$ and $S_{full}$ are standard and calculated by averaging along diagonals and putting elements to zero. The structure $S_{full}^{corr}$ is particular due to its $\mathbb{R}$-linear structure. However, the projection on $S_{full}^{corr}$ can be still simply computed as follows:

1) compute the orthogonal projection of a matrix $X$ on $S_{full}$, and let $\hat{\gamma}_{ij}, i, j = 1, 2$ be the four estimated polynomials coefficients for the corresponding blocks;

2) symmetrize polynomials' coefficients as:

$$\hat{\gamma}_{ij} \leftarrow \frac{1}{2}(\hat{\gamma}_{ij} + \text{flip}(\hat{\gamma}_{ji}^*))$$

where the flip(·) operator reverses the order of polynomial coefficients.
**Algorithm 2 ODE-based algorithm for the SLRA problem (8)**

**Input:** Sylvester matrix $\mathbf{S}$ with $2N$ columns and corresponding structure $\mathcal{S}$, bounds for perturbation norm $[\epsilon_{\min}, \epsilon_{\max}]$, maximum iteration number $k_{\max}$, tolerance $T$. Estimate the starting values $\epsilon_0, \mathbf{E}_0$;

Set $\mathbf{u}, \sigma, \mathbf{v}$ as the smallest singular triplet of $\hat{\mathbf{S}} = \mathbf{S} + \epsilon_0 \mathbf{E}_0$; for $k = 1 : k_{\max}$ do

1. Compute $\epsilon_k$ by a Newton-bisection method:

if $\sigma > T$ then

$$ \epsilon_k = \epsilon_{k-1} + \frac{\sigma}{\|\mathcal{P}_S(\mathbf{u}\mathbf{v}^*)\|_F}, $$

else

$$ \text{Set } \epsilon_{\max} = \epsilon_k \text{ and } \epsilon_k = 0.5\epsilon_{\min} + 0.5\epsilon_{\max}. $$

end if

2. Compute the perturbation $\mathbf{E}$ corresponding to $\sigma_{\min}(\mathbf{S} + \epsilon_k \mathbf{E})$ by integrating the ODE (9)

3. Compute the smallest singular triplet $(\mathbf{u}, \sigma, \mathbf{v})$ from $\hat{\mathbf{S}} = \mathbf{S} + \epsilon_k \mathbf{E}$

if $\sigma < T$ then

stop

end if

end for

Apply Algorithm 1 to $\hat{\mathbf{S}}$ with $\hat{\mathbf{C}} = \sqrt{\hat{\gamma}_{11}^2[0] + \hat{\gamma}_{22}^2[0]}$

**Output:** Estimates $\hat{\mathbf{x}}_1$ and $\hat{\mathbf{x}}_2$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Average error comparison (in logarithmic scale) among different computational method for the estimation of a signal from correlation polynomials</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Average computational times comparison (in logarithmic scale) for the problem of estimating signals of increasing length.</figcaption></figure>

## IV. NUMERICAL EXPERIMENTS

This section deals with the evaluation of the numerical performances of the proposed Sylvester-like matrix structures and associated algorithms to solve problem (1). To illustrate the respective role of each matrix structure and each algorithm, we consider multiple combinations. In total, we may consider three different matrix structures $S_{\text{partial}}, S_{\text{full}}$ and $S_{\text{full}}^{\text{corr}}$ of increasing complexity. Algorithm 1 can be used with matrices in $S_{\text{partial}}$ and $S_{\text{full}}$; corresponding strategies will be labeled as rker($S_{\text{partial}}$) and rker($S_{\text{full}}$). Algorithm 2 can itself be used with any of the three matrix structures considered in this paper; corresponding strategies are therefore labeled as Sylv-ODE($\mathcal{S}$), where $\mathcal{S} = S_{\text{partial}}, S_{\text{full}}$ or $S_{\text{full}}^{\text{corr}}$ is the considered structure. For comparison, we also implement the classical Cadzow's algorithm [23] using the full correlation structure $S_{\text{full}}^{\text{corr}}$, where the final step of the algorithm calls the right kernel method (Algorithm 1) on the estimated Sylvester-like matrix.

All numerical experiments are performed on a 2024 Apple M4 MacBookPro with 24GB RAM using Matlab 2024b. Due to the inherent global phase ambiguity in (1), reconstruction error is defined as follows. Letting $\hat{\mathbf{x}} = [\hat{\mathbf{x}}_1, \hat{\mathbf{x}}_2] \in \mathbb{C}^{2N}$ be an estimation of the (vectorized) ground truth signal $\mathbf{x} = [\mathbf{x}_1, \mathbf{x}_2] \in \mathbb{C}^{2N}$, the relative estimation error $r(\hat{\mathbf{x}})$ is defined as

$$ r(\hat{\mathbf{x}}) = \min_{\Phi \in [0, 2\pi)} \frac{\|\hat{\mathbf{x}}e^{i\Phi} - \mathbf{x}\|_2^2}{\|\mathbf{x}\|_2^2} \quad (10) $$

which is the standard metric in, e.g., phase retrieval [3], [4].

Figure 1 depicts the average reconstruction error as a function of the noise standard deviation $\alpha$ for the observation model (6). The same ground-truth signals $\mathbf{x}_1, \mathbf{x}_2 \in \mathbb{C}^{32}$

(generated as i.i.d. complex circular Gaussian vectors) were used in the 50 independent Monte-Carlo runs for each value of $\alpha$ ranging logarithmically from $10^{-6}$ to $10^{-1}$. First observe that the use of the structure $S_{\text{full}}$ or $S_{\text{full}}^{\text{corr}}$ systematically leads to an improvement in the estimation error with respect to $S_{\text{partial}}$: this happens both by comparing rker($S_{\text{partial}}$) and rker($S_{\text{full}}$) as well as by comparing Sylv-ODE($S_{\text{partial}}$) with Sylv-ODE($S_{\text{full}}$) and Sylv-ODE($S_{\text{full}}^{\text{corr}}$). The second observation is about the comparison between the full rank and the rank deficient matrices. Indeed, for all noise levels all approximated matrices yield a smaller error, independently on the considered Sylvester structure with 2 or 4 blocks. Thus the SLRA approach for the signal estimation problem is beneficial (at least for the considered algorithms), from the point of view of the error on the computed solution.

Figure 2 displays computational timings in the case $\alpha = 10^{-3}$ for increasing size of signals from $N = 4$ to $N = 64$. Computed times were averaged over 10 independent reconstructions. The complexity of all the proposed methods is
dominated by the computation of the (full) SVD, thus it can be
estimated from the number of SVD factorizations performed.
Methods rker(S<sub>partial</sub>) and rker(S<sub>full</sub>) only need one SVD for
the direct error estimation. Cadzow's algorithm uses one SVD
per iteration, while the ODE-based methods integrate one ODE
(9) at each iteration, and such iterative integration needs a
SVD per iteration (see [11, Algorithm 1]). The interesting
observation (that is only a numerical evidence, up to now)
is the fact that Sylv-ODE(S<sub>full</sub><sup>corr</sup>) runs much faster than Sylv-
ODE(S<sub>full</sub>).

Overall, Cadzow’s iteration with the projection on S<sub>full</sub><sup>corr</sup> looks as a good trade-off between the computational complexity and the reconstruction error. Note that its performance can probably be further improved by exploiting variations of this classical algorithm, such as stochastic optimization techniques [24].

V. CONCLUDING REMARKS

In this paper, we studied the recovery of two signals from
their auto- and cross-correlation vectors, which is an important
problem arising in blind channel identification and phase
retrieval. This problem is notoriously difficult due to its non-
convex nature. The originality of the proposed approach lies in
exploiting a polynomial representation of the original problem,
enabling solutions through approximate GCD computations.
To address limitations of current methods – notably their
poor robustness to noise – we proposed a structured low-
rank approximation (SLRA) approach using Sylvester-like
matrices. By leveraging all available correlation polynomials
and preserving the matrix structure, we demonstrated signif-
icant improvements in reconstruction performance. Numeri-
cal experiments confirmed the effectiveness of the proposed
approach, particularly when using the full block Sylvester
structure with correlation constraints. Future work will focus
on 1) improving the computational performances of ODE-
based algorithms and 2) extending these algebraic GCD-based
approaches to multidimensional signals (e.g., images) or to
other measurement settings such as spectrogram-like measure-
ments [25] in audio signal processing [26] or ptychographic
imaging [27].

REFERENCES

[1] G. Xu, H. Liu, L. Tong, and T. Kailath, "A least-squares approach to blind channel identification," *IEEE Transactions on Signal Processing*, vol. 43, no. 12, pp. 2982-2993, 1995.

[2] L. Tong, G. Xu, B. Hassibi, and T. Kailath, "Blind channel identification based on second-order statistics: A frequency-domain approach," *IEEE Transactions on Information Theory*, vol. 41, no. 1, pp. 329-334, 1995.

[3] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, "Phase Retrieval with Application to Optical Imaging: A contemporary overview," *IEEE Signal Processing Magazine*, vol. 32, no. 3, pp. 87-109, May 2015.

[4] T. Bendory, R. Beinert, and Y. C. Eldar, "Fourier Phase Retrieval: Uniqueness and Algorithms," in *Compressed Sensing and its Applications*. Cham: Springer International Publishing, 2017, pp. 55-91.

[5] A. H. Sayed and T. Kailath, "A survey of spectral factorization methods," *Numerical linear algebra with applications*, vol. 8, no. 6-7, pp. 467-496, 2001.

[6] K. Jaganathan and B. Hassibi, “Reconstruction of signals from their autocorrelation and cross-correlation vectors, with applications to phase retrieval and blind channel estimation,” *IEEE Transactions on Signal Processing*, vol. 67, no. 11, pp. 2937–2946, 2019.

[7] J. Flamant, K. Usevich, M. Clausel, and D. Brie, “Polarimetric Fourier phase retrieval,” *SIAM Journal of Imaging Sciences*, vol. 17, no. 1, pp. 632–671, 2024.

[8] O. Raz, N. Dudovich, and B. Nadler, “Vectorial phase retrieval of 1-D signals,” *IEEE Transactions on Signal Processing*, vol. 61, no. 7, pp. 1632–1643, 2013.

[9] K. Usevich, J. Flamant, M. Clausel, and D. Brie, "On factorization of rank-one auto-correlation matrix polynomials," 2023. [Online]. Available: https://arxiv.org/abs/2308.15106

[10] R. Beinert and G. Plonka, “Ambiguities in One-Dimensional Discrete Phase Retrieval from Fourier Magnitudes,” *Journal of Fourier Analysis and Applications*, vol. 21, no. 6, pp. 1169–1198, Dec. 2015.

[11] A. Fazzi, N. Guglielmi, and I. Markovsky, "An ODE based method for computing the Approximate Greatest Common Divisor of polynomials," *Numerical Algorithms*, vol. 81, pp. 719–740, 2019.

[12] N. Guglielmi and I. Markovsky, “An ODE based method for computing the distance of coprime polynomials,” *SIAM J. Numer. Anal.*, vol. 55, pp. 1456–1482, 2017.

[13] K. Usevich and I. Markovsky, "Variable projection methods for approximate (greatest) common divisor computations," *Theoretical Computer Science*, vol. 681, pp. 176–198, 2017.

[14] K. Nagasaka, "Toward the best algorithm for approximate gcd of univariate polynomials," *Journal of Symbolic Computation*, vol. 105, pp. 4–27, 2021.

[15] J. J. Sylvester, "On a theory of the Syzygetic relations of two rational integral functions, comprising an application to the theory of Sturm's functions, and that of the greatest algebraical common measure," *Philosophical Trans.*, vol. 143, pp. 407–548, 1853.

[16] W. Qiu, Y. Hua, and K. Abed-Meraim, "A subspace method for the computation of the gcd of polynomials," *Automatica*, vol. 33, no. 4, pp. 741–743, 1997.

[17] D. A. Bini and P. Boito, “Structured matrix-based methods for polynomial epsilon-gcd: analysis and comparisons,” in *Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation*, ser. ISSAC ’07. New York, NY, USA: ACM, 2007, p. 9–16.

[18] D. Christou, M. Mitrouli, and D. Triantafyllou, “Structured Matrix Methods Computing the Greatest Common Divisor of Polynomials,” *Special Matrices*, vol. 5, no. 1, pp. 202–224, 2017.

[19] I. Markovsky, "Structured low-rank approximation and its applications," *Automatica*, vol. 44, no. 4, pp. 891–909, 2008.

[20] A. Fazzi, N. Guglielmi, and I. Markovsky, “Generalized algorithms for the approximate matrix polynomial GCD of reducing data uncertainties with application to MIMO system and control,” *J. Comput. Appl. Math.*, vol. 393, 2021.

[21] ——, "Computing common factors of matrix polynomials with applications in system and control theory," in *Proc. of the IEEE Conf. on Decision and Control*, Nice, France, 2019, pp. 7721–7726.

[22] A. Fazzi, N. Guglielmi, I. Markovsky, and K. Usevich, "Common dynamic estimation via structured low-rank approximation with multiple rank constraints," in *IFAC - Papers Online*, vol. 54, 2021, pp. 103–107.

[23] J. Cadzow, "Signal enhancement-a composite property mapping algorithm," *IEEE Transactions on Acoustics, Speech, and Signal Processing*, vol. 36, no. 1, pp. 49–62, 1988.

[24] J. Gillard and A. Zhiglavsky, "Stochastic algorithms for solving structured low-rank matrix approximation problems," *Communications in Nonlinear Science and Numerical Simulation*, vol. 21, no. 1, pp. 70–88, 2015, numerical Computations: Theory and Algorithms (NUMTA 2013), International Conference and Summer School.

[25] K. Jaganathan, Y. C. Eldar, and B. Hassibi, "Stft phase retrieval: Uniqueness guarantees and recovery algorithms," *IEEE Journal of selected topics in signal processing*, vol. 10, no. 4, pp. 770–781, 2016.

[26] P. Magron, R. Badeau, and B. David, "Model-based stft phase recovery for audio source separation," *IEEE/ACM Transactions on Audio, Speech, and Language Processing*, vol. 26, no. 6, pp. 1095–1105, 2018.

[27] O. Melnyk, "Phase retrieval from short-time fourier measurements and applications to ptychography," Ph.D. dissertation, Technische Universität München, 2023.