# Time and covariance smoothing for restoration of bivariate signals

Yusuf Yigit Pilavci, Pierre Palud, Julien Flamant, Pierre-Antoine Thouvenin,

Jérémie Boulanger, Pierre Chainais

► To cite this version:

Yusuf Yigit Pilavci, Pierre Palud, Julien Flamant, Pierre-Antoine Thouvenin, Jérémie Boulanger, et al.. Time and covariance smoothing for restoration of bivariate signals. 2025 IEEE Statistical Signal Processing Workshop, Jun 2025, Edinburgh, United Kingdom. hal-05125331

HAL Id: hal-05125331

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

Submitted on 24 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
# Time and covariance smoothing for restoration of bivariate signals

Yusuf Yiğit PİLAVCI¹, Pierre PALUD², Julien FLAMANT³, Pierre-Antoine THOUVENIN¹, Jérémie BOULANGER¹, and Pierre CHAINAIS¹

¹Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000, Lille, France

²Universite Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris, France

³Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France

**Abstract**—In many applications and physical phenomena, bivariate signals are polarized, i.e. they trace an elliptical trajectory over time when viewed in the 2D planes of their two components. The smooth evolution of this elliptical trajectory, called polarization ellipse, is highly informative to solve ill-posed inverse problems involving bivariate signals where the signal is collected through indirect, noisy or incomplete measurements. This work proposes a novel formulation and an efficient algorithm for reconstructing bivariate signals with polarization regularization. The proposed formulation leverages the compact representation of polarization through the instantaneous covariance matrices. To address the resulting quartic optimization problem, we propose a well-suited parameter splitting strategy which leads to an efficient iterative algorithm (alternating direction method of multipliers (ADMM)) with convex subproblems at each iteration. The performance of the proposed method is illustrated on numerical synthetic data experiments.

**Index Terms**—bivariate signals, denoising, polarization, covariance, ADMM

## I. INTRODUCTION

A myriad of examples involving electromagnetic waves, gravitational waves, seismic waves, underwater acoustics signals or EEG signals falls into the realm of bivariate signal processing (BSP), i.e., the study of signals with two correlated components [7]. In multivariate signal processing, BSP is a special case where many useful representations are accessible. Several seminal works [9], [8], [14] have studied the definition and practical use of quaternions, covariance matrices and Stokes parameters (widely used in optics) to interpret the polarization properties, i.e. the properties related to the elliptical trajectory described by the bivariate signal in the 2D plane. Solving inverse problems involving bivariate signals, such as signal restoration, relies on these representations to formulate polarization priors into regularization [17], [3]. They often lead to non-linear, typically quartic, formulations. As a result, the associated regularizations are often non-convex functions.

The authors acknowledge support from the ANR project "Chaire IA Sherlock" ANR-20-CHIA-0031-01, the programme d'investissements d'avenir ANR-16-IDEX0004 ULNE and Region Hauts-de-France. This research has been supported in whole or in part, by the Agence Nationale de la Recherche (ANR) under project ANR-21-CE48-0013-03. In the interest of open access publication, the author/rights holder applies a CC-BY open access license to any article/manuscript accepted for publication (AAM) resulting from this submission. The corresponding author: firstname.lastname@univ-lille.fr; firstname:yigit, lastname: {pilavci}.

Previously, the Jones polarizer matrix has been used to impose a predefined polarization for reconstructing gravitational waves [4]. This method is very effective when the polarization is known, e.g., circular or linear. However, knowing the polarization drastically limits the range of application and excludes cases of time varying polarization. Closer to this work, [17] proposes two formulations to impose smoothness both in the time and Stokes domains. The numerical results show that the best reconstruction results are obtained by regularizing in both domains. However, the solutions to these non-convex formulations are given by gradient descent based algorithms, which are not necessarily well adapted for the problem at hand. In this work, we propose a new formulation of the signal reconstruction problem from noisy measurements. In doing so, we use the instantaneous covariance matrices [7], and propose a parameter splitting strategy handled with ADMM optimization algorithm [1] to solve the corresponding quartic optimization problem.

The rest of the paper is organized as follows. Section II gives the formulation of the problem with various regularizations. Section III presents the proposed method based on an ADMM optimization algorithm. Section IV illustrates and discusses the performance and runtime of the method. Concluding remarks are reported in Section V.

## II. PROBLEM DESCRIPTION

### A. Forward model and problem formulation

Given a bivariate signal $\mathbf{X} = [u[n], v[n]]_{1 \le n \le N} \in \mathbb{R}^{N \times 2}$ composed of $N$ samples, the observation model for the noisy signal $\mathbf{Y} = (y_d)_{1 \le d \le D} \in \mathbb{R}^{N \times D}$ measured in $D$ channels is

$$ \mathbf{Y} = \Phi(\mathbf{X}) + \varepsilon \in \mathbb{R}^{N \times D}, \quad (1) $$

where $\Phi$ is a known linear operator, and the noise $\varepsilon = [\varepsilon_d]_{1 \le d \le D}$ has stationary Gaussian components $\varepsilon_d \sim N(0, \Xi_d)$ with known covariance matrices $\Xi_d \in \mathbb{R}^{N \times N}$. For simplicity, one can consider $\Xi_d = I_N$. The generic linear transformation $\Phi$ appears in various applications where the bivariate signal is not directly observed and reconstructed from measurements over $D$ channels, such as gravitational signal processing [3]. To solve this reconstruction problem, one needs to impose certain priors on the reconstructed signal
X to obtain better reconstruction performance. These priors can be generically formulated as

$$ \mathbf{X}^* \in \arg \min_{\mathbf{X} \in \mathbb{R}^{N \times 2}} f(\mathbf{X}, \mathbf{Y}) + \lambda_1 g_1(\mathbf{X}) + \lambda_2 g_2(\mathbf{X}), \quad (2) $$

with data-fidelity term $f(\mathbf{X}, \mathbf{Y}) = \|\Phi(\mathbf{X}) - \mathbf{Y}\|_F^2$. The regularization functions $g_1$ and $g_2$ respectively reflect the prior on the 2 components of the bivariate signal and the geometric properties, with corresponding hyperparameters $\lambda_1$ and $\lambda_2$.

## B. Towards relevant regularizations

A typical prior on individual components of X is that each component has smooth variations over time [20], [11]. This is often modeled by penalizing the amplitude of the time gradient of X, which can be formulated $g_1(\mathbf{X}) = \lambda_1 \|\mathbf{D}\mathbf{X}\|_F^2$ with D the finite difference operator. The squared Frobenius norm reflects the assumption of Gaussian increments between two successive samples. This regularization leads to smooth evolution of each component, but does not ensure that both the signal and the polarization evolve smoothly over time. In bivariate signals, the polarization ellipse typically changes much slower than the frequency content of the signal [3]. Thus, the reconstructed signal is assumed to have a smooth polarization evolution. In this work, we propose to use the instantaneous covariance matrices as a representation of the polarization ellipse.

Starting from the analytic signal $\mathbf{X}_a = \mathbf{H}\mathbf{X} = [u_a[n], v_a[n]]_{n \in \{1, \dots, N\}} \in \mathbb{C}^{N \times 2}$ where $\mathbf{H} \in \mathbb{C}^{N \times N}$ filters out the negative frequencies [12], the instantaneous covariance matrix can be defined as

$$ \begin{aligned} \Sigma[n] &= \mathbf{X}_a[n]^\mathrm{H} \mathbf{X}_a[n] = \mathbf{X}^\mathrm{H} \mathbf{H}^\mathrm{H} \mathbf{e}_n \mathbf{e}_n^\mathrm{T} \mathbf{H} \mathbf{X} \\ &= \begin{bmatrix} |u_a[n]|^2 & u_a[n]^\mathrm{H} v_a[n] \\ v_a[n]^\mathrm{H} u_a[n] & |v_a[n]|^2 \end{bmatrix}, \end{aligned} \quad (3) $$

where $\mathbf{e}_n$ is the n-th canonical basis vector of $\mathbb{R}^N$. Here, $u_a[n]$ and $v_a[n]$ are the analytic versions of $u[n]$ and $v[n]$. The eigen-decomposition of $\Sigma[n]$ is parametrized by the instantaneous polarization ellipse parameters such as its size, eccentricity and orientation [7]. As a result, smoothing $\Sigma[n]$ leads to smooth evolution of the polarization ellipses through time.

To promote a smooth polarization evolution, we propose to penalize the evolution of the covariance matrices using finite differences with

$$ \begin{aligned} g_2(\mathbf{X}) &= \sum_{n=2}^{N} \|\Sigma[n] - \Sigma[n-1]\|_F^2 \\ &= \sum_{n=2}^{N} \|\mathbf{X}^\mathrm{T} \mathbf{H}^\mathrm{H} \mathbf{e}_n \mathbf{e}_n^\mathrm{T} \mathbf{H} \mathbf{X} - \mathbf{X}^\mathrm{T} \mathbf{H}^\mathrm{H} \mathbf{e}_{n-1} \mathbf{e}_{n-1}^\mathrm{T} \mathbf{H} \mathbf{X}\|_F^2 \\ &= \sum_{n=2}^{N} \|\mathbf{X}^\mathrm{T} \mathbf{H}^\mathrm{H} \mathbf{J}_n \mathbf{H} \mathbf{X}\|_F^2, \end{aligned} \quad (4) $$

where $\mathbf{J}_n = (\mathbf{e}_n \mathbf{e}_n^\mathrm{T} - \mathbf{e}_{n-1} \mathbf{e}_{n-1}^\mathrm{T})$. The Euclidean distance does not necessarily take the geometry of covariance matrices into account, yet remains as a good approximation locally [16].

The proposed method, called Time And COvariance Smoothing (TACOS), consists in solving the optimization problem in (2) with these proposals for $g_1$ and $g_2$. This problem is differentiable, quartic due to the covariance regularization term $g_2$ so potentially not convex, and does not admit a closed-form solution. In the next section, a parameter splitting strategy is used to solve the problem while exploiting its structure.

## III. PROPOSED ALGORITHM FOR TACOS

### A. Parameter splitting and ADMM

We propose to apply parameter splitting to the quartic term (4) so that alternating the optimization between the primal variable and a splitting variable involves easier-to-solve convex sub-problems. To do so, we define another parameter $\mathbf{Z} = \mathbf{H}\mathbf{X}$ and rewrite

$$ \mathbf{X}^* \in \arg \min_{\substack{\mathbf{X} \in \mathbb{R}^{N \times 2} \\ \mathbf{Z} = \mathbf{H}\mathbf{X} \in \mathbb{C}^{N \times 2}}} f(\mathbf{X}, \mathbf{Y}) + \lambda_1 g_1(\mathbf{X}) + \lambda_2 \tilde{g}_2(\mathbf{X}, \mathbf{Z}), \quad (5) $$

where $\tilde{g}_2(\mathbf{X}, \mathbf{Z}) = \sum_{n=2}^N \|\mathbf{X}^\mathrm{T} \mathbf{H}^\mathrm{H} \mathbf{J}_n \mathbf{Z}\|_F^2$. The problem (5) can be addressed with an ADMM algorithm [1]. The constraint $\mathbf{Z} = \mathbf{H}\mathbf{X}$ is imposed by minimizing

$$ L(\mathbf{X}, \mathbf{Z}, \mathbf{U}) = f(\mathbf{X}, \mathbf{Y}) + \lambda_1 g_1(\mathbf{X}) + \lambda_2 \tilde{g}_2(\mathbf{X}, \mathbf{Z}) \\ + \frac{\rho}{2} (\|\mathbf{H}\mathbf{X} - \mathbf{Z} + \mathbf{U}\|_F^2 - \|\mathbf{U}\|_F^2), \quad (6) $$

where $\rho > 0$ adjusts the balance between the constraint and minimization of the objective function. Each iteration $l \in \mathbb{N}$ of ADMM is defined by the update steps

$$ \mathbf{X}^{(l+1)} = \arg \min_{\mathbf{X} \in \mathbb{R}^{N \times 2}} L(\mathbf{X}, \mathbf{Z}^{(l)}, \mathbf{U}^{(l)}) \quad (7a) $$

$$ \mathbf{Z}^{(l+1)} = \arg \min_{\mathbf{Z} \in \mathbb{C}^{N \times 2}} L(\mathbf{X}^{(l+1)}, \mathbf{Z}, \mathbf{U}^{(l)}) \quad (7b) $$

$$ \mathbf{U}^{(l+1)} = \mathbf{U}^{(l)} + \rho(\mathbf{H}\mathbf{X}^{(l+1)} - \mathbf{Z}^{(l+1)}). \quad (7c) $$

The typical convergence criteria for such an algorithm are the primal or the dual error being under a certain threshold [1]. The parameter $\rho > 0$ can have a significant effect on the convergence of this algorithm. It can be chosen either empirically or with automated strategies [1], [6], which goes beyond the scope of this work.

### B. Detailed update steps

To illustrate each update steps, we consider the following generic linear model

$$ \forall d \in \{1, \dots, D\}, \quad \mathbf{y}_d = \mathbf{T}_d \mathbf{X} \mathbf{r}_d + \varepsilon_d, \quad (8) $$

that is, the observation operator $\Phi$ in (1) takes the form $\Phi(\mathbf{X}) = (\mathbf{T}_d \mathbf{X} \mathbf{r}_d)_{1 \le d \le D}$, with $\mathbf{T}_d$ and $\mathbf{r}_d$ generic left-hand and right-hand side transformations. Many applications fall under this generic linear scheme, notably the reconstruction of gravitational waves from interferometric detectors [3].