Je Tiens D'abord À Remercier

Eric Chassande-Mottin

Pierre Chainais

Astrid Lamberts

Florent Robinet

Natalia Korsakova

Keywords: Gravitational waves, Galactic binaries, Sparse modeling, Missing data, Inpainting, Learned representation, Auto-encoder 1. c

de me faire l'honneur de participer à mon jury de thèse. Merci plus particulièrement à Eric Chassande-Mottin et Pierre Chainais d'avoir accepté la lourde tâche de rapporteur.

J'aimerais exprimer toute ma reconnaissance à mes deux directeurs de thèse Hervé Moutarde et Jérôme Bobin. Merci de m'avoir forcée à aller en conférence, à rédiger des articles, de m'avoir amenée à me poser les bonnes questions et de m'avoir appris ce qu'est la recherche. Les trois ans que j'ai passé à travailler avec vous m'ont énormément appris, que ce soit sur le plan scientifique ou humain. Un grand merci aussi pour le soutien dont vous avez fait preuve pendant que je rédigeais mon manuscrit, que ce soit par vos retours ou pour le chocolat que vous m'avez offert dans les moments difficiles.

Je tiens aussi à remercier de manière plus générale le Département de Physique Nucléaire du CEA Saclay sous la direction de Franck Sabatié et la codirection de Christophe Theisen, au sein duquel j'ai effectué mon travail de thèse. Vous m'avez appris que les échanges entre des communautés très différentes sont possibles et peuvent être très porteurs si l'on garde l'esprit ouvert ! Je vous suis aussi très reconnaissante pour m'avoir inculqué les bases de la physique hadronique -dont j'ignorais tout jusqu'à mon séjour dans ces locaux. Merci à Nicolas Pailley d'avoir assuré le support technique quand j'ai pu en avoir besoin. Les secrétaires du département Danielle Coret et Isabelle Richard ont joué un rôle central dans le bon déroulement de ma thèse au quotidien -toutes les démarches auraient été beaucoup plus compliquées sans votre aide précieuse ! Un grand merci à la communauté d'analyse de données LISA, au sein de laquelle j'ai présenté mes travaux à diverses reprises. En particulier, Antoine Petiteau, Quentin Baghi, Natalia Korsakova et Stas Babak ont, à travers les échanges que j'ai pu avoir avec eux, grandement contribué à enrichir ma thèse. Merci à Thierry Foglizzo d'avoir organisé le journal club autours de LISA ainsi que pour les discussions que nous avons pu avoir.

A Hervé Dutrieux, j'exprime ma gratitude pour les trois ans pendant lesquels nous avons partagé le même bureau. Cette collaboration nous a, je pense, beaucoup apportée à tous les deux

Résumé

Le futur observatoire spatial d'ondes gravitationnelles LISA ouvrira une nouvelle fenêtre de mesure des ondes gravitationnelles, permettant d'observer des systèmes émetteurs difficilement visibles avec les observatoires terrestres actuels.

Parmi ces sources, la détection de binaires galactiques comme sources d'ondes gravitationnelles promet une richesse d'informations sans précédent sur ces systèmes, mais soulève également plusieurs défis dans le traitement du signal. En particulier, le grand nombre de sources attendues et la présence à la fois de bruit instrumental complexe et d'artefacts entachant les données nécessitent le développement de méthodes robustes.

Nous décrivons ici une méthode originale de reconstruction non paramétrique de l'empreinte des binaires galactiques dans les mesures affectées par le bruit instrumental typique de l'observatoire spatial d'ondes gravitationnelles LISA. Cette méthode, basée sur une procédure de débruitage, vise à séparer du bruit la somme de tous les signaux provenant des binaires galactiques. Ses performances sont solidement évaluées et son applicabilité est démontrée avec succès sur un exemple simple impliquant des binaires galactiques de vérification récemment proposées dans les défis d'analyse de données LISA (LDC 1-3).

Nous montrons ensuite comment tirer le meilleur parti de ce modèle simple pour atténuer l'impact des interruptions de prises de données. Qu'elles soient prévues ou non, les données manquantes réduisent notre capacité à détecter les signaux faibles comme les binaires galactiques et augmentent le risque de mauvaise détection et d'identification biaisée. Inspirés par les progrès du traitement du signal, nous introduisons un algorithme d'inpainting non paramétrique basé sur la représentation parcimonieuse du signal binaire galactique dans le domaine de Fourier. Contrairement aux approches traditionnelles d'inpainting, la statistique du bruit est théoriquement connue uniquement sur des mesures non lacunaires. Il est donc nécessaire de récupérer conjointement le bruit et le signal binaire galactique. Ce processus limite les fuites de bruit et récupère avec précision le signal manquant, au point de pouvoir débiaiser un processus d'identification mené sur les données corrigées, même lorsque jusqu'à environ 30% des données sont perdues. En outre, nous soulignons que l'algorithme proposé produit une estimation statistiquement cohérente du bruit sans lacune. Nous évaluons ensuite les performances des méthodes d'inpainting proposées pour récupérer le signal des ondes gravitationnelles sur le défi LDC1-3.

Enfin, une caractérisation plus robuste des signaux binaires galactiques nécessite une modélisation plus précise de ces signaux. Nous proposons une approche nonparamétrique basée sur l'apprentissage pour la modélisation des signaux binaires galactiques, ce qui donne lieu à plusieurs prometteuses applications à LISA. Parmi cellesci, on trouve la représentation du signal dans une base restreinte et de faible dimension, mais aussi le débruitage du signal, la détection des binaires galactiques, la séparation des sources et l'estimation partielle rapide des paramètres. Les performances du modèle sont évaluées sur ces différentes applications.

Mots clés :

Ondes gravitationnelles, Binaires galactiques, Modélisation parcimonieuse, Données manquantes, Inpainting, Représentations apprises, Auto-encodeur

Introduction

This thesis deals with the characterization of galactic binary systems by gravitational waves in the framework of the Laser Space Interferometer Antenna (LISA) mission. LISA is a space-based gravitational wave detector scheduled for launch in 2034. By its design it will be able to detect gravitational wave sources of much lower frequencies than terrestrial observatories. These include all the galactic binaries in the galaxy -a few million systems. However, LISA has only a limited number of information channels and will have to limit itself to measuring the sum of these signals. In this context, characterizing them individually is a real challenge, in particular because other sources of gravitational waves are also present, but also because LISA is subject to a complex instrumental noise and presents many artifacts (interruptions of data taking, glitches).

LISA data analysis must take into account the characteristics of noise and artifacts and prove that it is robust against them. It is in this context that I propose to introduce non-parametric methods of data analysis ; these methods, if they do not allow the direct identification of the systems, are very complementary with the Bayesian approach which is most represented today within the consortium. This is a framework that also allows to integrate easily the presence of noise and artifacts in the data.

In a first step, I show how the classical framework of sparse signal representation can be adapted to the gravitational wave sources we are looking for. Based on the simple assumption that the signals of binaries are quasi periodic and therefore sparse in the Fourier basis, I show that we can have a good detection and reconstruction of these signals. We easily include the presence of a colored Gaussian noise in the modeling of the data. This first model allowed us to submit a solution to the LISA Data Challenge 1-3 during which we detected all the sources, and only the injected sources.

I have pushed the capabilities of this model further by showing the benefits it could bring when interruptions of data taking occur. This particular case can be directly related to compressed sensing methods, which I used to develop algorithms to compensate for the impact of interruptions on the data -and on the parametric analysis that can be done with it. These algorithms are able to cope with losses of up to nearly 30% of the data -which could be a realistic order of magnitude for the mission.

Finally, I propose an alternative approach to sparse modeling that allows to push the process of galactic binary signal characterization much further. This new approach, based on learning, allows to accurately represent the galactic binary signals. Such a model shows very good performances in terms of detection capacity, but also signal estimation and source separation. It presents interesting possibilities for the LISA mission because it relies heavily on the compressibility of the signal and on an underlying low dimensionality representation. I have conducted an extensive benchmark of its performances and propose many interesting avenues in the context of LISA data analysis.

Chapter 2

The two (or more)-body problem in gravitation : from Newton to gravitational waves

From Newton to Einstein to gravitational waves

Gravity is a key phenomenon to understand the world around us. Even if it is easily observable, it took centuries to reach the current description that we make of it. The development of General Relativity -the theory of Gravitation -in the course of the 20th century has literally transformed our vision of the world and has not only contributed to explain many phenomena but also led to many discoveries.

Where Einstein finds a solution to Newton's model limits

Newton theory of gravity has long been the best to describe the phenomenon of gravitational attraction as the equations were simple and it could describe very complex behaviors, such as Kepler's orbits, ballistics or even the three-body problem.

In this representation, gravitation is seen as a conservative central force (i.e. deriving from a central potential) -→ F A→B that a body A with gravitational mass m A exerts on a body B with gravitational mass m B :

- → F A→B (t) = -G • m A m B -→ AB(t) 2 2 • -→ AB(t) -→ AB(t) 2 .
(2.1.1)

This force derives from the central potential generated by the mass distribution ρ describing the massive bodies :

∇ 2 Φ = 4πGρ . (2.1.2)
The Poisson field equation (2.1.2) depends linearly on the mass distribution ρ [1].

However, in 1859 an important deviation from the model was observed when Mercury's expected perihelion did not match with the observation [2] by an important margin.

An other phenomenon also challenged this representation : with the rise of atomic clocks, precise measurements showed that clocks at different gravitational potential were experiencing different time flows.

In fact, Newton's theory was relying on several implicit axioms : the existence of an absolute time, the existence of Newtonian potential (2.1.2), the motion equation that ensues for massive bodies and the equality of inertial mass and gravitational mass.

Moreover, a byproduct of Eq.(2.1.2) is that any movement or change of body A is instantaneously propagated to body B with infinite speed [1] : there is no delay between the movement of A and its effect on B.

Einstein designed the theory of general relativity with these elements in mind. Changing the spacetime geometry was the simplest way to explain the clock deregulation depending on the gravity potential. One element he could still not explain and that was considered as an axiom of his theory is the apparent equality between inertial mass and gravitational mass.

Equivalence principle : (also called "Principle of Relativity") : "Experiments in a sufficiently small freely falling laboratory, over a sufficiently short time, give results that are indistinguishable from those of the same experiments in an inertial frame in empty space." [3]).

The very principle of General Relativity is to consider both time and space (i.e. spacetime) as parameters of a geometry which evolution is governed by mass distribution :

Spacetime curved geometry : "Space acts on matter, telling it how to move. In turn, matter reacts back on space, telling it how to curve" [1].

Based on these principles, the new unknown that has to be understood to describe the universe is its geometry, and more precisely its metric. Key element of this theory is Einstein's field equation :

G αβ = 8πT αβ .
(2.1.3) (stated here under classical convention c = 1, G = 1 1 ).

It explains how matter through the energy-momentum tensor T αβ acts on spacetime curvature and metric, as Einstein tensor G αβ can be expressed as a non-linear function of the metric g αβ and its derivatives of order 1 and 2 [1; 3].

Analogy with the electro-magnetic field : Einstein's equation Eq.(2.1.3) describes the behavior of spacetime curvature as a function of a source term T αβ . It is really similar to Maxwell's electro-magnetic (EM) field equations. In EM, a charged particle moving in a field is subjected to a force (the Lorentz force) which binds its trajectory. Then retroactively the charged particle and its movement act as a source term (charge and current density) for Maxwell's equations, and thus result in variations of the EM field.

Einstein's equations abide by a similar phenomenon. A free-falling massive body's movement is subject to gravitational field through spacetime metric and curvature that bind its trajectory along the spacetime geodesics. Retroactively, the massive body acts as a source term for Einstein's equation and thus plays a role in spacetime's curvature and metric variations.

Gravitational waves

As previously stated, Einstein's equation Eq. (2.1.3) is a non-linear equation which unknown is the spacetime metric g αβ . Confronted with such a complex equation, scientists adopted two strategies. The first strategy was to look for explicit solutions of the problem with a very simple source term : this resulted in Schwarzschild geometry ( [1], chapter 31) and Kerr-Newman geometry [1], chapter 33). The second one was to try to linearize Einstein's field equation. But linearize with respect to what ? And under what hypothesis ?

As the metric g αβ is the unknown, the linearization was done under the assumption that far from any matter source, this metric should be close to Minkovski's metric η αβ . This is directly justified by the equivalence principle ; this type of development is called "Post Minkovskian" 2 .

The linearized metric writes :

g αβ spacetime metric
= η αβ flat metric + h αβ small variation , (2. 1.4) with |h αβ | << 1 (the strain). Using Lorentz Gauge (thus performing the corresponding change of variable) :

h αβ = h αβ - 1 2 η αβ h , (2.1.5)
with h = h α α ,

(2.1.6) using Einstein's summation convention. Performing this change in Eq. (2.1.3), developing and keeping only first order terms with regard to h αβ , we obtain the following equation :

□h αβ = -16πT αβ . (2.1.7)
This is a wave equation with source term T αβ that admits solutions called gravitational waves (GWs). Still, they were obtained after so many approximations that the whole scientific community wondered : are these waves real ? Or are they only a mathematical artefact coming from the way the theory was stated ? These questions literally split the scientific community for decades between those who believed in their existence and those who did not -even Einstein was doubtful toward this result, as too many assumptions were made to obtain a proper, simple linearized equation.

It is only after their first indirect detection in 1982 (Taylor, Weisberg and Hulse linked a a pulsar orbit energy loss to the energy of the GWs they hypothetically emitted [4]) that this was accepted as one of the outcomes -and a proof of validity -of general relativity. GW first direct detection in 2015 by LIGO-Virgo collaboration (LVC) [5] then opened a new window to probe the universe around us.

Gravitational waves signal in a very simplified case : the galactic binaries

In this section, we show how to compute a first solution of Eq.(2.1.7) based on several more approximations : the GW source moves slowly (long wavelengths), GW amplitude is weak and the observer is far away from the source. Even if the solution is not exact, it is good enough so that we can rely on it to give several semi-quantitative properties of the expected GW signals.

The computations presented here follow the argument of [3] (chapters 22 and 23). We tried to detail them as much as possible to to make the text self-contained. More than an exact result, what we want to outline here is the physical origin of the phenomenon and how it affects its surroundings. GW equation Eq.(2.1.7) is a traditional wave equation with a source term. Using Green functions we know that the general solution at position -→ x of space and time t can be expressed as an integral of the delayed source (here the energy-momentum tensor T αβ ) :

h αβ (t, - → x ) = 4 d 3 x ′ T αβ (t -- → x - - → x ′ 2 , - → x ′ ) - → x - - → x ′ 2 ,
(2.2.1) using the convention c = 1. In absence of other interaction than gravitation and under the assumption that the source has a low velocity (compared to the velocity of light), the energy-momentum tensor is dominated by the rest mass density µ [3] :

T tt = µ(t, - → x ) , (2.2.2)
where µ(t, -→ x ) gives the mass density at time t and position -→ x of space for the system under study. Let us consider the case of compact binary stars : two bodies with identical shapes and masses m = M are rotating one around another with null eccentricity.

Choosing the origin of space coordinates as the center of mass of the system and an orthogonal coordinate system (x, y, z) in which the compact binary stars orbit in the (xy)-plane, one of the bodies has time coordinates :

     x(t) = R cos(Ωt) y(t) = R sin(Ωt) z(t) = 0 , (2.2.3)
the other body having opposite coordinates. We understand easily that the mass density here is only made of the two bodies (and their position in space at time t) : µ has finite support in space. Thus, the integral described in Eq.(2.2.1) also has finite support for -→ x ′ within the binary plausible positions. For an observer at position -→ x far from the source (i.e. ∥ -→ x ∥ = r >> R), we have -→ x --→ x ′ 2 ≃ r. Consequently :

h αβ (t, - → x ) ----→ r→+∞ 4 r d 3 x ′ T αβ (t -r, - → x ′ ) .
(2.2.4)

Now, taking into account the energy conservation for the energy-impulsion tensor for flat spacetime (which is, to linear order, equal to our spacetime), we have :

∂T αβ ∂x β = 0 , (2.2.5)
for any α ∈ {t, 1, 2, 3}. This is similar to the principle of charge conservation within a closed system, but generalized to the energy-momentum tensor in a flat spacetime.

Considering first :

∂T tβ ∂x β = 0 , (2.2.6)
and then time-derivating, yields :

∂ 2 T tt ∂t 2 = - j∈{1,2,3} ∂ ∂x j
∂T tj ∂t .

(2.2.7)

Now, the energy-momentum tensor being symmetric, we can swap the indices : T tk = T kt . Applying one more time the conservation law Eq. (2.2.5), we get :

∂ 2 T tt ∂t 2 = i∈{1,2,3} j∈{1,2,3} ∂ 2 T ij ∂x i ∂x j .
(2.2.8)

We multiply both sides of the equation by x k x ℓ and we integrate by part :

1 2 d 2 dt 2 d 3 x ′ x ′k x ′ℓ T tt (t, - → x ′ ) = d 3 x ′ T kℓ (t, - → x ′ ) .
(2.2.9)

We recognize here the second mass moment of classical mechanics :

I kℓ (t) = d 3 x ′ x ′k x ′ℓ T tt (t, - → x ′ ) = d 3 x ′ x ′k x ′ℓ µ(t, - → x ′ ) .
(2.2.10)

and get the final formula :

h kℓ (t, - → x ) ----→ r→+∞ 2 r d 2 dt 2 I kℓ (t -r) . (2.2.11)
Before going further, we can notice two important points : first, GW emission is directly related to the mass distribution and its trajectory with time. Second, as the term r -t only marks a delayed response and thus has no impact on the signal's amplitude, the GW amplitude decreases as 1/r. In comparison, EM waves decrease as 1/r 2 : for the same energy density of emission, an observer at distance r from the emitting body will see that GW are less attenuated than their EM counterparts.

We can push a bit further the computations to get a more precise modeling of the signal. For this, we use the trajectory of the binary system presented in Eq.(2.2.3). Let 2M = M + M be the total mass of the system (i.e. the sum of the two body masses).

Using Eq.(2.2.10), the second mass moments (taking into account the two stars) are :

     I xx (t) = 2M R 2 cos 2 (Ωt) = M R 2 (1 + cos(2Ωt)) I xy (t) = 2M R 2 cos(Ωt) sin(Ωt) = M R 2 sin(2Ωt) I yy (t) = 2M R 2 sin 2 (Ωt) = M R 2 (1 -cos(2Ωt)) (2.2.12)
All other momenta are zero.

h kℓ (t, - → x ) ----→ r→+∞ - 8Ω 2 M R 2 r      cos 2Ω(t -r) sin 2Ω(t -r) 0 sin 2Ω(t -r) -cos 2Ω(t -r) 0 0 0 0      (2.2.13)
Even if this formula is approximate, it provides insightful features on the GW signal emitted by a two-body system. Firstly, we get a GW frequency that is twice the orbital frequency of the emitting bodies. This is still verified in more precise computations. The amplitude of the GW decreases as 1/r. Secondly, we can now compute an order of magnitude of the wave's amplitude, as it is here defined by h 0 = 8Ω 2 M R 2 /r. Using data with representative order of magnitude for galactic binary systems, we get :

8Ω 2 M R 2 r ∼ 10 -21 .
(2.2.14)

GW are ripples propagating through spacetime. These ripples directly affect the spacetime metric and can be measured. They originate from an asymmetry of the quadrupolar moment, meaning that any system breaking spherical symmetry can emit GWs. Among these systems, there are of course the binary systems, but also others like the asymmetric explosions of supernovae. Similar computations can be performed to determine the signal they emit, but they require different approximations taking into account the system's specifics. However, spacetime is very rigid : these perturbations have a really low amplitude.

Motivations for GW observation

GW measurement is a new information channel that can be used to probe the Universe. It will critically complement the EM observations already in use in Astrophysics and Cosmology. Among the many elements motivating the observation of GWs, let us list the following :

Multi-messenger astronomy : using GWs, we can detect in advance signals predicting events with important EM emissions, among which we find some black hole mergers and supernovae. We can even use GWs to spot systems with a lower imprint in the visible spectrum, such as binary neutron stars. Moreover, we can organise space GW detector -ground GW detector collaborations to observe on longer runs events which frequency increases with time. Moreover, there will be synergies between all detector -both ground and space based, detecting GW or EM waves -such as with SKA [6], the Athena mission [7], LVC and the LISA mission [8; 9].

A window to unexplored time :

As we explained in Section 2.2, GW dampening is slower than EM wave damping. We could potentially observe phenomena that occur everywhere in the universe. Not only could we look further, we could also look "earlier" : no EM observations can be done before the apparition of the cosmic microwave background about 380.000 years after Big Bang. However, there is a possibility that GW signals emitted before that moment could be seen [10].

A new description of our galaxy : some objects are not (or hardly) EM detectable but could be GW detectable, like galactic binaries. Some objects are hard to detect using EM waves if we are not already aware of their presence. GW detection can improve cosmological models : the LVC has already started to characterize the population of compact binary objects emitting GWs in our Galaxy [11].

Discovery potential :

We could measure new types of signals corresponding to new physics. We could also probe several theoretical extensions of GR or perform precision tests of GR in the strong coupling regime [10; 12]. (source : [13]) Right : Test mass motion in spacetime in absence of GWs corresponding to a flat spacetime (orange, labelled as "Flat") and when a GW is passing (blue, labelled as "GW")

Gravitational waves observation : how to observe the non-visible ?

GWs may propagate far, but they still have an amplitude so low that it is a real challenge to measure their imprint. Many experiments were designed to directly detect them [13], but only few are actually able to do so.

Indirect detection

The existence of GWs was first confirmed in 1982. In their article [4], Taylor and Weisberg explain that the decay of the orbital period of the binary pulsar PSR 1913+16 can exactly be explained by the quantity of energy radiated through GWs. The Nobel Prize in Physics of year 1993 was dedicated to this discovery and the possibilities it brought to science.

A few decades later, the Pulsar Timing Array (PTA) [14] method is now presented to the community as a "natural detector" for GW. It relies on the idea that pulsar emit extremely periodic EM pulses than can be observed with ground-based detectors like SKA. Their time of arrival could be perturbed by incoming GW signals. By looking at correlations between arrival delays of the pulses coming from different pulsars, scientists expect to detect GW signatures indirectly. This work is still ongoing [15].

Direct detection

The first successful direct detection of GWs was realized in 2015 within the LVC collaboration, and was awarded the Nobel Prize in Physics of year 2017. The detectors of the collaboration are based on a laser Michelson interferometer [16] measuring with great precision the distance between two "free-falling" test masses separated by a wide distance, as shown in Fig. 2.1 (left). The test masses are free from any interaction except from gravity, and a distance variation between the two can be interpreted as the effect of a passing GW. Fig. 2.1 (right) shows the variation of distance between two test masses A and B (taking the position of A as a reference) in a flat spacetime when a GW passes.

The amplitude of the variation remains very small : for a wave propagating along the z-axis, we can estimate the variation of the distance between two test masses separated by L * in a flat spacetime. Taking one of the test masses as a reference, and considering the other one along the x-axis, we get ( [3], chapter 16) :

L(t) = L * 0 dx 1 + h xx (t, 0) 1/2 ≃ L * 1 + 1 2 h xx (t, 0) . (2.4.1)
Thus, writing L(t) = L * + δL(t) :

δL L * = 1 2 h xx (t, 0) ∼ 10 -21 , (2.4.2)
using the order of magnitude Eq.(2.2.14) : the wider the arms of the interferometer, the greater the distance variation δL that should be measured.

For LIGO, the armlength is 4 km (the effective armlength is wider thanks to a complex assembly of mirrors increasing the laser beam optical path length), thus can only detect phenomena with amplitude that is high (for instance, binary black holes mergers). The detection capacity of this detector is limited in lower frequencies by the Earth seismic noise. The detected sources have a frequency around 10 2 Hz restricting the observation to certain types of sources [17].

Having an observatory with wider arms and freed from seismic noise would open many doors in terms of observation. The joint ESA-NASA project "Laser Interferometer Space Antenna" (LISA for short) is an answer to both requirements, as it will be an interferometer in space with armlength of 2.5 million km.

Chapter 3

The LISA project LISA is a space-based GW observatory which has been developed since 1997 as a collaboration between ESA and NASA. The project experienced a critical phase when NASA retracted from it in 2011, optimistic prospects were simultaneously triggered by the first direct detection of GW in 2015 by LVC (announced at the beginning of the year 2016) and the success of the demonstrator LISA PathFinder (LPF) in december 2015.

LPF was a demonstrator aiming at validating the technologies that will be used during LISA mission. The results were impressive, as final noise levels were much lower than the requirements. Following these two events, NASA re-entered the project in 2016. Finally, in June 2017 LISA was selected as a ESA L3 mission ; phase A was completed in December 2021 and we are now in phase B1 that will last until mission validation mid 2024.

The following description of the LISA project is data-analysis oriented, and thus is non-exhaustive. More details about the mission can be found in [1; 2]

LISA Objectives and Scientific objectives

LISA measurements are based on the principle of interferometry that described in Section 2.4.2. Nonetheless, LISA's armlengths were desigend much wider than the ones of LVC to observe ranges of phenomenon inaccessible to the LVC. Fig. 3.1 displays the known GW sources with regard to their expected emission frequency range. Some events, like merging black hole binaries or supermassive black holes, are visible by both LISA and Earth-based GW observatories. As many of these events have EM counterparts, this should lead to fruitful cooperations with the different kinds of observatories to conduct multimessenger astronomy [3].

A detailed description of sources and associated scientific objectives can be found in [1; 4; 5]. The main observable sources are :

Massive Black Holes Binaries (MBHBs) : Two massive black holes rotating one around another during in-spiral phase, ending up merging. The amplitude of the expected transient signal can become strong at merging time, making it one of the most visible sources for this mission.

Galactic Binaries (GBs) :

Two bodies with similar masses rotating one around another, far enough from merging time producing a low amplitude, periodic continuous signal. It includes white dwarfs, neutron stars but also stellar-origin black holes. Tens of millions of such systems are present in our Galaxy ; among them, it should be possible to identify the sky position and the physical characteristics of about 20.000 of them. The GW signals coming from unidentified systems will form a non-negligible background noise (also called "Confusion Noise"). A more precise description can be found in [6].

Extreme Mass Ratio In-spirals (EMRI) :

A very light body orbiting around a heavy body. After a long in-spiral phase where the two bodies get closer, the light body should merge with the heavy body. These systems are expected to produce complex signals that are hard to model. Nonetheless, they could help putting to test general relativity as the movement of the light body gives a precise cartography of the gravity field around the heavy body. It is hard to estimate the detectability rate of these sources -between 0 and 1, 000 of such events could happen during the whole mission.

Stellar Origin Black Holes (SOBH) :

They are the sources the most detected by LVC. They should also be visible by LISA in their early in-spiral phase. This calls for synergies between the two devices : LISA could detect the merger long before it happens and produce a warning to alert LVC that an event is likely to happen in a near future.

Stochastic Background : Of cosmological origin, this background was created in the early universe. Its detection could help understand better the changes that happen at that time and give new directions regarding cosmological phenomena that happened in the early Universe. It is one of the few experimental possibilities to probe this epoch in the history of the Universe.

GW bursts and unforeseen sources :

As the first survey observatory scanning the gravitational universe in this frequency range, LISA may detect unforeseen sources and trigger unexpected exciting developments in physics.

Emitted GW signals depend on the sources physical characteristics. If we can measure the signals, we might be able to trace the characteristics of the bodies, such as their positions in the Galaxy, their orientation and their masses. LISA greatly differs from the LVC observatory not only in the frequency range of observation, but also in the amount of sources that are expected to be seen. When the LVC reports detecting about 1 source per week [7], LISA will observe thousands of continuous phenomena superposed to transients GW bursts similar to that of LVC observations. Many of these sources produce an EM counterpart at critical moments of their life that are often preceded by GW emission over long time periods. LISA could play the role of a trigger to forecast EM bursts visible by EM observatories. LISA is a constellation of three satellites separated by 2.5 million km one from another. As shown in Fig. 3.2 (left), this constellation rotates around its center while following the Earth's orbit around the Sun. The three satellites play symmetrical roles, and have the same structure illustrated in Fig. 3.2 (right).

Measurement principle : Time Delay Interferometry

One satellites is made of two main parts [2] : the interferometer measurement system and the disturbance reduction system -insuring the stability of the system even in presence of disturbances. We introduce here the measurement principle, and we tackle the noise questions in the next section. 1 and the test-masses movements ⃗ δ 1 , ⃗ δ ′ 1 must be taken into account. The real system relies on more than one laser per optical bench and a more complex optical path. (Source : [9])

Interferometric measurement system

LISA is a constellation of three satellites, also called "spacecrafts". Each spacecraft (S/C) includes two optical benches. Each one is composed of a laser, a test mass and two photodetectors [9] performing various interferometric measurements between the local and incoming laser beams. The test masses are continuously free falling along a geodesic of the gravitational field. Thus, any distance variation between the two masses can only be interpreted as the passage of a GW. Fig. 3.4 shows the effect of a passing GW on the LISA constellation for the two polarizations h + and h × : the variation observed for each arm is different and depends on the polarization of the incoming GW.

On one optical bench, three interferometric measurements based on Doppler effect are done [2] : a laser beam is emitted by a spacecraft at a nominal frequency, and any passing GW signal modulates the frequency of this beam. This modulation is equivalent to a Doppler shift of the laser frequency ; it can be measured by comparing the modulated beam with a laser at nominal frequency. However, this calls for a precise control of lasers nominal frequency that is hard to reach because many non-coherent laser sources are used.

To palliate this issue, we compare through interferometric measurements the lasers emission frequencies with a reference. The Main Interferometer compares the received laser beam frequency to the bench reference laser frequency (serving as a local oscillator). The Test Mass Interferometer compares the transmitted laser beam frequency to the bench reference laser frequency. The Reference Interferometer compares the transmitted laser beam frequency to the other bench reference laser frequency. A simplified version of the assembly is shown in Fig. 3.3 : the two optical benches perform interferometric measurements between the incoming laser beam and the local laser. Interferometric measurements are also realised between the two benches. The specificity of LISA as an interferometer comes from the way the measurement itself is done. Indeed, this interferometer has three arms that are unusually wide with variable lengths. On the one hand, the inter-satellite distance being 2.5 million km, the received light power is about a few hundreds pW whereas it was about 1 W at emission [2]. As only few photons reach the other spacecraft after being emitted, this prevents the light beam from going back and forth between two spacecrafts : the beam power would not be bright enough to be detected.

One laser beam can only go from one spacecraft to another : this makes it unavoidable to use at least one laser per spacecraft. But knowing that the constellation rotates around its center, the Sagnac effect must be accounted for. The propagation time of light between two satellites will not be the same if the light propagates in the direction of rotation or in the opposite direction. Thus, between spacecraft i and spa-cecraft j, we must measure the propagation time of light from i to j, but also from j to i. Thus, there is a total of 6 lasers exchanges. These lasers were built to emit beams around a nominal frequency ν 0 . However, each of them produces independently its own fluctuations around the nominal frequency. This is called laser noise -it has an amplitude that is 8 orders of magnitude above any detectable GW signal [9]. Furthermore the distance between satellites (i.e. the interferometer armlength) vary with time and differ from one arm to another as the satellites are not rigidly linked together. Satellites trajectories can deviate from the geodesic they follow because of external events (solar wind, etc). Even without these external solicitations, the theoretical trajectory along a geodesic implies for the inter-satellites distance to change slowly over time. This create an additional effect that must be considered to obtain precise interferometry. This does not fit Michelson's traditional frame : classically, equal armlength ensures laser noise suppression as only one laser is present. Therefore, it was necessary to find a way to remove the laser noise while maintaining the signal of interest.

Time Delay Interferometry

The Time Delay Interferometry (TDI) [9] was the first method developed to deal with laser noise. We will provide here the simplest example given in [9] as we believe this is enough to understand its underlying principle. More details can be found in [9] and references therein. Fig. 3.5 shows a simplified layout for a two-arm Michelson interferometer with different armlengths L 1 , L 2 . The laser emits a beam at frequency ν 0 with time fluctuations C(t) that we want to suppress. Given two interferometric measurements y 1 (t), y 2 (t), and under the hypothesis that we know exactly the armlengths L 1 and L 2 , it is possible to completely remove the laser noise.

The signal y i (t) measured at time t measures the difference between the laser frequency at time t and the incoming laser beam (emitted by the same laser at time t -2L i /c). The laser frequency is defined by ν 0 + C(t). Meanwhile, the laser beam frequency that was emitted at frequency ν 0 + C(t -2L i /c) is tainted by GW signature h i (t) and residual noise n i (t). Therefore, the interferometric measurements write :

y 1 (t) = ν 0 + C t -2 L 1 c + h 1 (t) + n 1 (t) -ν 0 + C(t) , y 2 (t) = ν 0 + C t -2 L 2 c + h 2 (t) + n 2 (t) -ν 0 + C(t) . (3.2.1)
The same laser fluctuation C(t) is measured twice through y 1 and y 2 . However, the delayed fluctuation C(t-L i /c) also appears in the measured signal : simply subtracting one measurement to the other is not enough to cancell all the laser noise. Yet, if instead we build a combination of the two measurements at different times :

X ≡ [y 1 (t) -y 2 (t)] -y 1 t -2 L 2 c -y 2 t -2 L 1 c . (3.2.2)
The laser noise is cancelled and we get :

X = h 1 (t) -h 1 t -2 L 2 c -h 2 (t) -h 2 t -2 L 1 c Signal + n 1 (t) -n 1 t -2 L 2 c -n 2 (t) -n 2 t -2 L 1 c Residual Noise . (3.2.3)
The signal is still there in a form mixing the signals present on the two arms taken at different times. This new GW signature is now detectable if its amplitude is high enough to surpass that of the residual noise.

The case of LISA is much more complicated, as we have 6 laser beams and thus 6 independent laser noises to cancel using all the interferometric measurements. Moreover, the armlengths are unknowns that also have to be estimated in the process. In the end, in a similar fashion, we obtain three "pseudo-Michelson" measurements, usually referred to as X, Y and Z. Nevertheless, the residual noises in these three measurements are correlated as the delayed residual noises from different lasers can appear in several of the measurements.

For statistical purposes, it can be easier to work on data with de-correlated noises. For TDI of first generation 1 , they are defined as :

A = Z -X √ 2 , (3.2.4) E = X -2Y + Z √ 6 , (3.2.5)
T = X + Y + Z √ 3 . (3.2.6)
Channel T is known as "null channel", because it is less sensitive to low-frequency GW [9; 10]. It can be used for other purposes than GW search, such as noise monitoring.

Many more TDI combinations are possible [11], depending on the assumptions we rely on. Even more, several versions of TDI are available, each of them responding to more and more realistic modeling of the instrument (armlength rigid or not, equal or not, large wavelength approximation). All the different versions are reported in [9]. [12] also proposes a generalized version of TDI based on principle component analysis that could bypass the armlengths estimation.

It is still possible to work on data with de-correlated noises for TDI later versions : [5] states the existence of de-correlated TDI as a mission requirement. In the following work, we often rely on de-correlated data as it is statistically easier to process (no crosscorrelation terms have to be taken into account when computing a likelihood for instance) in a joint analysis of TDI channels.

Disturbances & noise sources

Noise sources

Many events can impact the viability of the measurements ; they are classified under the "noise" or the "artifact" labels. Fig. 3.6 indicates the dominant noise contributions as a function of Fourier frequency range. They originate mostly from either physical origin or from the optical metric system.

Noises of physical origin mainly correspond to forces acting on the test-masses, be it magnetic forces, thermal fluctuations or others. They create test-mass position perturbations that do not correspond to GW signal. On Fig. 3.6, they mostly correspond to low-frequency noises (f < 3 mHz).

Noises originating from the Optical Metric System directly happen during the measurement process : among them we find the laser noise, the photon shot noise misleading the phase measurement, clock noise, spacecraft jitters, misalignment issues on the optical path, temperature effect on electrical components (that could also taint the data with correlated noises), but also the impact of switching on and off any device on a spacecraft. This last noise source was actually observed during the LPF mission [14]. On Fig. 3.6, they mostly correspond to higher frequencies noises (f > 3 mHz).

Taking into account all the different types of noise, we can give a "current best estimate" of instrumental noise performances. Keeping scientific objectives in mind, they should remain as much as possible below the requirements stated in the LISA Science Requirement Document [5]. The LPF mission, that was build as a demonstrator for LISA instrument, showed that the test masses feedback control system can efficiently reduce the level of the noise related to that origin [15; 16]. Fig. 3.7 shows this current best estimate compared to various source types that are expected : MBHBs, GBs and EMRIs. The computations related to the instrumental sensitivity are detailed in [17]. As a first approximation, most of these noises can be modeled as a zero-mean Gaussian colored noise in Fourier domain which estimated power spectral density (PSD) is given by the noise current best estimate. Yet, several other disturbances cannot be modelled as such. We gather them under the "artifacts" label.

Artifacts

In order to compute LISA sensitivity, the maximum interferometric noise level and test mass acceleration noise level are assumed to be known. Based on arm response simulations, they are propagated to compute the final instrumental sensitivity [17]. In reality, these noise sources could be greater than expected and other noise sources could deteriorate the performances. This is why they have to be dealt with before (as a pre-processing step) or during data analysis.

First and foremost, the laser noise has to be cancelled immediately as it represents the noise source with the highest amplitude among all. This is done trough the application of the TDI process that was presented in Section 3.2.2. More generally, the Initial Noise Reduction Pipeline (INReP) aims at suppressing Optical Metric System originated noises as much as possible (laser noise, clock noise, etc). This is the first step -and one of the most important -of data pre-processing.

After applying the TDI, many noise sources are left. Any noise that is not included in the description made in Fig. 3.7 is assimilated to a noise artifact ; an exhaustive listing of all of them can be found in [18]. Two types of artifacts are of interest here : Glitches : Their occurrences were observed during the LPF mission [15; 16]. A study was lead to determine their characteristics [14] in order to build a realistic simulator for the LISA mission. They have either electronic or mechanical origins (test-mass accelerations that are damped by the Drag-Free system)

Gaps : Data taking interruptions are likely to happen as it was observed during LPF mission. However, if they are too frequent, they will impact data analysis not only because a large amount of data will be missing, but also because several types of studies are directly conducted in Fourier domain. For LISA, data gaps could be scheduled (spacecraft maintenance), unscheduled (if the spacecraft goes into safety mode for instance), or we could even "gap" the data (as putting aside data sequences) that are too noisy to be used, if too many glitches are present for instance. The duty cycle (ratio of usable data measurement) could drop down to 75% [19].

Other types of artifacts exist, like the spectral lines (noises appearing at a specific frequency) or non-stationary noises. They are not investigated here. Since these non-gaussianities could affect the production of LISA science from the raw measurements, they have to be carefully described and integrated in the data analysis. In the following work, we will mainly focus on data gaps. However, a part of the work we present could presumably be adapted to deal with glitches (be it for detection or subtraction).

Chapter 4

Data analysis methods for gravitational wave signals

In this chapter are gathered the necessary elements to understand the work I did during my thesis. The object of study (galactic binaries) is presented and we define the data models (in terms of signal and noise) that will be used throughout the manuscript. We also introduce the Bayesian approach to LISA data analysis -this is the most commonly used approach in the community. Its functioning is demonstrated on a simple example based on a sampling code that I developed. The non-parametric approach to data analysis is then introduced as a complementary approach to Bayesian analysis. Most part of my work focused on GB sources from modeling them to dealing with them in presence of data gaps. As described in Section 3.1, we call a system of "galactic binary" two bodies with similar masses, rotating one around another, far from coalescence and thus emitting a quasi-periodic GW signal. For instance it can be made of black holes, neutron stars or white dwarfs. We illustrated in Section 2.2 that in the source frame such a system without eccentricity emits a sinusoidal signal. In reality, the phase has to take into account an additional term [2] :

Galactic Binaries : parametric representation

Φ(t) = ϕ 0 + 2πf 0 t + π ḟ0 t 2 , (4.1.1)
with Φ 0 the initial phase and f 0 the emission frequency that is twice the orbital frequency. The term π ḟ0 t 2 conveys that the system is not strictly periodic : its emission frequency can change over time due to e.g. tidal effects or the two stars getting closer. We want to deduce the GW signature that such a system will leave on the TDI channels. This is done in several steps taking into account both the source position and the instrument characteristics : ■ First, the signal has an amplitude h 0 and is decomposed into its two polarizations [1] :

h S + = h 0 (1 + cos 2 ι) cos [Φ(t)] , (4.1.2)
h S × = -2h 0 cos ι sin [Φ(t)] . (4.1.3)
■ Then, we apply the change of frame from the Source frame to the Solar System Barycenter (SSB) frame detailed in Fig. 4.1 (more details can be found in [1]). The frame change introduces the sky position of the source in terms of the ecliptic latitude λ and longitude β, as well as the system's polarization angle ψ. ■ The third step is to combine the GW signal in the SSB frame with the LISA orbit to compute the instrumental response : This can be done for an analytical simple orbit, as in [2]. With this step, one is able to compute the different interferometric measurements.

■ The final step consists in computing the signal in the TDI channels used in data analysis. An example of computation is detailed in [2]. In time domain, we observe that the "high frequency" signal, which is the one that we seek, admits a very low frequency envelop that is due to LISA constellation movement around the Sun. In Fourier domain, this Doppler effect causes the signal to spread over neighboring Fourier frequencies.

After applying all these steps, we obtain the signature of a GW signal emitted by a GB directly visible on the TDI channels. A typical example is represented in Fig. 4.2 : in time domain, the measured signal (also called waveform) locally appears as sinusoidal. However, when considering data on a larger timescale, we can see that the signal's amplitude is modulated by the instrument's orbits. In Fourier domain, the signal appears from afar as a single peak (which is expected for a monochromatic signal). However, when zooming on the peak, we observe that it has a width : the signal spreads over several frequencies because the instrument's motion generates a Doppler shift that depends on the position of the GB system.

Intrinsic Parameters Extrinsic Parameters

Frequency Therefore, we have a computable formula for GB expected TDI signal based on the 8 parameters gathered in Table 4.1. They are often split into two categories : intrinsic and extrinsic parameters, as the two categories play asymmetric roles in Bayesian parameter estimation.

f 0 Hz Amplitude h 0 Strain Frequency Derivative ḟ0 Hz -2 Inclination ι ∈ [0, π] Rad Ecliptic Latitude β ∈ [-π 2 , π 2 ] Rad Polarization ψ ∈ [ 3π 4 , 5π 4 ] Rad Ecliptic Longitude λ ∈ [-π, π] Rad Initial Phase ϕ 0 ∈ [ π 2 , 3π 2 ] Rad
These parameters can be related to quantities of interest in astrophysics and cosmology. The amplitude h 0 and the frequency derivative can be related to the system chirp mass M c and its distance to the Sun D L [2] :

M c = m 3/5 1 m 3/5 2 (m 1 + m 2 ) 1/5 , (4.1.4) ḟ0 = 48 5 × 2π GM c 2c 3 5/3 (2πf 0 ) 11/3 , (4.1.5)
h 0 = 4(GM c ) 5/3 c 4 D L f 0 π 2/3 , (4.1.6)
where m 1 , m 2 are the respective masses of the two stars, G is the gravitational constant and c is the light celerity 1 . The chirp mass M c and the distance to the Sun D L are quantities that characterize the position and composition of the observed GBs [3]. We will denote a set of parameters by :

θ = (f 0 , ḟ , β, λ θ int , h 0 , ι, ψ, ϕ 0 θ ext ) = (θ int , θ ext ) . (4.1.7)
At this point, we take the opportunity to distinguish between two types of parameters : extrinsic parameters θ ext and intrinsic parameters θ i nt. This distinction is linked to observations that were made in order to speed up parameter estimation [2] where they play asymmetric roles. See Section 4.3 and then App. A where this separation is explained in detail.

For any type of system emitting GW, we call waveform the function generating the LISA measurements corresponding to the physical parameters θ. For some systems like GBs, [1] provides a fast algorithm to produce these signals in the case where the corresponding system has a non-eccentric orbit.

As GB signals are quasi-stationary with low amplitudes, it is easier to spot them directly in Fourier domain because all the information is gathered around the emission frequency f 0 . This is why the search and identification of GB signal is usually conducted directly in this domain. It enables to localize the search over restricted frequency ranges (instead of considering all of the time measurements) which contributes greatly to speeding up detection and identification. From its observed signature, one of the stakes of data analysis is to trace back the source and its corresponding physical parameters ; Bayesian parameter estimation provides a framework that can address this inverse problem. . The total signal is plotted in blue. The galaxy is made of millions of GB systems, all emitting GWs in a restricted frequency band. Some signals are brighter than the others, which will facilitate their identification. Others will be assimilated to a continuous background noise, the "confusion noise". (Source : [4], challenge 2a) Reality is much more complex than the case of a single GB : scientists estimate that GB population in our Galaxy could reach about tens of million systems [3]. A realistic simulation of the signal produced by such a population is presented in Fig. 4.3 : all the systems are emitting in the same frequency range 1 -20 mHz, many are overlapping and mixing up making the identification really hard. Among all these systems, the scientific objective is to identify about 20, 000 -the leftovers will be assimilated to a Gaussian, non-stationary noise [3].

Realistic data : the Galaxy challenge

Data & Noise modeling 4.2.1 Data modeling

We adopt the following notations : the data will be observed on a set of C Michelsontype channels denoted by the index I. For instance, we can have I ∈ {X, Y, Z} (C = 3) or I ∈ {A, E} (C = 2) if one prefer to consider channels with de-correlated noises -see Section 3.2.2 for their definition. The observed data on channel I will be noted d I , and

the overall data d =      d I 1 . . . d I C     
will gather the data on all channels.

The measurements are done on an observation duration of T obs with time step ∆t. This represents N = T obs ∆t data points, taken at time t n = n∆t. The Fourier transform of the data is defined for for frequencies

f k = kδf , with δf = 1 T obs and k ∈ [-K, K], K = N 2 .
To keep notations compact, we will use the same symbols d I , d for functions or their Fourier transform, and use the argument to distinguish both :

d I [n], d[n] (resp. d I [k], d[k]
) for the time (resp. frequency) dependence. The discrete Fourier transform convention used here is the following :

d I [k] = N n=1 d I [n]e -2iknπ N . (4.2.1)
The inverse discrete Fourier transform is defined as :

d I [n] = 1 N K k=-K d I [k]e -2iknπ N . (4.2.2) d[k]
is the concatenation of the data in frequency domain :

d[k] =      d I 1 [k]
. . .

d I C [k]     
The measurements on channel I are the sum of the GW signal h I and the instrumental noise n I :

d I = h I + n I , (4.2.3) d = h + n . (4.2.4)
h I is the signal measured on the TDI channel I, and thus directly depends on the system parameters θ. We will directly highlight this dependence by writing h I (θ) when needed. In fact, h I (θ) designate the waveformi.e. the function that produces the time measurements on the TDI channel I from the physical parameters θ.

Noise modeling

In a first approximation, the noise n I is supposed to be Gaussian and colored in Fourier domain. Thus, its distribution is characterized by its mean (assumed to be 0-valued) and its PSD, defined in Fourier domain. The expected noise PSD for LISA is the one presented in Fig. 3.7. The PSD associated to the channel I is denoted by S I . It is assumed to be diagonal (noises at different frequencies are de-correlated). The cross-correlation PSD between TDI channels I and J will be denoted by S IJ . The corresponding PSD on d reads :

S =      S I 1 . . . S I 1 I C . . . . . . . . . S I C I 1 . . . S I C      . (4.2.5)
In case where the TDI channels are de-correlated in noise, this matrix is diagonali.e.

S I k is diagonal for any k and S I j I k = 0 for any j ̸ = k -because all the cross-correlation terms vanish.

The noise PSD on channel I is defined by :

S I [k] = E [n I [k] • n I [k] * ] . (4.2.6)
More details about the PSD can be found in Appendix D.4 and in [5]. Under the hypothesis that the noise is zero-mean, stationary, Gaussian colored in Fourier domain, the noise model writes :

n I [k] ∼ N 0, S I [k] 2 + i • N 0, S I [k] 2 , (4.2.7)
where N (µ, σ 2 ) denotes the Gaussian law of mean µ and variance σ 2 .

Bayesian estimators for the signal and the parameters 4.3.1 The Bayesian vision of data analysis

Bayesian inference is a framework often used to deal with data analysis based on statistical models. Suppose that the observed data is given by :

d = h * + n , (4.3.1)
with h * the "true" signal to estimate and n the noise with statistical description presented in Section 4.2.2. The Bayesian framework allows to estimate the probability of the underlying signal being h if the observation is d, and then set an estimator h of h * as the maximum a posteriori of the estimated probability distribution.

For any signal h, the probability density of observing the data d knowing that the underlying signal is h for sure, noted p(d|h), is given directly by the distribution of the residual dh = n, which exactly corresponds to the noise distribution. Assuming the noise follows the distribution that we presented in the Section 4.2.2, its probability density is the one of a Gaussian noise realization :

p(d|h) = 1 2π det(S) CN exp - 1 2 (d -h) † S -1 (d -h) ⟨d -h, d -h⟩ S , (4.3.2)
where † denotes the conjugate transpose of a vector and as mentioned above, N (resp. C) stands for the number of data points (resp. channels). This involves the inner product ⟨•, •⟩ S defined as :

⟨u I , v J ⟩ S IJ = k u I [k]S -1 IJ [k]v J [k] ⟨u, v⟩ S = I̸ =J ⟨u I , v J ⟩ S IJ + I ⟨u I , v I ⟩ S I (4.3.3)
Now, we would like to solve the opposite problem, i.e. to estimate the probability of presence p(h|d) of a signal h knowing the measurements d. This is possible using Bayes' theorem :

p(h|d) = p(d|h)p(h) p(d) , (4.3.4)
where p(h) contains the prior information on the signal h and p(d) is the model evidence associated to data d. p(d|h) is the likelihood of measuring d when the underlying signal is h. The probability of presence p(h|d) is called posterior distribution ; for simplicity and when there is no ambiguity about the data vector d, we will denote it L(h).

Usually, we work with the logarithm of the posterior distribution :

log p(h|d) = log L(h) = log [p(d|h)] + log [p(h)] -log [p(d)] ∼ - 1 2 ⟨d -h, d -h⟩ S + log [p(h)] , (4.3.5)
as we only keep the terms depending on the waveform h2 . The maximum a posteriori estimator of the signal is defined by maximizing the log-posterior :

h = Argmax h log L(h) .
θ = Argmax θ log L p (θ) , (4.3.7)
using the associated log-posterior :

log L p (θ) = - 1 2 ⟨d -h(θ), d -h(θ)⟩ S + log [p(θ)] . (4.3.8)
Two remarks are in order. First, we define estimators Eq.(4.3.6) and Eq.(4.3.7) based on the hypothesis that, when using a waveform representing exactly the real signal, then the maximum is unique and unbiased. But even if there is no evidence that this hypothesis is verified, this approach is successful [6] in presence of non-artifact instrumental noise.

Second, we would like to highlight that Eq.(4.3.5) and Eq.(4.3.7) are two distinct estimators. The first one relies on the signal morphology (and the prior is set directly on the signal : p(h)) whereas the second one relies on the parametrized waveform (with a prior set on the parameters : p(θ)). The two types of prior can be fundamentally different ; most of the data analysis done so far has been using the estimator Eq.(4.3.7), whereas our work can be recast as a variante of estimator Eq.(4.3.5) as we will illustrate in Section 5.1.3 and Section 7.3.1. with regard to two of its variables, marginalizing the others. From up to down (and from left to right) : log 10 h 0 , f 0 (mHz), log 10 ḟ0 , sin(β),λ,cos(ι),ψ,ϕ 0 . Estimating the true values of the parameters (blue line) can be hard because the log-posterior is multi-modal and has not the same sensitivity with regard to all parameters.

Sampling the posterior distribution

Even if we defined the estimator as the argument maximizing the posterior, looking for the extremum of the log-posterior is not enough. It is more interesting to estimate the behavior of the log-posterior around the maximum since it can provide a way to compute error bars on the estimator.

Markov Chain Monte Carlo algorithms

This is usually done using an algorithm of type Markov Chain Monte Carlo (MCMC), such as [8; 9; 10]

θ t+1 = θ ′ .
-If u > α : reject the new sample :

θ t+1 = θ t .
This basic algorithm can be declined on much more complex processes, involving for instance interactions between the different chains, parallel-tempered chains or many more [10].

The final result is obtained by looking at the distribution of the Markov chains over the parameter space : if the chains are long enough, they converge to a stationary distribution p(θ|d) with no residual effect from the starting point.

Contribution & Typical results

I developed my own parameter estimation code based on the principles described in Section 4.3 and App. A. The sampler we use is PTEMCEE [9] with 4, 000 walkers and chains with 3, 000 steps. The number of temperatures was set to 5. The computing time needed to obtain a sampling with fair quality remains high, even if it can be decreased by using multiple cores or CPU. It takes about 30 minutes on 96 cores on our calculator to estimate the parameters of a single GB.

In order to sample properly the log-posterior, we did not directly work on the physical parameters described in Section 4.1. Instead, we worked with transformations of the parameters detailed in Table 4.2 on which we applied very simple priors p(θ).

This choice of prior is justify by the following approach : ■ Uniform prior over λ, sin(β) : this distribution results in a sampling that is uniform over the sphere, i.e. over the sky location 3 .

■ Uniform prior over ψ, cos(ι) : this distribution results in a sampling that is uniform over the sphere, i.e. over all the possible local orientations.

■ Uniform prior over ϕ 0 : we do not have any information about the initial phase. ■ Uniform prior over mf 0 : without any more physical information over the GB population, this prior enables to search over all the possible main frequencies uniformly.

■ Uniform prior over log 10 (A), log 10 ( ḟ ) : prior used for a fast first search to have an estimate of the order of magnitude of these parameters. While it might not be mathematically correct, it makes it possible to probe several order of magnitude at once -the first guess hence obtained can be used for a more detailed search with uniform priors directly over the parameters (A, ḟ ). We stopped at the rough

Parameter Prior

Frequency We can give a visual representation of this posterior distribution such as Fig. 4.4. As we are sampling a function of parameters that are not necessarily independent, we have to consider the joint impact of the parameters. The corner plot [7] Fig. 4.4 shows the behavior of the log-posterior L p (θ) with regard to two of its variables, marginalizing the others. Several observations can be done : first, the log-posterior is a multi-modal function, especially with regard to the variables ψ and ϕ 0 . The estimation of some parameters is particularly challenging, such as the frequency derivative ḟ0 . Some parame- ters can be highly correlated, for instance the amplitude h 0 and the inclination ι, or the frequency and its derivative. Finally, other parameters such as f 0 , β and λ seem to be independent from the others, and are the easiest to estimate.

mf 0 = 10 -3 f 0 ∼ U([f min , f max ]) Frequency Derivative log 10 ( ḟ0 ) ∼ U([-19, -15]) Ecliptic Latitude sin(β) ∼ U([-1, 1]) Ecliptic Longitude λ ∼ U([-π, π]) Amplitude log 10 (h 0 ) ∼ U([-24, -20]) Inclination cos(ι) ∼ U([-1, 1]) Polarization ψ ∼ U([0., 2π]) Initial Phase ϕ 0 ∼ U([0., 2π])

Advantages & limits of the Bayesian approach

This approach has many advantages : for one, it gives an estimation of the physical parameters of the system emitting the GW signal, as well as the corresponding waveform. Moreover, using the results of the sampling, it is possible to assess the reliability of the results, for instance by giving error bars on the parameter estimates. Finally, this framework can be adapted to the case where more than one signals are present (simply by changing the model and the parameter space) or if the total number of sources is unknown [6]. It can be straightfowardly extented to handle different types of sources as long as one has a parametric representation of the sources. More details about the sampling algorithms developed for LISA can be found in [11].

Nevertheless, the log-posterior function is very ill-conditioned. Proof is that several parameters are highly correlated and the multi-modality of the log-posterior with regard to other parameters. This results in a high computational cost, even if it can be reduced by using multiple cores or CPUs. MCMC are subject to the curse of dimensionality : looking jointly for the parameters of several sources corresponds to such a wide space that it becomes arduous to probe. Particularly, not only are all the parameters potentially correlated, but also the regions corresponding to highly probable parameters become more and more distant one from another. Of course, the same principle applies when searching for an unknown number of sources as the model dimensionality has also to be assessed.

Moreover, this technique is extremely model dependent : the results are good as long as the parametric model somehow fits reality. But as we announced in Section 4.1, waveform parametric models are developed under simplified assumptions : they might not fit exactly the "true" signals. For instance for GBs, waveforms taking into account different kinds of effect like tidal [12] or magnetic [13] effects are under development ; these effects create waveforms with three strongly correlated peaks instead of a single one resulting in a real change in the harmonic structure of the sought signal. Furthermore, some sources might not have a parametric model -we must always keep in mind that unexpected signals could be seen by LISA.

Finally, we have to highlight that the sampling is highly impacted by the parameters we work with and the related priors that we consider. The one we used to compute the results shown in Fig. 4.4 are gathered in Table 4.2 : they did make the sampling more efficient, but working with the log of a parameter does not give the same final estimation precision than working with the parameter itself.

Artefacts : data gaps, glitches

Some recent studies [14; 15] showed that missing data can have a huge impact on Bayesian data analysis. In [15], the authors detail the impact of a gap on a MBHB parameter estimation. As the information related to these systems is localized around the merging time, the impact of a data gap depends on its position and its duration. Complementarily in [14] the authors investigate the impact of missing data on GB parameter estimation. In this case, classical parameter estimation is realized in Fourier domain : the occurrence of data gaps in time domain results in a deformation of the noise PSD that has potential repercussion on the detectability of GBs and their identification. This phenomenon is detailed in Chapter 6. A generalized sampling process is proposed where the missing data is considered as a parameter and thus estimated along the physical parameters.

Glitches also have an impact on the noise distribution. But since they are transientlike perturbations (similar to a Dirac) on the time measurements, they will have a lesser effect than the gaps on the Fourier transform because their Fourier transform is spread over the whole frequency spectrum ; however their occurrence rate could be mis-interpreted as a signal : Since we did not investigate this issue, we will not develop any further on this subject.

Truth is that there is a lot of work to be done concerning artefacts. LISA Data Challenges are encouraging the community to start taking them into account through Data Challenge 2b [4].

Non-parametric analysis : a complementary approach

In data analysis, any type of GW source calls for three courses of action : Detection : first step towards analysis, we have to be able do recognize an event happening and to categorize it as a known -or unknown -source type.

Subtraction : Cleanly subtracting the signal related to one source type can help recovering the signals related to other sources types, but also to estimate the noise PSD.

Identification : If the source waveform is known, it might be possible to estimate the physical parameters of the source from the measurements.

We would like to stress out the importance of the two first points, as the third one has an obvious outcome in the scientific community. Indeed, the detection step has its importance as it must be able to "see" events that may have not been foreseen by the scientific community. LISA has a high potential for new discoveries, and this should be taken into account from the start as to not automatically assign a detection with a known source type. The key of a detection process is to be able to confront different models and to propose several assumptions -which "unknown" and "glitch" are also part of -to explain a given GW event.

Secondly, we have to consider our ability to efficiently subtract any type of signal from the data. Two major applications can be thought of : noise PSD estimation and parameter estimation. At the beginning of data taking, the noise distribution will still be unknown and will have to be estimated. A proper way to do that would be to subtract any type of event from the data -be it GW event or glitches or other systematicsand then estimate the noise distribution. Concerning parameter estimation, as only the sum of all GW signals are observed, signals of different source types could potentially contaminate each other, making the identification process harder or even biased. A basic solution to this issue would be to subtract the signals coming from all types of sources other than the one we are working on to facilitate the identification process.

At the same time, one has to think of dealing with the artefacts : should it be done as a pre-processing ? During the detection step ? Considering a parametrizable glitch, should its parameters be estimated at the same time as true GW source parameters ? Or should it be identified and dealt with before starting the analysis ? The same questions can be asked for gapped data : should gaps be dealt with before performing the analysis, or during the analysis ?

The natural answer would be to do a joint analysis of artifacts and data in order to simultaneously identify all of the components of the measured signal. In practice, that will likely be too complex to be performed properly. The question of performing a pre-processing step before the analysis is done is worth asking as it could prove to be a good approach to analyze such complex data.

In this context, we want to explore the relevance impact of non-parametric analysis for LISA physics. Non-parametric modeling is a framework relying on a signal representation that is not based on the physical parameters θ, but on what we expect to observei.e. the signal's morphology. Such methods were already introduced in the GW community [16; 17; 18] -essentially within LVC -but their use remained limited.

[18] states clearly the stakes of non-parametric search : "Transient GW searches can be divided into two main families of approaches : modelled and unmodelled searches, based on matched filtering techniques and time-frequency excess power identification respectively. The former, mostly applied in the context of compact binary searches, relies on the precise knowledge of the expected GW phase evolution. [...] The goal of non-parametric modeling is to obtain the best possible representation of the sought-for signal without using the physical parameters, while still ensuring the low-dimensionality of the representation. This low-dimensionality principle is crucial as it proves to be helpful when dealing not only with instrumental noise but also with the fact that many different types of sources are present in the data, each of them presenting a specific morphology. A non-parametric framework can lead to a wide range of applications, such as detection (Chapters 5, 7), source separation (Chapter 7), signal estimation (Chapters 5, 7) and a way to deal with data gaps (Chapters 6) -to cite the one closest to the above mentioned issues. Moreover, as it relies on signal's compressibility, it opens a world of possible with regard to fast waveform production and waveform approximations (Chapter 5,7).

Contributions

My PhD focuses on exploring non-parametric data analysis methods for LISA data analysis. By "non-parametric ", we mean that the model should not rely explicitly on the GB system physical parameters. Instead, we build a signal estimator based on finding a representation that is "adapted" to its expected shape. All the algorithms were developed while working on GB signals, but in principle could be generalized to other types of signals -being GW signals or artifacts such as glitches -to perform efficient data pre-processing.

From a very simple non-parametric model, we prove that we are able to detect GB signals and deal efficiently with data gaps. We demonstrate the performances of the proposed algorithms over an extensive benchmark. This lead to two publications [19; 20].

Then we present an innovative way to build non-parametric non-linear models for GW signals with promising results in terms of representation quality, source separation and fast parameters estimation. We explore several applications possible for these types of models and show how they could be used by the LISA community. A publication corresponding to this work is currently in preparation.

Chapter 5

Sparse signal modeling

The Sparse Way

Robust source detection and signal subtraction are the connecting threads of my work. In a context where identification relies on waveforms designed under simplifying assumptions (absence of eccentricity, etc) that still potentially require noticeable computing cost1 , it is interesting to provide fast accurate detection algorithms that do not rely too much on parametric models (that could be biased since they are simplified). Similarly, quickly extracting all signals coming from a specific source type without identifying them presents an interest : for instance to estimate the noise PSD ; or to study another type of source which signal could be impacted by the presence of the aforementioned source type.

In the field of signal processing, the sparsity framework is particularly well adapted to address these questions. Sparse representations are already used in the LVC pipeline Coherent Waveburst [1] to improve detection. In this pipeline, the Wilson-Daubechies-Meyer transform is applied to the measurements. This transform represents the data in a redundant dictionary ; sparsity is used to focus the signal's power over the fewest coefficients possible to improve its detectability. Then this sparse representation is compared (using a graph) to the sparse representation of a bank of templates in order to identify the physical parameters of the source [2]. In [3], the authors directly try to recover the polarizations h + , h × of GW signals.

Sparsity was proved to be an efficient tool that is well adapted to GW context. We now wish to introduce this framework to the LISA community as it is flexible enough to deal with many problems faced by data analysis today. Among these, we can mention source detection, source separation and data gaps.

Sparsity Framework

Sparsity is highly linked to the concept of signal compressibility. Let us consider a dictionary T = t j which atoms t j are typically elementary templates representing adequately the signal. Then a time signal x ∈ R N can be expressed with in the dictionary T :

x = Tα = N T j=1 α j t j . (5.1.1)
We say that x is "encoded" by α because it can fully be recovered from the coefficients α j . Compression takes place when only few α j carry maximal information, i.e. few α j have a large amplitude, and the others have small amplitude in comparison. Two distinct configurations can be thought of.

Exact sparsity

Exact sparsity (or strong sparsity) is the idea that the signal x can be encoded exactly by L elements of dictionary T. Let Λ L be the set of the L indices required to encode x ; then we have :

x = j∈Λ L α j t j , (5.1.2)
since α j = 0 for any j / ∈ Λ L . The fewer elements are needed, the more the signal is compressible and the more the signal is sparse in T. In this context, we say that x is sparse in T if it is encoded only by few non-zero coefficients, i.e. if α = {α j } j has many null components. In many situations x is not strictly sparse, i.e. it mostly has non-zero components. However if the dictionary T is well chosen, then the α coefficients coding x will have a behavior characterizing what we call "weak sparsity" : a small proportion of coefficients will have a large amplitude, but this amplitude will fastly decrease to zero. A typical example of this phenomenon that we can mention is the power law-like decrease of Fourier coefficients associated to a smooth periodic signal in Fourier domain. We observe similar behaviors when applying bi-orthogonal wavelet transforms to images [4]. Fig. 5.1 gives an example of this phenomenon : the bi-orthogonal wavelet transform provides a dictionary that is well-adapted to representing pictures as the coefficients' amplitudes decrease relatively fast.

Weak sparsity

In this case, a natural move is to try to approximate x based on its L largest coefficients in the dictionary T :

x

≃ j∈Λ L α j t j , (5.1.3)
with Λ L the set containing the indices of the L largest coefficients α j . The bigger the considered set Λ L , the smaller the approximation error :

lim L→N T x - j∈Λ L α j t j 2 2 = 0 . (5.1.4)
This is the principle of signal compression : approximating a signal x by a L-sparse signal in a dictionary Ti.e. a signal with only L non-zero components in the dictionary T.

To go back to the example of Fig. 5.1, we could try to approach the cat picture by few elements of the bi-orthogonal wavelet transform : this is exactly the principle of the JPEG 2000 image compression format [5]. Such an approach results in an approximation error that also decreases as a power law of the number of non-zero coefficients used to approximate the signal [4].

More generally, we seek to give the best representation possible of x that will use the fewest components of T :

Argmin α ∥α∥ 0 such that ∥x -Tα∥ 2 ≤ σ , (L 0 )
where ∥α∥ 0 (known as "norm 0") counts the number of non-zero components of the vector α, ,and σ is the approximation error.

What is the interpretation of (L 0 ) ? (L 0 ) states that we are looking for an approximation of x in T such that the approximation error does not exceed σ. Moreover, this approximation Tα must rely on as few non-zero coefficients ∥α∥ 0 as possible. It is "the best" approximation of x based on "the fewest" atoms of T.

It has been shown [6] that Eq.(L 0 ) can be relaxed as :

Argmin α ∥α∥ p such that ∥x -Tα∥ 2 ≤ σ , (L p ) with p ≤ 1 and ∥α∥ p = j |α j | p 1/p
. In this case and under certain conditions, (L 0 ) and (L p ) have the same solution. This relaxation is interesting because contrary to (L 0 ), the problem (L p ) is continuous and can even become convex if we set p = 1. This is a noticeable result, even more since both continuity and convexity contribute to facilitate the search for the minimum. For more details, we refer the interested reader to [7; 4; 8; 9].

Sparse modelling for Galactic Binaries

The sparsity framework suits particularly well GBs. As mentioned in Section 4.1, GBs emit smooth, stationary, nearly sinusoidal signals. Based on this observation, the best (To increase the number of coefficients used to approximate the signal, we increased the input signal's amplitude while keeping the threshold at constant value). Middle : The number of components used for the approximation will mainly impact the estimation of the signal time envelop. The approximation gets better with increasing number of atoms considered for the approximation. Right : Zooming over a short time period, both estimates (with 34 and 80 components) are in phase with the real signal they try to approximate. solution is to proceed to their study directly in Fourier domain where all the information is gathered mostly on a few Fourier atoms. Hence, we can say that they are naturally weakly sparse in the Fourier atoms dictionary. This statement can be illustrated by a very simple example. Fig. 5.2 shows the approximation of a GB waveform by its L largest coefficients and the impact of this approximation on the time-representation of the waveform. Only considering the few Fourier atoms with the largest amplitude (left picture), the recovery in time domain already has a perfect phase match -even if the amplitude is poorly recovered. Now, doubling the number of coefficients improves considerably the amplitude recovery while keeping the phase match.

More generally, Fig. 5.3 shows how the number of Fourier coefficients (and thus Fourier atoms) used to approximate the signal impacts the quality of approximation Q dB (the quality factor Q dB will be defined in Section 5.5). The evolution of the quality as a function of the number of coefficients used for the approximation evolves in two stages. At first, the progression slope is steep : each added coefficient contributes greatly to the improvement of the signal quality. But after a certain number of coefficients (the limit is around 40 coefficients in this example), the slope breaks and the improvement brought by the addition of a coefficient becomes much weaker. Nonetheless, in presence of noise the recovery plot is degraded and the task of estimating the signal despite the noise becomes a real challenge.

Practical Resolution

Two elements must be accounted for : first, our GB signal x in only weakly sparse in T. Second, it is tainted with complex instrumental noise. These two elements are enough to understand that the recovery scheme Eq.(L 0 ) will not be perfect. However, using the sparsity denoising principle [4], we can build a good approximation of the signal underlying the noisy data. This denoising principle is founded on a thresholding process applied to the atoms coefficients α j . There are many possibilities ; the one we selected is obtained solving the LASSO convex problem [8] :

Argmin α γ ∥α∥ 1 + 1 2σ 2 ∥x -Tα∥ 2 2 , (LASSO)
where σ stands for the noise standard deviation in case where the noise is Gaussian.

What does (LASSO) do ? It can be interpreted as a tradeoff between "finding the sparsest solution possible", embodied by the ∥•∥ 1 term on the left, and "finding the solution that fits the measurements the best", embodied by the data fitting term on the right, that evaluates the square distance between the measurement x and its representation in the dictionary Tα. This tradeoff is balanced through the choice of the regularizing parameter γ.

(LASSO) has been used in many contexts because of its efficiency on denoising problems (i.e. the problem consisting in recovering a signal from a noisy measurement) [3; 4; 8]. Going back to GBs and the realistic case presented in Section 4.1.2, we understand that only a small part of the GB signal will exceed noise level. Moreover, many peaks (corresponding to GB signals) are likely to rise above noise level. Ultimately what is proposed here with (LASSO) is to approximate the sum of all GB signals at once using the sparsest representation possible in Fourier dictionary. We investigate the viability of this strategy, first by using a problem directly similar to (LASSO), and then by introducing a notion of "information cluster" [7; 4] that will help to improve the overall signal estimate.

In the case of GBs, if we work with the dictionary of Fourier atoms, then α exactly coincides with the signal's Fourier coefficients. This means that T can be assimilated to the Fourier transform operator, which is an invertible isometry -with T -1 the inverse Fourier transform. In this context, Parseval identify yields :

|| x Time data -Tα|| 2 2 = || T -1 x
Frequency data -α||2 2 .

(5.1.5)

Thus we can express (LASSO) directly in Fourier domain :

Argmin α γ ∥α∥ 1 + 1 2σ 2 ∥x[k] -α∥ 2 2 , (REF-L 1 )
where x[k] is the Fourier transform of the time signal x[n], as introduced in Section 4.2.

We will work directly with (REF-L 1 ) as it admits the same solution as (LASSO) but with the great advantage of not requiring to apply any transform. This makes the solution computation faster.

Notations

Using formulation (REF-L 1 ), we directly work in Fourier domain. We introduced in Section 4.2 the corresponding data d

[k]. For any k ∈ [-K, K], the k-th coefficient of data d[k]
is the Fourier coefficient associated to the Fourier atom :

t k = e -2πikn N (5.2.1)
We will also refer to t k as the k-th frequency (since it corresponds to frequency f k = k • δf ). We consider the de-correlated TDI channels (see Section 3.2.2) I ∈ {A, E}. We do not consider the channel T as it has by design a low sensitivity to low-frequency GW.

We will denote by v the signal variable expressing the optimisation problem (REF-L 1 ).

We introduce the vector of relative signal amplitude for every Fourier atom on channel I with regard to the noise level :

A S I (v I ) = {A S I (v I )[k]} k = v * I [k] • v I [k] S I [k] k = |v I [k]| S I [k] k . 2 (5.2.2)
This quantity is interesting because it directly shows if there is a power excess of the data with regard to the noise level. Moreover,

k A S I (v I )[k]
can be interpreted as the norm 1 of v I re-weighted (or whitened) by the noise standard deviation :

∥A S I (v I )∥ 1 = k |A S I (v I )[k]| = S -1/2 I v I 1 = k |v I [k]| S 1/2 I [k]
.

(5.2.3)

Note here that we assumed that the noise power spectrum is colored in Fourier domain.

Concretely, this means that the noise associated to two different atoms will not have the same standard deviation : this has to be taken into account to properly recover the signal. It is all the more true as the noise PSD can vary over several orders of magnitude (See Fig. 3.

7).

As we decided to consider the dictionary made of Fourier atoms, we impose a "sparse distribution of signal amplitudes" in this dictionary as a prior,i.e. we want the coefficients of the identified signal to be spread over the fewest Fourier atoms possible. The "sparse" variable will therefore be the signal amplitude.

The advantage of Eq. (5.2.2) is that it can easily be generalized to the multi-channel case. As it should be done for uncorrelated channels, the total quadratic amplitude associated to the atom k will be given by the sum of the quadratic amplitudes of individual channels :

A S (v) = I v * I [k] • v I [k] S I [k] k =    I A 2 S I (v I )[k]    k .
(5.2.4)

In this case, the sum over atoms still enforces the sparsity with a 1-norm, but this time over the overall amplitude. This means that the coefficients related to the k th Fourier atom are either zero for all channels, or non-zero for all channels (see the reasons for introducing this operator in Section 5.3.2 ).

Last, we will use the Hadamard product ⊙ such that for two vectors u, v of identical size :

u ⊙ v [k] = u[k]v[k] ∀k ∈ [-K, K] .
(5.2.5)

A separable sparse model : unstructured sparsity

The work presented in the following sections was published in Physical Review D [10]. We consider the following data model :

d I = h I + n I , (5.3.1) d = h + n . (5.3.2)
h I is the sought signal and is assumed to be sparse in Fourier basis. n I is the Gaussian noise contaminating the data.

Single channel model

Let us now rewrite the general reference problem (REF-L 1 ) taking into account the specificities mentioned in Section 5.2. Those include the colored noise distribution over Fourier atoms, and the notion of "amplitude" that we expect to sparsify.

For a single channel I, the data d I is expressed by its coordinates in Fourier dictionary. This means that we directly look for its sparsest approximation h I in the Fourier dictionary :

h I = Argmin v γ ∥A S I (v I )∥ 1 γ k A S I (v I )[k] + 1 2 ⟨d I -v I , d I -v I ⟩ S I ∥d I -v I ∥ 2 2,S I . (5.3.3)
We recognize a problem with a structure similar to the one described in Eq.(REF-L 1 ). The data fitting term is weighted by the expected noise PSD, and the sparsity term is applied as mentioned before directly to the estimated signal amplitude. Nevertheless, the reweighted L1 formulation [11] advises to consider a model a bit more general to limit the bias induced by the coefficient γ on the recovered solution :

h I = Argmin v k γ[k]A S I (v I )[k] + 1 2 ⟨d I -v I , d I -v I ⟩ S I , (5.3.4)
where the regularizing parameter γ of Eq. (5.3.3) is replaced by an atom-dependent positive real regularizing parameter {γ[k]} -K≤k≤K . This choice will be justified when considering the impact of the parameter γ on the corresponding analytical solution ;

this will also help setting the value of this parameter. Following the notations introduced in Section 5.2, we can rewrite the single channel problem as :

h I = Argmin v ∥γ ⊙ A S I (v I )∥ 1 + 1 2 ⟨d I -v I , d I -v I ⟩ S I , (L 1 , I)
Fourier atoms are orthogonal ; therefore the problem is entirely separable 3 . The estimator h I can now be computed analytically atom by atom. The solution is given by the so-called soft-thresholding operator :

h I [k] =    A S I (d I )[k] -γ[k] A S I (d I )[k] d I [k] if A S I (d I )[k] ≥ γ[k] , 0 otherwise . (5.3.5)
The regularizing parameter γ acts as a threshold since the measurement d I [k] corresponding to atom k will be discarded if its amplitude is smaller than γ[k], and will be replaced by its excess to γ[k] otherwise. We will label the atom k as active if the signal amplitude exceeds the threshold γ, and inactive otherwise.

The soft thresholding process is illustrated on Fig. 5.4 on a very simple sinusoidal signal : soft thresholding consists in zeroing all sub-threshold amplitudes in Fourier domain, whereas the amplitudes above threshold are decreased by a factor γ. This explains why the sparse estimate in time domain does not have the right amplitude. We remark that the result would have been the same if we had considered Eq.(5.3.3) instead of Eq.(L 1 , I), only that γ[k] would have been a constant.

In the following sections, we will lay down a joint analysis of the A and E channels, and further elaborate on the choice of the threshold γ. However, the structure of the considered problem will always remain the same throughout the chapter. Only the sparsity-enforcing term and the minimizing algorithm will be sophisticated in order to adapt them to data processing situations of increasing modeling complexity.

Combining different information channels with joint sparsity

While the solution stemming from (L 1 , I) provides important insights into the algorithmic nature of the problem, it does not benefit from the redundancy of the physical content over multiple channels, which have been so far treated separately. Indeed the information contained in TDI channels are highly correlated (the imprint of the same GW is encoded in all channels), whereas the noise is totally uncorrelated (see Section 4.2.2). This leads to the following two observations : ■ A signal emitted by a physical source at a given frequency f 0 will be measured in all channels, although with different amplitudes. For instance, this signal could be visible in one channel, but more difficult to detect in the other one.

■ The noise impacting the coefficient of the k th Fourier atom could be exceptionally high on a channel (i.e. be an outlier), but low on the others. An independent analysis of the two channels could output a false positive signal on one (and only one) channel for this noise realisation. Jointly processing all the channels addresses both problems and is the natural extension of the preceding development.

The very structure of the problem of constructing the sparse estimator h I is preserved by making use of the compact notation h gathering all information channels (see Section 4.2.1 ). Using the definition of the joint amplitude Eq. (5.2.4) over multiple data channels, a joint sparse representation over multiple channels is achieved through :

h = Argmin v ∥γ ⊙ A S (v)∥ 1 + 1 2 ⟨d -v, d -v⟩ S , (Joint L 1 )
where γ = {γ[k]} -K≤k≤K still denotes a positive real threshold that is shared by all channels. As in Eq. (5.3.5), the analytical solution is given by :

h[k] =    A S (d)[k] -γ[k] A S (d)[k] d[k] if A S (d)[k] ≥ γ[k] , 0 otherwise . (5.3.6)
The problem is still separable and can still be solved atom by atom, even if we now jointly examine all the channels. The role played by γ is still clearly that of a threshold, but now it is the combined amplitude A(v)[k] of the channels that should excess the threshold to tag a Fourier atom as active.

Similar approaches for multi-channels data processing with sparse signal representations can be found in [12].

Fixing the threshold as hypothesis testing

Whether based on (L 1 , I) (see Section 5.3.1) or (Joint L 1 ) (see Section 5.3.2) sparsityenforcing terms, building estimators for GB events eventually led to a thresholding operation over Fourier atoms' coefficients. Going back to the simple example shown in Fig. 5.4, we understand that we want a threshold high enough so that all Fourier coefficients corresponding to noise will be put to zero, while keeping it as low as possible so that the amplitude is not too biased.

It turns out that fixing the threshold γ can be recast as an hypothesis testing problem, which provides important insights into the construction of the estimator itself. This threshold choice can be performed independently for each atom k.

We classically define the two hypotheses H 0 and H 1 :

H 0 : there is no GW signal at the k-th Fourier atom ;

H 1 : there is a GW signal at the k-th Fourier atom.

We assumed (Section 4.2.2) that the real and imaginary parts of n I [k] obey independent standard normal distributions. Consequently, under H 0 (only noise) A S I (d I ) 2 [k] admits a chi-square distribution with 2 degrees of freedom. Now if we consider C uncorrelated channels for joint resolution, then A S (d) 2 [k] admits a chi-square distribution with 2 • C degrees of freedom.

Elaborating on Eq. (5.3.5) and Eq. (5.3.6), this provides a criterion to set the threshold value based on a p-value test. Adopting a priori a rejection rate ρ and defining the real number x 0 by the value of the cumulative distribution function :

P(χ 2 2C ≥ x 0 ) = ρ , (5.3.7) the hypothesis H 0 (resp. H 1 ) is adopted for frequency k if A S (d) 2 [k] ≤ x 0 = γ 2 [k] (resp. A S (d) 2 [k] > x 0 = γ 2 [k]
). The same test can be adopted for a single channel I by setting

C = 1.
If the noise distribution is not assumed to be standard normal, then A S (d) 2 [k] will not obey a χ 2 distribution. However, the reasoning above can be adapted to the actual distribution of A S (d) 2 [k], and the principle of fixing the threshold as hypothesis testing remains.

Reweighted minimization

As discussed in Section 5.3.1, the estimator h I [k] derived from Eq. (5.3.5) or Eq. (5.3.6) is a rescaling of the noisy signal d I [k] by a multiplicative factor in [0, 1], which makes it intrinsically biased. This is understandable considering the case where the data is exactly the signal : d = h. Then, the estimator writes :

h[k] =    A S (h)[k] -γ[k] A S (h)[k] h[k] if A S (h)[k] ≥ γ[k] , 0 otherwise .
(5.3.8)

For a threshold γ[k] that is non-zero, the estimated amplitude of the signal is biased by a factor

A S (h)[k] -γ[k] A S (h)[k] ∈ [0, 1]
. This bias becomes even more noticeable as γ increases.

This was already shown in Fig. 5.4.

In the general case of noisy data d = h + n, we are faced with a dilemma : ■ Either we set a low γ : the bias is low, but there is a high probability of detecting an atom as active whereas it only corresponds to noise (false positive).

■ Or we set a high γ : in this case, the probability of false positive is low, but the bias on the solution is important. In the best case scenario, we want to have both a low false-positive rate and a low bias on the estimated signal amplitude.

We present here an iterative process, named reweighting and described in [11], to correct this bias. It consists in iteratively alternating at each step m ∈ N between a resolution phase, yielding a signal estimator h m I [k], and a threshold determination phase, producing a regularizing parameter γ m [k]. The idea is, after a first selective step of detection, to set a low threshold for active atoms and keep a high threshold for the others. Therefore the bias correction only acts on active atoms. The convergence of this algorithm has been established in [9] and in practice only few (about 3-4) iterations are needed.

As in Section 5.3.3, we set a rejection rate ρ and define a (square of) threshold x 0 through Eq. (5.3.7). We initialize the iterative procedure with the following threshold :

γ 0 [k] = √ x 0 ∀k ∈ [-K, K] . (5.3.9)
We target a residual ϵ γ > 0 which will act as a stopping criterion : the algorithm will stop at the first step m such that max

k |γ m+1 [k]-γ m [k]| < ϵ γ .
Each iteration of the reweighting procedure performs as follows :

Phase 1 : Using γ m [k], compute the estimator h m at step m from Eq. (5.3.6).

Phase 2 : Using h m , evaluate the threshold γ m+1 [k] as described in [11] :

γ m+1 [k] = (γ 0 [k]) 2 κA S ( h m )[k] + γ 0 [k] , (5.3.10)
where κ is a positive real parameter that is set either to amplify (κ > 1) or to reduce (κ < 1) the basic reweighting (κ = 1). Using κ ∈ {1, . . . , 10} yields quantitatively similar results, and we have set κ = 3 throughout this study.

Note that the reweighting procedure impacts only active atoms (i.e. the atoms k satisfying A( h m )[k] > 0). The higher the norm of the solution, the greater is the correction and the lower becomes the bias, as shown in Fig. 5.5.

Adding physical information to the model through structured sparsity

Whether dealing with one or several channels, none of the methods we introduced so far make use of any physical input on the shape of the GW signal beyond its quasimonochromatic nature. Since these methods use no information about the structure of the GB waveform, they will be further dubbed unstructured sparsity-based methods. To provide more insight into these methods, we will more specifically consider the following two phenomena : ■ A Fourier atom for which no signal should be detected has been tagged as active.

We will refer to it as false positive (FP). This situation is illustrated in Fig. 5.6. ■ A Fourier atom for which a signal should be detected has not been tagged as active. We will refer to it as false negative (FN). This situation is illustrated in Fig. 5.6. In particular, estimating the coefficient attributed to a Fourier atom independently from the values taken by the others like in Eq. (5.3.5) or Eq. (5.3.6) leads to a high rate of both FPs and FNs.

If we consider the problem of detecting GB systems to create a catalogue of sources, then its viability will be limited if there is a high number of detected FPs. To be able to detect a low amplitude peak, the chosen initial threshold should not be too high. However, the lower the threshold, the higher the probability to get a FP signal. Therefore, an efficient and robust method to separate genuine signals from FPs is much needed.

This separation can be achieved by implementing the simple yet important remark : a GB signal is not exactly monochromatic. Instead, because of the constellation movements, the measured signal has a frequency that vary with time due to a Doppler effect. In practice, a measured GB signal will spread over few frequencies neighboring the main emission frequency f 0 . The signal power will be distributed over the corresponding Fourier atoms, as shown in Fig. 5.6 -upper row. The presence of a signal induces a significant correlation of the coefficients corresponding to atoms that are close in frequency (see upper figure). On the contrary, a FP is often the manifestation of a rare powerful noise realization impacting an isolated Fourier atom (see Fig. 5.6, lower figure) on a single channel. Considering the large number of records (the original time series contain tens of millions of data points), the presence of FPs is naturally expected.

Note that there are also simple mechanisms providing FNs. For example, choosing a atom-independent threshold leads to the rejection of all atoms with amplitudes lying too low, such as Fig. 5.6 (upper figure) whereas the signal is clearly spread over neighboring Fourier atoms, forming a peak.

Nevertheless, for FPs as well as FNs, scrutinizing the content of the signal over blocks of neighbouring atoms appears as a natural solution to solve both problems at once. More precisely, it is unlikely that a noise-related FP activates several contiguous Fourier atoms. In contrast, a GB signal activates a dozen of Fourier coefficients. Therefore averaging over the neighborhood of atoms should allow the discrimination of a FP from genuine signal peak.

Why should we use a block decomposition ? Unstructured approaches are known to provide biased estimators of the signal with sub-optimal convergence rate in terms of quality recovery [4]. Treating clusters of atoms jointly instead of treating them individually relieves partially the bias and improves the recovery quality. Upper row : When the signal is present, there is a power excess over all Fourier atoms with frequency close to f 0 , forming a peak. However, within the peak, some harmonics have amplitudes that are lower than the chosen threshold, and thus will not be interpreted as signals, generating FN. Lower row : An isolated noise realisation jumps above threshold and is likely to be detected as a signal, generating a FP.

Notations

This time, instead of treating independently the Fourier atoms, we will jointly process all the Fourier atoms that are "close" one to another. These neighborhoods of Fourier atoms will gather only few elements -even if it means aggregating several of them to cover a GB signal.

In practice we will consider fixed disjoint sets of consecutive frequencies B = k k min ≤k≤kmax (and thus the corresponding Fourier atoms) such that :

{-K, . . . , +K} = 1≤j≤J B j with B j ∩ B j ′ = ∅ for j ̸ = j ′ , (5.4.1)
any neighborhood B j being referred to as a block which number of elements (or cardinal) is designated by |B j |.

The set of all the neighborhoods is therefore a partition of the atoms, and is noted :

B = {B j } 1≤j≤J . (5.4.2)
This partition contains |B| = J disjoint neighborhoods, and includes all the Fourier atoms. Of course, several partitions are possible : this will be investigated later on. Now, the idea is to consider that the signal is modeled not only over each Fourier atom, but also over each neighborhood B. The sparsity prior that was "finding the sparsest solution possible over the dictionary Fourier atoms" has to be adapted to the new decomposition B because the signal has to be modeled both on Fourier atoms and on Fourier neighborhoods. The updated prior is directly stated on B : "the sought signal can be represented over few Fourier neighborhoods".

In the end, this new approach remains very similar to the previous one. Introducing the joint signal amplitude over the block decomposition :

A B S (v) =    k∈B I v * I [k] • v I [k] S I [k]    B∈B =    k∈B A 2 S (v)[k]    B∈B , (5.4.3)
where the squared joint amplitude is the sum of the squared amplitudes of the atoms belonging to the block, we will seek the signal the closest to the observed data with the sparsest joint amplitude distribution over B.

Problem formulation over a block-structured domain

Exactly as in Eq. (L 1 , I) and Eq. (Joint L 1 ), we can generalize the problem REF-L 1 to the block distribution by simply modifying the term responsible for sparsity.

Being provided a partition B of the measured frequency range {-K, . . . , +K}, the block-sparse estimator of the signal requires to solve the following optimisation problem :

h = Argmin v B∈B γ B [B]A B S (v)[B] + 1 2 ⟨d -v, d -v⟩ S , (5.4.4)
which is simplified using the notations of Section 5.2 by :

h = Argmin v γ B ⊙ A B S (v) 1 + 1 2 ⟨d -v, d -v⟩ S . (L B 1,2 )
Analogously to the problem (Joint L 1 ) which was separable by atom, the problem (L B 1,2 ) is separable by block. Moreover, within each block, the problem can be solved atom by atom. For an atom k ∈ B, the solution indeed writes :

h[k] =    A B S (d)[B] -γ B [B] A B S (d)[B] d[k] if A B S (d)[B] ≥ γ B [B] , 0 otherwise , (5.4.5)
which displays the same pattern as the unstructured solution Eq. (5.3.6). We thus foresee the need for reweighting elaborating on the discussion of Section 5.3.4 ; this issue will be addressed below in Section 5.4.4. This leaves us with two questions :

1. For a given partition B, what is the best choice for the threshold ?

2. What is the best choice for the partition B ?

Setting the threshold

The process is very similar to what was done in the unstructured case in Section 5.3.3 except that this time we look for power excess on Fourier atoms neighborhoods B ∈ B instead of power excess on single Fourier coefficients. We now consider a χ 2 -test by block and will assess :

H 0 : there is no GW signal in block B ; H 1 : there is a GW signal in block B.
Since we assumed that the real and imaginary parts of n I [k] obey independent normal distributions (Section 4.2.2,5.3.3), under H 0 (noise only)

A B S (d)[B] 2 admits a χ 2 - distribution with 2C • |B| degrees of freedom χ 2 2C•|B| .
We can adopt a hypothesis test similar to the one proposed in Section 5.3.3.

Adopting as above a rejection rate ρ B , and again defining the real number x 0 by the value of the cumulative distribution function :

P(χ 2 2C•|B| ≥ x 0 ) = ρ B , (5.4.6) the hypothesis H 0 (resp. H 1 ) is adopted for block B if A B S (d)[B] ≤ √ x 0 = γ B [B] (resp. A B S (d)[B] > √ x 0 = γ B [B]).

Improving the domain decomposition : BlockTree algorithm

The remaining question deals with the choice of the partition B. This choice is sensitive as it can impact greatly the solution : there is no reasons for all partitions B of {-K, . . . , +K} to yield the same sparse signal estimators. Optimal solutions would stem from an optimization of the overall domain decomposition, which is a NP-hard problem : the computing cost of its resolution with current computers is prohibitively expensive. And yet, a suboptimal approximation of the best decomposition can already improve the signal estimation. We hereby propose a suboptimal, yet pragmatic and efficient, solution to the domain decomposition problem.

Tree-based block decompositions were introduced in sparsity-based signal processing methods so as to adapt them better to the structures of the signals to be recovered. In different contexts, tree-based block decompositions have been used in signal and image denoising (see [13] and references therein). The approach we will develop below has been inspired by [14] where the authors build an adaptive block decomposition directly from the observed data through a dyadic process. We adapted the advocated top-down building of the dyadic tree to a bottom-up process binding the fate of adjacent Fourier atoms depending on the presence or absence of signal.

GB signal bandwidth

Following the terminology introduced in Section 5.3.1, we will say that a block B ⊂ {-K, . . . , +K} is an active block (resp. inactive block) if the joint amplitude of the signal over the block is larger (resp. lower) than a given threshold.

If series of measurements contain only noise, summing over connected blocks will act as an averaging process. If the block is big enough (i.e. includes the right amount of atoms), outliers -which scarcely occur on two neighboring Fourier atoms -will vanish and the FP rate will decrease. On the contrary, if the block is too big (i.e. includes too many atoms), a nearly monochromatic signal may be drowned into the noise and there is the risk that the signal cannot be detected anymore. Since the amplitude of GB signals can be quite close to that of the noise, this case cannot a priori be excluded. At last, if a large block is tagged as active whereas it is not, then its impact on the sparse solution will be proportionally important to its size. Thus we intuitively understand that active blocks should neither be too small nor too large.

To begin with, let us consider a uniform domain decomposition, i.e. a partition of {-K, . . . , +K} made of blocks of the same size. It seems natural to select a block size similar to the width of the peak of the expected GW signal. According to [15], the TDI signal instant phase is given by :

ϕ(t) = 2π f 0 t + 1 2 ḟ0 t 2 + (f 0 + ḟ0 t)R cos(β) cos(Ωt + η 0 -λ) , (5.4.7)
where we note ϕ the instant phase, f 0 the emission frequency, ḟ0 the frequency deriva- tive, β, λ the ecliptic coordinates, η 0 the position of LISA around the sun at time t = 0, Ω = 2π/(1 year) LISA's orbital frequency and R = 1 astronomical unit. Under the assumption ḟ0 t << f 0 , the instant frequency is approximated by :

φ(t) ≃ 2πf -2πf cos(β)RΩ sin(Ωt + η 0 -λ) , (5.4.8)
which means that, for an observation duration long enough (T obs > 1 year), the observed half peak width in Fourier domain at first order is given by |2πf 0 cos(β)RΩ|. It means that a wide range of signal morphologies can correspond to a given emission frequency f 0 . The BlockTree algorithm Such a small block size |B| = 10 does not allow a efficient enough averaging process.

This calls for the possibility of merging adjacent blocks to foster this averaging. For that purpose we developed the BlockTree algorithm, which is a bottom-up approach to frequency domain decomposition :

Initialisation Start from a uniform decomposition of blocks with minimal size.

Iteration m Try to merge adjacent blocks 2 by 2 (only if they are of similar sizes) : ■ If a signal was detected in one of the blocks at the iteration m -1, but is not detected at the current one : this signal was a FP, so we keep the block resulting of the merging.

■ If a signal was detected in one of the blocks at the iteration m -1, and is still detected at the current one : this signal was not a FP, so we do not merge blocks.

■ If no signal was detected in both blocks at the iteration m-1, we merge blocks.

The detailed implementation is described in Appendix B.3, Algorithm 3. The first iteration aims at cutting drastically the number of FPs. An illustration of the algorithm behavior is shown in Fig. 5.8. The parameters that are usually used for this code are summarized in Appendix B.3.2.

The BlockTree approach brings a major improvement : we can set a threshold linked to a probability much lower than the one chosen in the case of unstructured sparsity since we are able to discard FPs, paving the way to the detection of more GB signals. The easier detection of low amplitude signals enhances the robustness of signal detection. This point will be quantitatively assessed in Section 5.5.

Reweighted Block Sparsity

We now adapt the discussion of Section 5.3.4 to block-structured data. Reweighting will not participate to the tagging of a block as active or inactive but it will just correct the bias introduced4 by Eq. (5.4.5). There are two main options to extend the previous discussions to frequency blocks :

1. Set a unique threshold for each block.

2. Set a Fourier atom-dependent threshold within each block.

The second option arises from the observation of amplitude variations between consecutive peaks (see Fig. 5.6).

Block threshold reweighting : Similarly to Section 5.3.4, a reweighting procedure can be implemented to correct for the bias induced by the proposed block thresholding. This is described by following steps :

Phase 1 Using γ m B [B]
, compute the estimator h m at step m from Eq. (5.3.6). Phase 2 Using h m , compute the new threshold γ m+1

B

[B] as follows :

γ m+1 B [B] = (γ 0 B [B]) 2 κA B S ( h m )[k] + γ 0 B [B]
.

(5.4.9) Following Section 5.3.4, κ is a positive real parameter that allows to control the strength of the reweighting procedure. It is again fixed to κ = 3 for the block threshold procedure.

The initial block threshold γ 0 B [B] is chosen as a threshold for a χ 2 2C•|B| test as explained in Section 5.4.2. With the same notations, the choice of the threshold is given by

γ 0 B [B] = √ x 0 .
Atom reweighting for block sparsity : A potential drawback of the global block reweighting is that the correction factor, which appears in Eq. (5.4.9), is identical for all atoms of an active block. Therefore, this might not be as effective as the entrywise reweighting introduced in Section 5.3.4 to unbias the thresholding procedure. As shown there, an entrywise reweighting scheme allows to adapt to the amplitude of individual entries of the estimated signal h. In order to get the best of both the block-based approach and the atom-based reweighting, we further propose to perform both alternately. In brief, a first estimate h is computed using the block-based thresholding procedure. This allows to carefully account for FPs, and identify active/inactive blocks. The exact same atom-based reweighting introduced in Eq. (5.4.9) is further applied to the atoms of the active blocks. This allows to preserve the detection performances of the block-based procedure while significantly enhancing the quality of the detected signals.

Thus, two rejection rates are to be considered :

1. the first one is associated to the BlockTree, and therefore will be used to sort out real signals and FPs. Typically, we will choose ρ tree such that there is less than 1 FP for the considered input size.

2. The second rejection rate is used for the atom-based reweighting. Since we already selected the active atoms, we will chose a rejection rate rather high, typically ρ reweighting = 0.9 (we would like to keep active most of the atoms that have been pre-selected with the BlockTree).

Performances Benchmark

The code that we used for this study is released as an open source code. More information can be found in Appendix B.1.

Quality checking tools

Several criteria will be used to assess the quality of the solutions. The detectability of a signal, and its recovery, greatly depend on the noise realisation. Thus we will benchmark our algorithms with multiple noise realisations.

Rejection Rate : In the remaining of the article, we will refer to the rejection rate ρ defined by : P(H 0 is rejected under H 0 ) = ρ , (5.5.1) ρ is defined through Eq.(5.3.7) for the unstructured sparsity, and through Eq.(5.4.6) for the block-structured sparsity.

The rejection rate describes the probability to reject H 0 for a given frequency or a given block whereas H 0 is in fact true. For a uniform domain decomposition (i.e. a decomposition frequency by frequency or with uniform blocks), it corresponds to the expected FP rate. Its interpretation is a bit more involved in the case of an adapted domain decomposition, as we will show below.

Signal-to-Noise Ratio (SNR) : it estimates the signal power compared to the noise power. The greater it is, the more the signal is "visible". The definition of SNR for LISA can be found in [16]. When working on noise de-correlated TDI channels, the SNR writes :

SN R 2 (v) = 4δf I k |v I [k]| 2 S I [k] = 4δf k A S (v)[k] 2 .
(5.5.2)

More information about SNR computations can be found in Appendix D.4.

Quadratic relative error : for a signal h approximating a true signal h, we want to quantify the error of approximation. We define this error E 2 (h, h) by :

E 2 (h, h) = ⟨h -h, h -h⟩ S ⟨h, h⟩ S .
(5.5.

3)

The lower the relative error, the better the approximation.

Parameter Prior

Frequency

f 0 ∼ U([1 mHz, 10 mHz])
Frequency Derivative ḟ0 = 0 Ecliptic Latitude

sin(β) ∼ U([-1, 1]) Ecliptic Longitude λ ∼ U([-π, π])

Amplitude

Computed to obtain the target SN R

Inclination cos(ι) ∼ U([-1, 1]) Polarization ψ ∼ U([ 3π 4 , 5π 4 ]
)

Initial Phase ϕ 0 ∼ U([ π 2 , 3π 2 ])
Table 5.1 -The algorithms will be tested on a set of 1000 waveforms which parameters have been chosen according to the the laws reported here.

meaning that any wanted SN R wanted for signal v can be achieved through multiplying v by the amplitude h 0 computed by :

h 2 0 = SN R 2 wanted SN R 2 (v) .
(5.5.8)

We will often assess the performances of the algorithm over the set composed of waveforms at fixed SN R. The noise will be generated according to the process described in Section 4.2.2.

Individual versus Joint estimation performances

In order to show the advantage of a joint resolution over several channels compared to single channel resolution, we will consider the answer of the algorithms to the test set presented in Sec. 5.5.2 for signals with SN R = 50 5 . The fast waveforms produced by the LDC code [17] are restricted to a small frequency range of about 30, 000 frequency bins centered on f 0 (for the chosen simulation parameters). We assess the performances of the recovery over this restricted range. Two different solutions were considered : the first one is a separate resolution for channels I ∈ {A, E}, and the second one is a joint resolution. Both are using the bias correction algorithm (reweighting process) introduced in Sec.5.3.4.

In Fig. 5.9, we plot the average quality factor Q dB as defined in Eq. (5.5.4) over the test set for the solution obtained with a separate resolution (labelled as "mono"), and the one obtained with joint resolution (labelled as "multi").

For a high rejection rate ρ, the average Q dB is really low because the FP rate is important. Decreasing ρ greatly improves the results, with the best performance achieved for ρ = 10 -4 . Then, when ρ becomes even lower, the FN rate becomes prominent, which degrades the estimation quality -even with a bias correction such as the one introduced in Sec.5.3.4.

The joint resolution performs better than the single channel resolution as it decreases the atom FN rate. Around the value of interest ρ = 10 -4 (i.e. 1 FP at most for the considered input size 30, 000) where the best performance is reached, there is more than a 1 dB difference between the qualities of the two estimates, which makes more interesting to use joint estimation. We also observe that the multi-channel estimate is more robust to the choice of rejection rate ρ as the average solution quality is moderately impacted by ρ = 10 -4 being below the "best-case" scenario. For the following study, we will only focus on the joint resolution, leaving the separate resolution aside.

Estimation of the FP rate

We argued in Section 5.4.3 that the uniform domain decomposition (be it atom by atom or uniform blocks decomposition as in Fig. 5.8) presents a high FP rate compared to the adapted domain decomposition given by the BlockTree. To illustrate this point, we plot an estimation of the median FP rate as defined in Eq. (5.5.5) over 100 noise realisations over the whole frequency domain 6 , for various rejection rates and the three families of methods : unstructured sparsity, uniform block decomposition and BlockTree decomposition -for both structured and unstructured estimates. Fig. 5.10 shows the evolution of the FP rate with the rejection rate ρ.

For the unstructured sparsity and structured uniform block-based approaches (Fig. 5.10, left), the FP rate coincides with the rejection rate. This is expected since both domain decompositions are uniform. For the structured BlockTree domain decomposition, thanks to the averaging process, the effective FP rate is much lower than the rejection rate. Indeed, the probability threshold remains the same (i.e. the block threshold x 0 is always computed based on the same probability, its variations from one block to another are only due to the change in block size), but the final number of blocks is much lower than for a uniform decomposition, which explains that less blocks result in a FP. Unstructured block sparsity (Fig. 5.10, right) produces a lower FP rate for both uniform blocks and BlockTree, making it particularly advantageous in a context where we want to limit the FP rate as much as possible.

Evolution of quality factor Q dB with input SN R

In order to characterize the sensitivity of the proposed methods with respect to the input signal's amplitude, we investigated the evolution of the quality factor Q dB with the input SN R at constant noise level.

For a given test set as described in Sec.5.5.2, we adapt signal's amplitude according to Eq. (5.5.8) so that the signal has the target SN R. Then, we draw noise over all the atoms of Fourier domain and we add it to the signal. The measurements write :

d = h 0 (SN R) • h + n .
(5.5.9)

We apply the algorithm of reweigthed unstructured sparsity Eq.(Joint L 1 ) and block sparsity Eq.(5.4.5) combined with both a uniform block decomposition and a BlockTree decomposition, performing either block reweighting or unstructured reweighting inside the blocks. We set the detection threshold based on Fig. 5.10 such that the false positive rate is inferior to 1 for the considered input size (around 30.000 frequency bins).

This corresponds to :

Unstructured sparsity : ρ ∼ 10 -5

Uniform block sparsity : ρ unif orm ∼ 10 -4

BlockTree sparsity : ρ tree ∼ 10 -3 Fig. 5.11 shows the quality factor of recovered solutions for various rejection rates around these values : ρ ∈ {10 -5 , 10 -7 }, ρ unif orm ∈ {10 -3 , 10 -4 }, ρ tree = 10 -3 . For the lowest rejection rate (lower row), all methods did not present any FPs. In this case, we can see that the block decomposition methods are pretty much equivalent, and provide a better result than the unstructured sparsity method. 2. In absence of FPs, atom reweighting and block reweighting give results of similar quality.

Block Reweighting

3. In absence of FPs, uniform blocks decomposition has a recovery quality that is similar to the one obtained using BlockTree.

4. Even for a high rejection rate ρ tree = 10 -3 , the BlockTree decomposition never presented a FP.

In conclusion, the best solution is to use the BlockTree approach to estimate the signal. Indeed, since it is a block-structured approach, its results are better than when using an unstructured approach. Moreover, BlockTree enables to use higher rejection rates than uniform blocks decomposition for lower FP rate.

FN rate

In Fig. 5.12 is represented the evolution of the average FN rate (defined as in Eq. (5.5.6))

with the input SN R. The averaging is realised over a set of 1, 000 signals chosen as explained in Section 5.5.2 and amplitudes computed as in Eq. (5.5.8). One could be surprised by the FN rate of the BlockTree algorithm, that can seem rather high. While building the algorithm, we made the choice to favor having the less FPs possible instead of having the less FNs possible, since we would rather have fewer but safer data than more but more uncertain data. This explains that our FN rate is not the best among all techniques that we tried, whereas the FP rate shown in Fig 5 .10 goes in favor of our algorithm robustness.

Sensivity to the minimal block size

As explained in Section 5.4.3, the block size can have an impact on the signal recovery realised by the algorithm. Fig. 5.13 shows the average factor quality Q dB sensitivity with respect to the choice of block size |B|. The averaging is done over a set of 1, 000 signals chosen as in Section 5.5.2, with SN R = 50. We set the parameter ρ tree = 10 -4 according to the value found in Sec.5.5.4 for the considered input size. For a block size inferior to 32 atoms, the block size only has little impact on the quality of estimation. In order to limit the impact of a block if it where to be a FP, it is better to have a block size that is not too wide. But we still want it to be large enough to detect -through averaging process -low-SN R signals. For that purpose, we selected |B| = 10 as a minimal block size.

Application to LISA Data Challenges

The LISA Data Challenges (the datasets and their description) can be found at [17]. LDC1-3 deals with ten verification GBs. The corresponding LISA signal, as expected for the TDI channel A as a function of frequency, is displayed on Fig. 5.14.

PSD estimator

For the LISA Data Challenge 1 (LDC1), the noise was created using the LISACode noise simulator [17]. This code simulates a Gaussian noise in time domain, then realises a sub-sampling and a high frequency filtering in the spirit of what will be realised for LISA real data. The LDC code provides an estimation of the noise PSD (which we will refer to as theoretical PSD and denote S th in the following) that does not take the last two processing steps into account (i.e. sub-sampling and filtering).

The proposed methods are highly sensitive to the PSD (through both data whitening and threshold choices). Consequently, the straightforward use of the theoretical PSD did not give satisfying results. Indeed the real and imaginary parts of each signal in Fourier domain whitened by the theoretical PSD is expected to obey a standard normal law N (0, 1) (this is true because only few GBs are present : they do not impact much the global noise distribution). It turns out that at high frequencies the noise distribution of frequency blocks does not follow the expected χ 2 distribution as several instrumental effects are not taken into account.

The MAD (Median Absolute Deviation) estimator :

σ M AD = M edian(|X -M edian(X)|) s M AD , (5.6.1)
with X a window-sample of the initial data and s M AD ≃ 0.674, is an empirical estimator of the dispersion of stochastic distributions, which has mainly been advertised in robust statistics [18]. We use it here to estimate the standard deviation for the real and Without the aforementioned sub-sampling and high frequency filtering effects, this estimator should roughly be equal to 1. As made manifest by Fig. 5.15, this is not the case in practice, and the deviation from the expected value is significant enough to impact the behavior of our algorithm. Thus, for f < 0.027Hz we computed the MAD estimators of the standard deviation of

ℜ(d A [k]), ℑ(d A [k]), ℜ(d E [k]) and ℑ(d E [k]
) using a sliding sampling window of 5, 000 points. We restricted ourselves to f < 0.027Hz because higher frequencies are too much impacted by the high frequency filtering -this is not a problem since the signals we are looking for have lower frequencies. The size of the sliding window was chosen by finding a trade-off between a window big enough to have an important averaging and a window small enough to have an accurate representation of the correction.

Each of these four estimators can be fitted by a polynomial of degree 4 to enforce smoothness of the PSD, and each of these polynomials can be reliably approximated by their average σ pol described in App. B.4. The impact on the whitened noise distribution is illustrated in Fig. 5.16. To summarize, instead of using the theoretical PSD S th , we used the effective PSD S ef f defined by :

S ef f [k] = σ 2 pol [k]S th [k] , (5.6.2)
to obtain the results presented in the following sections. We designed our PSD estimator in this case taking into account that the amount of GB signals was negligible compared to the amount of noise. Our estimator of the noise PSD directly works in Fourier domain, and is hardly sensitive to the GB frequency peaks of the considered dataset. Only high frequencies are impacted by filter distortions. Here, we limited ourselves to f < 27mHz because the evolution of the PSD correction could be modeled rather precisely with a simple polynomial function -which was not the case if we took more frequencies into account.

LDC results

The algorithm parameters used in our code to compute the following results are summarized in Appendix B.3.2. Fig. 5.17 shows the results that can be expected from the structured sparsity method combined with the BlockTree algorithm. The left column features in upper row the solution found in the Fourier domain in red and the target (real signal without noise) in green : all frequency peaks were recovered, and no FP was detected. The lower plot shows a zoom of the solution on the last frequency peak from the left. This gives an insight on the algorithm behavior : all signals above noise level were recovered. However, this also shows that we have a poor representation of the waveform tails 7 .

The right column upper row represents in the time domain the corresponding solution in red and the target in green. The lower plot is a zoom on the solution on a day portion. The phase of the solution matches really well with the one of the target, but 7. By "tail" we mean the signal's structure outside the main peak, which amplitude quickly goes to 0.

Peak SNR Q dB (ρ tree = 10 (10). For the chosen rejection rate ρ tree = 10 -3 , no FP was detected. We added the SN R as an indicator to match the results with what was found in Fig. 5.11. the amplitude of the extracted signal is still weaker than the one of the input signal. This is to be expected since the basis that we use do not represent the waveform tails in Fourier domain : the energy loss is unavoidable.

Peak-to-peak Q dB We assessed the performances of the algorithm through estimating a peak-to-peak and a global Q dB .

We define the peak-to-peak Q dB as follow : ■ We window both the noiseless signal and the extracted signal around the i th peak of interest in frequency domain : h I,(i) , h I,(i) . The window in Fourier domain is defined as follow :

1. We chose a probability threshold (p window = 0.9 for instance).

2. We compute the bandwidth over which the signal exceeds noise level with certainty p window 3. We select the regions of the whitened signal that rise above this threshold.

Each region corresponds to a frequency peak neighborhood. Thereafter, a frequency peak will be considered as "detected" only if there is at least 1 detected frequency in the detection region defined as above.

■ We compute the quality factor based on these partial signals :

Q dB (h I,(i) , h I,(i) )
The results are presented in Tab. 5.2 for ρ tree = 10 -3 according to what was found in Section 5.5.4. For the chosen settings, all the sources are detected without any FP. Our main objective is thus completed. Moreover, based on the sources SN R, the quality of estimation for each individual source matches with the prediction made in see Section 5.5.5 : the recovery algorithm worked just as planned.

Residuals study

The residuals for the channel I ∈ {A, E} is defined by :

r I = d I -h I .
(5.6.3) Fig. 5.18 represents the distribution of |A S (r)[k]| 2 among all frequencies k such that |kδf | < 0.027Hz 8 . The residuals seem to follow is in good agreement with a law of χ 2 4 , which is expected since it is the law followed by the noise (and we only subtracted a few number of frequencies).

Conclusion and Prospects

Based on a simple signal modeling, we are able to detect and retrieve the imprint of GBs at a very low computational cost. The proposed algorithm is robust to noise, with controllable FP rate and can detect rather low-SN R sources thanks to an averaging process. Moreover, it can deal with several GB sources superposition ; thus it is a direct answer to the topics of signal detection and subtraction. We proved on realistic data (LDC 1-3) that we could detect and extract signal even in presence of complex instrumental noise. This model, though very simple, has a direct application when part of the data is missing ; this will be developed further on in the next chapter.

8. We focused on a frequency range containing a priori all GB signals. [2]. We could also contemplate using this approach to massively remove glitches from the data by choosing an adapted dictionary -or even learning an appropriate dictionary based on the data collected during LISA PathFinder. Learning the best representation is still a possibility as this theory can also adapt to the case where the dictionary has redundancies [4]. Finally, we could study the usability of such a flexible framework to search for unmodelled or hard-to-model sources like EMRIs (see Section 3.1).

Yet, many elements could be improved : for one, the present modeling does not allow to separate GB signals that are close. Indeed, the model is not precise enough (as it does not represent the structures intrinsic to the signal, neither does it represent the waveform "tails"). Another element is that there is no visible way to identify GBs individually -this is also related to the simplicity of the model. We will show in Chapter 7 that we can build a reliable non-parametric model precise enough that it could be used in a context of identification. A third point is that the proposed algorithms rely on the knowledge of the noise PSD -which will be unknown in reality. A process alternating between noise PSD estimation and signal estimation could surely lead to promising results and would not change much the current algorithms.

Chapter 6

Inpainting for LISA gapped data

Introduction

In Chapter 5, we introduced a novel, general, non-parametric framework for the detection and recovery of GBs based on a sparse representation of the sought signal. Being non-parametric, it yields a fast low-bias estimate of the GBs signal, without prior knowledge of their number. This new tool permits a precise detection of GBs, with an accurate control of the false discovery rate, which makes it an effective approach to robustly deal with the noise that contaminates the LISA data.

In the field of computer science, dealing with gapped data has long been considered in the framework of sparsity-based signal processing, leading to sparse data inpainting methods (see [1; 2] and references therein). In this chapter, we therefore propose an innovative method that combines the non-parametric GBs recovery method we introduced in [3] and sparse data inpainting to mitigate the impact of data gaps on LISA science. Elaborating on our prior knowledge of the ungapped noise distribution, we estimate not only the missing signal but also the missing noise. We show that it can efficiently recover missing data ; based on a simple example we also show that it can correct posterior distributions for parameter estimation. The work presented here was published in MNRAS [4].

The general context of our work, the corresponding framework and data modelling are presented in Section 6.1. Section 6.2 describes the two algorithms that we developed to mitigate the impact of gaps on LISA data -detailed proofs and information can also be found in Appendices C.2,C.3,C.4. Section 6.3 focuses on assessing the performances of the two algorithms in various configurations. Last, Section 7.5 draws conclusions and prospects over the present study.

Gaps and data analysis

One of the main challenges we will have to face to perform data analysis is the presence of data gaps [5]. In this section, we present the characteristics of the expected gaps, and we show that they have an important impact on GB identification.

Different types of data gaps

The document [5] gives a detailed list of the elements that could impact the robustness of data analysis. Among them, data gaps will happen in various ways that can be sorted into three categories :

Long, unplanned gaps : During the LISA Pathfinder mission, a 5-day-long data gap occurred, during which the device went into safety mode and could not collect data. This can also happen during the LISA Mission : the device will not always be able to collect data. The mission duty cycle is estimated around 75%, which means that there will be around 25% of missing data in the form of long duration gaps.

Maintenance gaps :

Regularly the antennas need to be realigned. This process is expected to take place once every two weeks and last approximately 7 hours. No measurements can be taken in the meanwhile.

Short, unplanned gaps : It is still possible to have daily short events that will cause data loss. Such events are likely to happen every day, over a short time period (from a few seconds to a few hours). For instance, this could be the result of "gapping" the signal in presence of data glitches.

These gaps will impact data analysis both on the short and the long run. This will deteriorate the observability of any GW phenomenon. For instance, transient signals like MBHBs are mainly visible during merger ; yet it is very likely that a gap -even short -could happen at this very time. The SN R loss would be important and would degrade the identification quality [6]. For continuous signals like GBs, this is the quasi-periodicity of the gaps and the cumulated amount of lost data that will impact their SN R. Fig. 6.1 represents the noise amplitude in Fourier domain for the ungapped signal (labelled as 'real') and the gapped signal (labelled as 'gapped'). The presence of gaps created a distortion in the noise PSD and added correlation between the different frequencies (whereas the noise is initially supposed not to contain any correlation between frequencies in Fourier domain).

Impact of gaps on data Fourier transform

The presence of gaps in the data leads to a significantly more challenging problem :

Loss of signal power : In presence of gaps, the waveform is distorted and is not sparse anymore : Fig. 6.2 shows the waveform resulting from a gapped signal in the Fourier domain. Even if the main part of the waveform seems unchanged (Fig. 6.2, lower row), there is a gap-dependent power loss -not visible on the picture -that is due to the information originating from the incompleteness of the data. In practice, this can lead to inaccurate detection and identification. Noise leakage : Fig. 6.3 shows how gaps impact a signal simulated for the first LISA Data Challenge [7]. This experiment illustrates how gaps also result in a noise leakage, and radically change the statistics of the noise in the Fourier domain.

With respect to noise, gaps again hamper the detection and identification of GBs as a significantly larger amount of GB signals will be drowned in noise in the case of a gapped signal.

To mitigate the impact of gaps on data analysis, several strategies can be thought of : Do nothing : Fourier-transform the data as such and conduct the usual analysis based on fast waveforms. It is the worse strategy possible : we will show in the next paragraphs that gaps can bias parameter estimation.

Apodize the gaps : Smooth the gaps, Fourier-transform the data and conduct the usual analysis based on fast waveforms. Unfortunately, apodization does not completely mitigate the impact of gaps, it can still lead to biased parameter estimation.

Work in time domain : Conduct the analysis directly in time domain using waveforms in time domain. This way, taking the gaps into account becomes much easier (as we just have to apply a mask to the waveforms). It is the most accurate strategy among all since the missing data can be taken exactly into account during the analysis. However, conducting the analysis in time domain has a computing cost much higher than when doing it in Fourier domain -at least for GBs : this can be prohibitive.

Correct the data :

Correct the effect of missing data -for instance by filling the gaps with information. Then Fourier-transform the data and conduct the usual analysis based on fast waveforms. We chose to follow this approach because it allows a good compromise between analysis precision and computational cost.

Within the last category, the bayesian "data augmentation" algorithm [8] was developed in order to make it possible to proceed with data analysis even when the data is gapped. The core idea behind the method is to consider the missing data as a variable that has to be estimated just like the physical parameters. The missing data and the parameters are then jointly estimated using a MCMC process. The algorithm was assessed over gaps that can be assimilated to the small and medium-sized gaps we consider here, in presence of a single GB source.

The method we introduce here, even if it also relies on the estimation of missing data, is fundamentally different in that it is a "source-agnostic" alternative based on our non-parametric modeling of GB signals. The estimation process depends on an inpainting mechanism that radically differs from the Bayesian approach used in [8].

Gaps Apodization

To mitigate the impact of missing data on the Fourier transform, a classical approach consists in apodizing the gaps. The principle of apodization is presented in Fig. 6.4 : whereas the missing data can be modeled by multiplying the full data by a crenel mask of 0 (when the data is missing) and 1 (when the data is measured), the apodization consists in using a mask with smooth transitions between the sequences of missing points and the sequences of observed data. More data is impacted by the new mask, but this is negligible in comparison with the total amount of lost data. is, the slower its Fourier coefficient amplitudes |m[k]| will go to 0, and the greater will be the impact on observed gapped data d g [k] (See Section 6.1.3 for more precisions.). Fig. 6.5 shows the impact of crenel and smooth masks on noise distribution : whereas both present a huge discrepancy compared to the expected distribution at lower frequencies, applying the smooth mask makes the distribution closer to the expected one above 2 mHz compared to the crenel mask that distorts the whole millihertz frequency range.

On the other hand, Fig. 6.6 shows the impact of both masks on a waveform in the millihertz frequency band. In this frequency range, both masks produce the same deformation on the waveform (no visible difference), which information is spread over a wider frequency range than when the signal is ungapped. The chosen waveform has central frequency around 6 mHz, where Fig. 6.5 yet indicates that the impact of gaps is supposed to be low (no visible deformation of the noise PSD around this frequency).

Impact of gaps on parameter estimation

As explained in Section 4.1, GB data analysis is generally directly performed in Fourier domain, as a fast waveform was developed to produce directly GB signatures in frequency domain. On the one hand, this makes analysis faster as we can only focus on a few frequency bins of interest. On the other hand, taking into account phenomena that happen in time domain (data interruptions for instance, but also transient glitches) becomes much harder. In particular, trying to conduct classical parameter estimation on a gapped signal directly in Fourier domain could lead to bias the parameter estimation, as the gapped waveform does not match the full-data waveform (i.e. the theoretical waveform). This is illustrated in Fig. 6.6 : the missing data in time domain results in the waveform spreading in the Fourier domain in addition to the power loss due to the missing data points. The theoretical waveform does not match the observed waveform anymore.

Study [6] quantifies the impact of gaps on MBHB parameter estimations. Depending on when the data is missing, the source SN R can be heavily impacted making the identification harder. Nevertheless, the setup is different since for MBHBs the parameter estimation is directly conducted in time domain : it is easier to take the data mask into account than when the analysis is performed in frequency domain.

The impact of gaps on GB parameter estimations is studied in [8]. In this article, the author shows that if the analysis is directly conducted on gapped data, the estimation bias could become important, especially at lower frequencies where the data mask has more impact. Nonetheless, the study does not consider long gaps that have the most impact on the observed SN R in Fourier domain.

Here we want to stress the impact of mixed types of gaps on parameter estimation. We ran the sampling code that we developed to estimate parameters of a GB given

Crenel gaps

Smooth gaps Figure 6.7 -Posteriors obtained with MCMC parameter estimation for one GB with parameters presented in Appendix C.1, for a full, ungapped signal (green), a signal gapped with crenel gaps (blue,left) and a signal gapped with smooth gaps (blue,right) for parameters f 0 , β, λ. f 0 is estimated in milliHertz (as mf 0 ) and sin(β) is estimated instead of β. Priors are the same than in Table 4.2. All three types of gaps are combined resulting in about 28% of data loss. The real parameters are spotted with a red line and labelled as "Truth". The contour plots delimit the 95% presence probability. The presence of gaps -smooth or crenel -that are not taken into account for the analysis create a bias in the estimation. Moreover, the posteriors for gapped data are more spread than the one for full data. in LDC1-3, reported in Appendix C.1, with typical LISA noise with and without gaps. We present the final posterior found in the two cases for three parameters : the main frequency f 0 and the ecliptic coordinates β and λ in Fig. 6.7.

Using a MCMC algorithm we estimate the parameters values by the maximum of the joint posterior distributions. Fig. 6.7 shows that the posteriors found when sampling the full signal exhibit maxima a posteriori that are close to the true values. In contrast, the posteriors obtained after sampling gapped signals -be it with crenel or smooth gapsare biased. In this case, even if smoothing the gaps corrects a bit the bias on sky location (β, λ), it creates a bias much stronger on frequency estimate. Two conclusions can be drawn when working directly in Fourier domain : gaps must be accounted for in the analysis process, and smoothing them is not enough to correct their impact on the signal.

Gapped measurements modeling

The full data model was introduced in Section 4.2.1 and comes down to : Thus, the gapped data d g is defined, in time domain, by :

d = h + n .
d g [n] = m[n] h[n] + n[n] ∀n ∈ [0, N t ] . (6.1.3)
Conjointly we can define the missing data d m by : For more simplicity, we will note :

d m [n] = d[n] -d g [n] = (1 -m[n]) h[n] + n[n] ∀n ∈ [0, N t ] . ( 6
d g = M (h + n) , d m = (I d -M )(h + n) , (6.1.6)
with M the mask operator M associated in time domain to the mask time function m which effect is described in Eq. (6.1.3) and Eq. (6.1.4) and I d the identity operator. In particular, for any data u, the operator M is defined in time domain as : In this chapter, we will have to deal with both time domain and Fourier domain at the same time. In order to facilitate understanding, we would like to point out that : ■ The mask operator M is defined through time measurements, and thus operates on time variables v[n] ■ The amplitude operator introduced in Eq. (5.2.4) is defined over Fourier coefficients, and thus operates over frequency variables v[k]. ■ Identically, the inner product ⟨•, •⟩ S is weighted by a PSD S that is defined in Fourier domain, and thus operates over frequency variables v[k].

M u[n] = m[n]u[n] ∀n ∈ [0, N T ] , (6.1.7
Following this statement, we will not specify the variables when using these operators to alleviate the notations as much as possible.

In the Notation Table (p.17) we provide a synthesis of our writing conventions with the list of the most frequently used variables of this study. In particular, bold lower case letters are used for data, for estimators (denoted by a hat) and for the sequential solutions of optimization problems (denoted by sequence index p). Upper case letters (S, M ) are used for operators.

Dealing with incomplete data in sparse signal processing framework

Incomplete data and sparse modeling

Let us go back to the sparsity-based model that was introduced in Section 5.1.1. A signal x can be represented in a dictionary T if it can be written as :

x = Tα = N T j=1 α j t j . (6.1.8)
Additionally, if x is weakly sparse in T, then the information it carries can be compressed over few non-zero coefficients α j that can be determined by solving the (L 0 ) problem :

Argmin α ∥α∥ 0 such that ∥x -Tα∥ 2 ≤ σ , (L 0 )
Now we observe a gapped version x g of signal x. The mask M applies, and the measurements write :

x g = M x = M Tα = N T j=1 α j M t j .
(6.1.9)

Yet, the underlying signal x is still weakly sparse in T : there is only a change in the number of measurements that were done. Therefore, we would like to adapt Eq.(L 0 ) to the missing data context as :

Argmin α ∥α∥ 0 such that ∥x g -M Tα∥ 2 ≤ σ , (6.1.10)
Without any other hypothesis, this problem does not necessarily produce an accurate recovery of the ungapped data. However, under the additional strong assumption that the dictionary T is incoherent with the measurement domain, the compressed sensing theory [9; 10; 11; 2] states that the problem is well-posed if enough measurements are available 1 .

The incoherence principle -also known as Heisenberg uncertainty principle [12] -can be explained plainly as : "if x is sparse in T then it should be the least sparse possible in measurement domain". Applied to Fourier domain (where GBs are sparse) versus time domain (where measurements are performed and data is missing), there is an ideal maximal incoherence : a Dirac in Fourier domain -sparsest signal possible -corresponds in time domain to a sinusoid carrying information at any time.

If x is strongly sparse in T and T is strongly incoherent with the measurement domain, then a perfect recovery is possible [12]. When noise is present and x is still strongly sparse, then the recovery is not perfect anymore but the error can be bounded depending on the sparsity level of x [9]. In our case, x is only weakly sparse. We cannot say for sure that the error is bounded ; yet an approach similar to that of Section 5.1.3 experimentally conducts to a proper signal estimate. Like in Section 5.1.3 this problem can be relaxed using the norm ∥•∥ 1 instead of ∥•∥ 0 (see [13] and references therein) :

Argmin α γ ∥α∥ 1 + 1 2σ 2 ∥x g -M Tα∥ 2 2 .
(CS)

Comparing Eq.(CS) with Eq.(LASSO) announces forthcoming troubles : before, we could use the inverse Fourier transform T -1 to work solely with Fourier data x[k] instead of time data x[n]. However this is not the case anymore : the operator M acts in time domain and does not commute with T, neither with T -1 . Eq.(CS) is a formulation stuck between time domain and Fourier domain, which makes its resolution much harder. Another key element is that in CS context, σ is the noise standard deviation on the gapped signal. In our case, we have prior knowledge on the noise PSD for ungapped signal, the difference between the two is illustrated in Fig. 6.1. We cannot use Eq.(CS) as such.

LISA data specificities

As justified above, the LISA missing data model writes :

d g = M (h + n) . (6.1.11)
Applying (CS) to the recovery of GB signals from gapped data d g would write :

h = Argmin v γ ⊙ A Sgap (v) 1 + 1 2 ⟨d g -M v, d g -M v⟩ Sgap (dg-M v) † S -1 gap (dg-M v)
, (6.1.12)

where S gap stands for the gapped noise PSD in Fourier domain, as highlighted before (see Fig. 6.1 to understand how S gap differs from S).

Standard CS methods generally make the assumption that the noise statistics is known at the level of the measured data, i.e. that S gap is known. This puts to light the specificity of the problem of LISA gapped data : we have prior information on both the signal and the noise in absence of gaps. S gap could of course be computed based on the expected noise PSD S and the mask M , but the computation would be costly and the corresponding problem would still remain hard to solve. Moreover, it would be a full, dense matrix as the gaps introduce complex time correlations that have repercussions in Fourier domain.

Accounting for the proper noise distribution is particularly important as it allows one to carefully control the false positive rate of the sought-after GB signal [3]. This is however a key discrepancy with respect to standard methods since the theoretical noise statistics has to be modelled accurately in the Fourier domain. We therefore propose an innovative sparse inpainting algorithm that aims at retrieving both the ungapped signal and noise.

Inpainting LISA data

To that end, we introduce a new estimator d m of the missing data d m as defined in Eq. (6.1.4). d m includes both missing signal and missing noise. We have :

d m ∈ Ran(I d -M ) = Ker(M ) =⇒ M d m = 0 . (6.1.13)
Thus we will search for d m on the set of variables u such that M u = 0.

We rewrite the problem of estimating the gapped signal as a joint estimation of the full, sparse signal and the missing data. The new formulation remains close to that of Eq. (REF-L 1 ) and Eq.(6.1.12) :

( h, d m ) = Argmin (v, u) M u = 0 ∥γ ⊙ A S (v)∥ 1 + 1 2 ⟨d g + u -v, d g + u -v⟩ S . (6.1.14)
How does inpainting work here ? Eq.( 6.1.14) comes down to re-creating an ungapped data d g +u which noise PSD is similar to the expectation S. We then enforce sparsity of the full signal contained in the recovered full data d g +u through the ∥•∥ 1 term.

Estimating the missing data d m implies that we also estimate the missing noise. This may not be intuitive, but can be explained by the colored nature of the noise in Fourier domain. The smoothly frequency dependent PSD (Fig. 3.7) entails that the ungapped noise has an inherent structure2 in Fourier domain that has to be recovered in the process to limit leakage as much as possible between noise and signal. If the missing data is recovered, both the priors on the ungapped signal (i.e. sparsity) and the ungapped noise (i.e. the expected noise PSD) can be used. Let us notice that the resulting optimization problem now exhibits a mixed formulation with terms expressed both in time and frequency : the minimization constraint M u = 0 is in time domain, whereas the sparsity and data adequacy terms are written in Fourier domain. Consequently, it does not admit a closed-form solution and the minimizer needs to be computed numerically with an iterative algorithm, which is described in the next section.

Interpretation as a deconvolution problem

Inpainting is sometimes misconceived as a procedure that artificially creates new information beyond the available observational data. Here, we want to highlight that data ungapping can be equivalently formulated as an ill-posed deconvolution problem. Eq.(6.1.14) aims at deconvoluting the data from the mask in Fourier domain based on the aforementioned assumptions : signal sparsity and noise structure. Indeed, in time domain gapping the data using a mask m is equivalent in Fourier domain to convolving the Fourier transform of the mask m with the Fourier transform of the sought signal. The inpainting algorithms that we propose equivalently aims at correcting the observation from the impact of the mask. The huge advantage of thinking about this problem as a deconvolution one is that we can directly make use of the sparsity prior that we have on the sought signal. In fact, the application of the mask multiplicative operator in time domain, or its equivalent convolution operator in Fourier domain, results in a less-sparse observed signal. Enforcing the sparsity of the solution then allows to efficiently correct for the mask effect [1]. Considering colored noise and gapped data, and then Fourier transforming the gapped signal results in several effects : ■ a SNR loss due to the lost data and the convolution with the mask, ■ a deformation of the noise spectrum : gapping in time domain is equivalent to creating noise correlations in Fourier domain, ■ a loss of spectral resolution, also because of the convolution with the mask that spreads Dirac distributions. Solving the deconvolution problem becomes harder as the amount of missing data becomes more important and depending on the gap distribution we are facing, as shown in Fig. 6.8. It represents the superposition of the Fourier transform of a discrete sine with finite horizon and the Fourier transform of the same signal but gapped with the mask represented on the left-hand side column. The comparison between lines one and two shows that the spread of the signal increases as the gap size increases. Line two and three show the impact of the gap distribution on the final result (as there is the same amount of missing data but with different distributions in time) : for a random distribution (line three), the presence of gaps can even be perceived as an additional noise.

Algorithms & implementation

Within the scope of LISA data processing, the goal of data inpainting is twofold. On the one hand, and following the approach introduced in [3], it aims at providing an estimate of the total signal originating from GBs that is robust with respect to noise and gaps. On the second hand, data inpainting can be deemed a general pre-processing step for LISA data treatment, whose objective is to deliver estimated ungapped data. The latter more specifically emphasizes an accurate signal estimation with no power loss as well as an estimate of the noise that matches the expected statistics. Doing so, the resulting inpainted data can be used as inputs to classical identification techniques such as Bayesian inference. In that regard the sparsity framework has the advantage of being a non-parametric methodology, which works independently from the actual number of GW sources. To that purpose we hereafter introduce two algorithms that both tackle these two views of the inpainting problem.

Resolution Algorithm

The problem in Eq. ( 6.1.14) benefits from certain properties, which can be used to build an efficient minimization algorithm. Indeed, it is not only a strictly convex problem, but also a block-convex problem : it is convex with regard to the variables v while u is kept fixed and vice versa. Let us now denote the global cost function to be minimized as :

J (v, u) = ∥γ ⊙ A S (v)∥ 1 + 1 2 ⟨d g + u -v, d g + u -v⟩ S . (6.2.1)
Thanks to the block-convexity of the problem in Eq. ( 6.1.14) both variables u and v can be sequentially and iteratively updated. For that purpose, we make use of a block coordinate descent (BCD) algorithm [14] which can be summarized with the two following steps :

     d p+1 m = Argmin u M u=0 J (h p , u) , h p+1 = Argmin v J (v, d p+1 m ) , (BCD)
with initialization h 0 = 0. The sequence {(h p+1 , d p+1 m )} p converges to the sought estimators of signal and missing data ( h, d m ), the solution of Eq. (6.1.14) [15]. We now detail each of these steps.

Updating the noise

Instead of estimating the full signal and the missing data in the gaps, it can be more convenient to think in terms of full, ungapped signal and full, ungapped noise. The missing data is then the sum of full signal and full noise in the gaps. Let us introduce the estimator of the ungapped noise n. It is related to the other estimators by the following equation :

d g + d m Completed data = h + n . (6.2.2)
For the sake of simplicity, we can rewrite the proposed BCD-based algorithm (BCD) so that the global signal and the noise in the gaps are computed sequentially. To that end, From Eq. (6.2.6), we get that :

d p+1 m = (I d -M )d p+1 m because d p+1 m ∈ Ran(I d -M ) = (I d -M )(h p + n p+1 -d g ) = (I d -M )(h p + n p+1 ) because d g ∈ Ker(I d -M ) .
(6.2.9)

Then the full updated data writes :

d p+1 = d g + d p+1 m = d g + (I d -M )(h p + n p+1 ) . (6.2.10)
The full updated data d p+1 is exactly what we expect it to be : we fill the gaps of d g with the missing signal (I d -M )h p and the missing noise

(I d -M )n p+1 .
Consequently, the signal estimation step in Eq.(BCD) can be recast as the estimation of the signal over the ungapped data d p+1 :

h p+1 = Argmin v ∥γ ⊙ A S (v)∥ 1 + 1 2 ⟨d p+1 -v, d p+1 -v⟩ S . (6.2.11)
Quite remarkably, thanks to the particular structure of our inpainting algorithm, the update of the signal h p+1 can be done from pseudo measurements where the noise is inpainted. This allows to use the theoretical PSD rather than the more complex PSD of the gapped noise. The resulting update is then equivalent to the denoising problem investigated in Chapter 5. More precisely, the signal h p+1 admits a closed-form expression which takes the form of a specific thresholding applied to d p+1 in the Fourier domain.

Classical inpainting (C.I.) algorithm :

The overall algorithm of classical inpainting (C.I.) is described in Alg. 1.

Algorithme 1 : Classical data inpainting algorithm (C.I.)

Input : d g , M, S, ϵ Initialization : h 0 = 0; while max k |A S (h p+1 -h p )[k]| > ϵ do n p+1 gap = d g -M h p ; n p+1 = Argmin q n p+1
gap =M q 1 2 ⟨q, q⟩ S ;

d p+1 = d g + (I -M )(h p + n p+1 ) ; h p+1 = Argmin v ∥γ ⊙ A S (v)∥ 1 + 1 2 ⟨d p+1 -v, d p+1 -v⟩ S .
It was further highlighted that a block-sparsity regularization leads to significantly better detection and estimation precision [1; 3]. In the algorithms descriptions, we used the unstructured sparsity formulation because the notations were lighter than for the block formulation. However, the results shown later on were based on the block-sparsity formulation. We only recover noise frequencies for which the correlation length is superior to the gap size. This is not a shortcoming, as Fig. 6.1 shows that the PSD deformation is more important for low frequencies than high frequencies. Based on Section 5.4 and keeping here the same notations, the block formulation of the overall problem writes :

( h, d m ) = Argmin (v, u) M u = 0 γ B ⊙ A B S (v) 1 + 1 2 ⟨d g + u -v, d g + u -v⟩ S , (6.2

.12)

We can still split the resolution in two steps. The resolution with regard to noise is still given by Eq. (6.2.6), with the only change in the sparsity term of the signal update equation :

h p+1 = Argmin v γ B ⊙ A B S (v) 1 + 1 2 ⟨d p+1 -v, d p+1 -v⟩ S . (6.2.13)
The solution is then analytical, as in Sections 5.4.2 and 5.4.4.

Behavior of the inpainted noise in a gap

In this section, we illustrate the performances of the proposed inpainting algorithm especially on typical LISA noise. For that purpose, we consider simulated noise realizations for LISA (not containing any gravitational signal), whose statistics is described by the theoretical PSD of the LDC1-3 (see Section 4.2.2). Gapped data are generated with small unplanned gaps and planned gaps (see Section 6.3.1 for more details about gaps generation). Next, inpainting is applied to the resulting gapped noise according to Eq. (6.2.6). The solution is displayed in Fig. 6.9 with a zoom on a single gap in the time domain (upper row) ; the Fourier transform of the global solution is also represented (lower row).

In the time domain, the inpainted noise inside a gap only exhibits a low-frequency smooth variation, which already allows for a good signal extraction. Indeed, this algorithm is completely deterministic and cannot generate information lost in the gaps. More specifically, such an inpainting methods basically fills out gaps by exploiting correlations propagating farther than the size of the gaps. This entails that within a gap, the retrieved information content tends to have low frequencies, corresponding to wavelengths larger than the size of the gap (see Section 6.3.2 for more details). Higher frequency information cannot be recovered and is definitely lost. As we will see more precisely below, this results in an unavoidable power loss on the frequency spectrum of the recovered noise compared to the expected distribution, even if in Fig. 6.9 the inpainted noise seems to follow the expected noise PSD. This shows that Eq. (6.2.6) acts like a low-pass filter with regard to the ungapped data.

Finally, we have to highlight the behavior of the algorithm with the gaps size. The inpainting algorithm takes more time to converge when two phenomena occur : when the missing data rate increases and when the gaps lengths increases. It is expected as in these two cases, the deconvolution problem becomes harder to solve as more data are missing.

Even if this correction already seems to mitigate the impact of gaps in Fourier domain, Fig. 6.9 highlights that locally in time domain the noise does not follow the expected noise distribution. This is due to the fact that we can only recover lower frequencies of the missing noise. However, if one needed to perform part of the analysis in time domain after applying the correction, this last point could be problematic. In order to provide a more consistent noise inpainting procedure, we describe in the following section an extension of our algorithm.

Modified Sparse Data Inpainting

Let us recall that the objective of our inpainting approach is twofold : designing an inpainting algorithm that i) provide an efficient detection and reconstruction of the GB signal and ii) more generally yield a statistically consistent inpainted noise. In the previous section, we pointed out that a traditional sparsity-enforced inpainting does not reach the second objective as it leads to a gap-dependent noise power loss. In this paragraph, we propose extending sparse inpainting to further correct for this effect.

A straightforward approach would consist in adding a high-frequency stochastic term to this low-frequency solution within the gaps in time domain. However, this would only produce a poor solution, as the added high-frequency noise in the gaps would not be compatible with the one observed outside the gaps. This incompatibility at the boundaries of the gaps would create high frequency artifacts. The main challenge then boils down to finding a high frequency noise correction that matches both the expected distribution and the boundary condition on the border of every gap. To that purpose, we propose now a method that combines the approach of Alg. 1 with the use of a stochastic term in order to create an appropriate inpainted noise solution.

Generation of a compatible high frequency term :

To that end, for any gapped measurement v g we define the function f CI (v g ) as :

f CI (v g ) = Argmin q v g = M q 1 2
⟨q, q⟩ S .

(6.2.14)

Algorithme 2 : Modified data inpainting algorithm (M.I.)

Input : d g , M, S, ϵ Initialization : h 0 = 0 , draw n samp ∼ N 0, S 2 + i • N 0, S 2 ; while max k |A S (h p+1 -h p )[k]| > ϵ do n p+1 gap = d g -M h p + M n samp ; n p+1 CI = Argmin n p+1 gap =M q 1 2 ⟨q, q⟩ S = f CI (n p+1 gap ) ; n p+1 = n p+1 CI -n samp ; d p+1 = d g + (I d -M )(h p + n p+1 ) ; h p+1 = Argmin v ∥γ ⊙ A S (v)∥ 1 + 1 2 ⟨d p+1 -v, d p+1 -v⟩ S

Implementation

Compensating the power loss

Although Alg. 1 improves the effective PSD, the final estimated noise PSD does not exactly match the expected distribution. Indeed the estimated PSD (after noise inpainting) and the expected one are empirically identical up to a multiplicative coefficient matching the amount of data that was lost. This phenomenon is a manifestation of power loss, and is mainly related to the fact that the inpainting algorithm C.I. can only inpaint lower frequencies -this can be easily understood by looking at Fig. 6.9.

This power loss impacts the detection capacity of the C.I. algorithm since we set the threshold γ, as reminded in Section 6.1.3, with regard to the theoretical noise distribution whereas it should be set with regard to the effective noise distribution. Since the effective PSD is lower than the expected one due to power loss, it corresponds to choosing a threshold that is higher than what it should be, therefore increasing the risk of non-detection. We correct it by adjusting the noise level, computing an effective PSD S ef f as follow :

S ef f = r data S , (6.2.18)
with :

r data = 1 N T N T n=1 m[n] , (6.2.19)
the proportion of remaining information evaluated from the mask M. This first-order correction yields a decent description of the noise distribution after the use of the C.I. algorithm, but becomes imperfect as the proportion of lost data increases.

Improving the speed of convergence

From the viewpoint of optimization, the more data is missing the slower the convergence of the iterative minimization algorithms is. Limiting this phenomenon can be done by making use of the fixed point continuation or FPC [17], which advocates computing sequential estimates with decreasing regularization parameters. In the present context, instead of setting a fixed threshold γ as explained in [3], it is first set to a large value and then decreased towards the final objective threshold γ f inal . Implementing the FPC significantly improves the convergence speed.

Experimental Results

The performances of our algorithms are assessed on three main outcomes : (i) the detection capacity (comparatively to the ungapped case), (ii) the statistics of the recovered noise and (iii) the quality of extracted signal.

As mentioned above, this work could be used as a pre-processing step in the LISA pipeline, and we evaluate the performances of the algorithms with this potential application in mind. We also assess the impact of gaps by estimating the quality of the extracted signal compared to the signal we would have extracted without gaps. This part is entirely dependent on the chosen extraction method.

Finally the proposed algorithms are evaluated on the realistic simulations LDC1-3 [7] to which gaps have been added. This goes beyond the proof-of-concept stage where these algorithms are probed on test cases containing gaps and either noise only (like in Fig. 6.9 and Fig. 6.10) or both noise and one GB signal.

Gaps generation

Data gaps are characterized by two parameters : their duration L gap and their period T gap (time period over which we observe one gap). In Section 6.1.1 we described the three different types of gaps that we consider here. The numerical values used for this study are reported in Table 6.1.

Small and medium-sized gaps are frequent events that will happen for sure. To simulate them, we split the signal in consecutive blocks of length T gap . For each block, we randomly position a gap of length L gap with a uniform distribution over a time interval of duration T gap . Large gaps are related to the duty cycle -estimated to 75%. During LPF mission, the longest data taking interruption lasted 5 days [5]. Based on this observa- tion, the duty cycle is simulated as a series of long-lasting interruptions distributed over the whole observation duration. Two consecutive large gaps of length L gap = 3 days are separated by a duration given by a Poisson law 3 of parameter λ = 9 -making approximately a 25% data loss.

Noise inpainting evaluation

Independence of the solution with respect to the chosen sample

Since the M.I. algorithm Alg. 2 makes use of an extra stochastic term, we first evaluate its sensitivity to the particular drawing of n samp . The answer is partially contained in Alg. 2 itself. We explained in Section 6.2.3 that the lower-frequency part of the solution only depends on the measurements -meaning that it is independent from n sampwhereas the higher frequency component mainly depends on n samp . In order to assess this dependence, let us consider a single input of the form :

d I = M n I , (6.3.1)
with n I a noise sample following the expected ungapped noise distribution. Inpainted noise solutions are then computed with the M.I. algorithm Alg. 2 for various draws of n samp .

Fig. 6.11 displays the standard deviation of the difference between the real noise and the solutions obtained for each sample in Fourier domain, frequency by frequency, for different types of gaps. This difference has further been whitened with the inverse theoretical noise PSD. Fig. 6.11 quantitatively shows that inpainting leads to a low-pass filter effect confirming the qualitative features discussed in Sections 6.2.2 and 6.2.3. Additionally, we point out two more features of the algorithms : the presence of a cutoff frequency and the behavior of the maximum deviation to the real noise distribution.

The cut-off frequency, labelled as f c on the plot and defined as :

f c = 1 2L gap , (6.3.2)
matches the effective cut-off frequency of the filter-like behavior of Alg. 2. As the gap duration increases, the cut-off frequency decreases : the noise component can only be rightfully recovered when its half-wavelength is superior to the size of the gap. Below the cut-off frequency, the recovered noise is quite close to the real noise. Above this frequency, there is a deviation that becomes more important as the amount of lost data increases.

Recovered noise distribution

We emphasized previously that Alg. 1 and Alg. 2 could be used as a pre-processing step for filling out gaps prior to performing further analysis such as event identification.

To that end, we assess the impact of the inpainting on the noise statistics using the Kullback-Leibler divergence as a performance indicator.

Kullback-Leibler divergence (D KL ) :

The Kullback-Leibler (KL) divergence [18] measures a discrepancy between the expected and the recovered noise distributions. We normalize the noise estimator n (limit of the sequence {n p }) in frequency domain as :

η I = S -1/2 I n I . (6.3.3)
If the recovered noise follows the expected law, we expect both the real part and the imaginary part of η I to be drawn from a reduced centered normal law at any frequency f . Thus, we will measure the discrepancy of the law of the normalized noise compared to the reduced centered normal law. Under this approximation, the KL divergence writes [19] :

D KL (f ) = 1 2 σ(f ) 2 + µ(f ) 2 -1 -ln σ(f ) 2 , (6.3.4)
where σ(f ) is the measured standard deviation at frequency f and µ(f ) is the estimated expectation for the same frequency. As we only have one noise realization at each frequency f , we estimate the expectation and the standard deviation over a small neighborhood of frequencies around f . Technical details can be found in App. C.5.

For an input constituted of noise only, as :

d g = M n , (6.3.5)
we computed the solutions given by the two inpainting algorithms for various input noises n I and various samples n samp . We compared the final noise distribution to the expected distribution through the KL divergence. Fig. 6.12 represents the KL divergence evolution with frequency over 50 samples (M, n, n samp ), when different types of gaps are present, for the ungapped signal, the gapped signal and the inpainted signals. Looking at the discrepancy of the gapped signal, we note that the frequency band of interest (the mHz band) is impacted by the presence of gaps, whatever the type of gaps. Small but frequent gaps impact most the expected noise distribution (top plot). However, these are also the easiest type of gaps to deal with, as C.I. is enough to correct the noise distribution in Fourier domain. Planned gaps (middle plot), even though they are not that wide, show the limits of this inpainting algorithm. As the amount of data loss becomes more important, it cannot make up for the power loss : this explains the inconsistency between the expected distribution and the inpainted distribution. The M.I. algorithm can handle this power loss. Large gaps (bottom plot) associated with important amount of data loss, are the most difficult to mitigate as the corresponding power loss is even more sensible than for planned gaps. This case makes plain the necessity to use M.I. in order to compensate this huge power loss.

To summarize : the two algorithms help mitigate the impact of data gaps. Alg. 1 cannot compensate the power loss, but reaches a noise distribution that is -up to a multiplicative coefficient-similar to the expected noise distribution. However scrutinizing the local noise statistic in the neighborhood of a gap, one finds that it is very different from the expected distribution. On the contrary Alg. 2 enables to reach a final noise statistic matching the expected one both locally and globally, correcting the power loss along the way.

Impact on the detection capacity

The following experiment evaluates the global performance of the overall algorithm (combining inpainting and signal extraction) in terms of false positive detection rate, i.e. in proportion of signals that are detected where there is no corresponding input signal. In order to assess the detection capacity, we evaluated the number of false positive (FP) signals. We consider the same definition of FP rate as the one introduced in Section

d g = M n , (6.3.6)
and a fixed mask M , we computed the solutions given by the two inpainting algorithms for various input noises n and various samples n samp . We then computed the FP rate corresponding to these experiments for various rejection rate ρ. The results are presented in Fig. 6.13. For all types of gaps, the FP rate obtained for both inpainting algorithms is identical to the one obtained for ungapped data. Thus the combination of signal extraction and noise inpainting does not create any kind of detection artifact.

Quality of the extracted signal

We now assess the performances of the overall algorithms (combining both inpainting and signal extraction) in terms of the quality of the recovered signal. To that end, computing the quality factor Q dB (defined in Eq. (5.5.4)) between the expected signal and the recovered signal quantifies the loss that originates from the gaps and which cannot be recovered. By construction Q dB is large when the sparse estimate h is close to the real solution h. It provides insights into the global quality of the reconstructed signal. However it does not provide any information about the FP and FN rates.

Consider an input of form :

d = M (h 0 (SN R) • h + n) , (6.3.7)
where h is the signal emitted by the GB with parameters detailed in Table 5.1, received on channel I ∈ {A, E}. h 0 (SN R) is the signal amplitude computed such that h 0 (SN R) • h has the SNR we want -this process was already detailed in Section 5.5.2. We examined three fixed masks, each corresponding to a different gap situation. For a given SN R, the median Q dB is obtained as the median of the Q dB of the solutions found for a total of 10 noise and signal realizations. In Fig. 6.14 we presented the quality factor obtained for an ungapped signal and for a gapped signal using the C.I. and M.I. algorithms for the three types of gaps. As said before, each point corresponds to the median over only 10 different situations. We Right : Large gaps. This was realized over 10 samples for each point. We compare the recovery quality in absence of gaps (green, labelled as "REF") and in presence of gaps using classical inpainting (blue, labelled as "C.I.") and modified inpainting (orange, labelled as "M.I."). The envelops of corresponding colors show the repartition of all samples. In presence of small or medium gaps, the final Q dB is not affected. However, larger gaps deteriorate the quality of extracted signal.

could not afford to do more because of the high computational cost 4 . Therefore, the statistical relevance of the median can be questioned. This is why we display at the same time the envelop of the Q dB values reached for a given SN R input. Looking at small and medium-sized gaps, we see that even if the median results of Q dB seem lower for inpainted data, the envelops of the results are compatible. We can say that small and medium-sized gaps have little impact on the final Q dB : the recovered signal has the same quality for a gapped and an ungapped signal. This is true because only a small portion of the information was lost. With larger gaps (corresponding to a data loss of 25% in this case -which corresponds to realistic previson of mission duty cycle), the recovered results are impacted by the information loss. Indeed, we see that the envelop of recovered Q dB for inpainted signals are not compatible anymore with the reference.

The amount of recovered information is still acceptable ; nevertheless, the signal power loss is harder to compensate, which results in a deterioration of the extracted signal compared to the case without gaps. We can get a better idea of the efficiency of the algorithms by studying the quality of restitution for a specific signal (which parameters are given in App.C.6) in the presence of different types of gaps. Fig. 6.15 shows the median quality factor this signal set at different SN R, with different noise and gaps realizations. The phenomenon observed in Fig. 6.14 is confirmed : for small and medium-sized gaps, the signal quality is fully recovered. Large gaps impact more the quality of the result.

At this stage only, we can notice a difference between the two algorithms performances : Alg. 1 seems to perform better than Alg. 2. This performance spread can be explained by the fact that Alg. 2 reintroduced higher frequency noise in the gaps and thus results in a noisier extracted signal.

Application to LDC1-3 with gaps

We applied Alg. 1 and Alg. 2 to a more realistic case. The dataset was produced for the LISA Data Challenge LDC1-3 which can be found online [7]. It consists in 10 verification binaries that can be seen in Fig. 6.3. We modified it with the three types of gaps presented in Table 6.1, which corresponds to a loss ratio of about 27%, and we evaluated the quality of the signal extracted by both algorithms.

Results are presented in Table 6.2 and are compared to the quality of extracted signal for ungapped data (computed in [3]).

We computed a local Q dB corresponding to the quality of each extracted source, and a global Q dB corresponding to the quality of the total extracted signal. The Q dB depend on the chosen algorithm and noise configuration. The C.I. algorithm tends to better detect signals that are close to the noise level (peaks 3, 8 and 9), whereas the M.I. algorithm tends to yield better signals with high amplitude (peaks 5, 6, 7 and 10). Both missed peak 4 : too much information was lost, the remaining signal does not exceed noise level anymore. Thus, it could not be recovered. This was expected since its initial SN R was the lowest of all (SN R = 9.56) : according to previsions of Fig. 6.15, this is below detectability level in presence of large gaps.

On the overall, Alg. 1 and Alg. 2 extract signals of similar quality but that might not be the case for other types of sources.

Inpainting & parameter estimation

We now want to give an intuition of the impact that inpainting could have on parameter estimation. We do not provide a thorough study, but we will show the impact of inpainting on the example that we used in introduction in Section 6.1.1. We applied both inpainting methods to the crenel-gapped GB signal, and then estimated the parameters using the MCMC code that we developed.

We compared the sampling of the gapped and full data presented in Fig. 6.7 with the sampling of inpainted data. The results are presented in Fig. 6.16. The sampling realized on inpainted data, contrary to the one realized on gapped data, leads to posterior distributions that are very close to the one obtained by sampling the log-posterior of the ungapped signal. Looking at the sky position (β, λ), we observe no difference between the classical inpainting and modified inpainting results : both could correct the bias induced by the presence of gaps. The main difference between the two algorithms is observed on the frequency estimate : for modified inpainting, the posterior distribution has tails that are heavier than for classical inpainting. However, this effect remains marginal and cannot yet be held as an inpainting effect without being verified through a thorough study. In the end, what matters is that the maximum a posteriori is close to the true value.

Here, we understand the advantage of using inpainting algorithms to correct the impact of gaps on data analysis : using basic sparse signal representation, we are able to estimate the sum of all signals that are present without having to go through identification.

Finally, we showed that the gain can be huge on parameter estimation, as the bias that was present has been completely corrected. We can potentially correct at once all the biases caused by the gaps on every GB signal that are detected.

Conclusions

With a foreseen duty cycle of 75%, data gaps will constitute an important feature of realistic LISA measurements that may impact the scientific deliverables of the mission. To date only a handful of studies have addressed the problem of gap mitigation. We contributed to this effort by introducing a new non-parametric method, in the form of two complementary algorithms. Based on the sparsity framework introduced in [3] we showed that it is possible to fill the gaps with both signal and noise so that the recovered noise distribution matches the expected one and the signal power loss is compensated. More precisely, the proposed inpainting algorithms could mitigate the impact of gaps for all GB signals at once, without needing to identify them, for a computing time that is essentially related to the amount of missing data and not the number of sources.

We conducted extensive tests of this non-parametric approach and demonstrated the performances of both algorithms when confronted to different types of gaps. We also considered a more realistic case where more than 27% of the data was lost in presence of multiple GB sources. The algorithms achieve similar performances in terms

The algorithms depend on the noise modeling through the definition of the threshold γ with regard to the noise level. However the adaptation to the case where the noise distribution is unknown (but still supposed Gaussian in frequency domain) is straightforward and only requires an estimation of the PSD during the resolution process.

The current limitations of the proposed algorithms are mainly related to the signal extraction component. As we adopted here a model-independent approach through a representation of the signal in Fourier domain, there is natural room for improvement in the matching of the solution that we find with the expected waveforms. We anticipate a marked improvement of the recovered signal quality in presence of large gaps with the use of an adequate representation. For instance, the inpainting algorithms could be combined with other types of representations like the learned representation we present in Chapter 7.

the whole bank of templates, the reduced basis has generally a dimension that is much higher than the manifold actual dimension.

Non-linear Models & Learned Representations

Retaining the distinctive feature of the non-linear dependence of GW waveforms on their parameters, the investigation of the applicability to this field of learned non-linear representations comes as a natural approach. They intersect non-negligibly with models based on learning ; they have also demonstrated significant performance in this area [9]. In LVC, [10] already proposed to learn a model starting from the physical parameters to produce faster waveforms and thus speed up parameter estimation. Here, we consider a different approach that focuses on learning from the signal's morphology independently of the physical parameters.

Like in Section 7.1.2, we want to create a model of elements lying over a low-dimensional manifold. In this context, the most adapted network is the auto-encoder (AE) [9]. AEs are networks used in deep learning to build data representation models. They are known to successfully build low dimensional representations for complex data -like the GB signal -while ensuring an efficient recovery. More generally they have been used for data compression, data denoising and many more applications [9].

Classified as unsupervised learning, the AE principle is based on learning an encoder Φ and a decoder Ψ such that : ■ For any signal s, the low-dimensional encoded value is c = Φ(s). c is also called the latent variable or code.

■ For any encoded information, we can recover an estimate of the corresponding waveform simply by applying the decoder Ψ : s = Ψ(c). ■ Φ and Ψ are learned so that they minimize the global reconstruction error :

E 2 rec = s s -Ψ(Φ(s) c ) 2 2 . (7.1.2)
The low-dimensionality is generally directly engraved by the network "bottleneck" architecture illustrated in Fig. 7.2 : the dimension of the code c is usually much lower than the one of input s.

Nevertheless, two main drawbacks could deter their use in GW waveform representation case : these networks are known to need a huge training set (i.e. a large set of waveforms) to learn a low-dimensional model generalizing well. In GW case, the waveforms are potentially very costly to produce and to store, the size of the required training set could be prohibiting. Moreover, the more layers the network has, the heavier the training is (both in terms of complexity and computing resources) -this can quickly become prohibitive.

In the midst of all existing AE models, only few present the capacity to learn on a small-sized training set. Among them, the interpolatory auto-encoder [11; 12], referred to as IAE, is an AE designed specifically to address this last point. This is the network architecture that we selected to demonstrate how non-parametric methods can produce precise models for GW signals.

Contribution :

We developed a non-linear, non-parametric model for GB signals based on the IAE principle. We provide with a sparsity-based extension of this network the data s is encoded to a low-dimensional latent space using the encoder Φ. The low-dimensionality is ensured by the architecture which brings large inputs to small-sized encrypted data c. This specific network design is widely known as the "bottleneck" architecture. The encoded information c can then be decoded using the decoder Ψ to recover an approximation s of the input. that promotes low-dimensionality for the model. We show that this new model has its sensitivity increased compared to the sparsity-based model presented in Chapter 5, not only in terms of detection but also in terms of recovery accuracy. We illustrate many potential applications for LISA from simple detection to source partial characterization. We outline how a detection pipeline could be based on this model.

We present the IAE concept in Section 7.2 and how the new model is built based on this concept. Section 7.3 provides a thorough investigation of the model's performances and the perspectives it offers. Conclusions are drawn in Section 7.5.

Learning a representation with an interpolatory autoencoder

We describe here the training set and the principle of IAE. We demonstrate its performances in terms of waveform approximation on a training set and show that it performs better than a linear model.

Data overview : building the training set

The training set T train is composed of samples from the manifold that the model must learn to represent. To ensure that the model is generalizable, it has to cover a wide range of plausible signal morphologies that is representative of all possible signal's morphologies. To cover the morphological diversity as best as possible, we propose to build the training set based on a sampling of the physical parameters space Θ GB . The waveform selection was made according to the distributions presented in App.D.3.1.

A typical input for IAE is presented in Fig. 7.3. It is built according to the following process : 1. Setting the observation duration T obs , the time step ∆t and sets of waveform parameters θ, we produce the waveforms using the LDC-Sangria code [13] (fast waveforms). The sets of parameter θ is chosen according to a sampling process detailed in App.D.3.1. We choose to work with the TDI channels A and E.

2. For each signal, we only keep N input points of the signal. We select these data points such that the selected signal is centered, as in Fig. 7.3

3. We separate the real and imaginary parts to work on real-valued data. As this is only a demonstrator, we proceed with disjoint analysis for the different channels of information : for now, we do not force any special relation between the different channels 1 . We only learn the network on the real part of channel A and apply it to all information channels (real and imaginary parts of A and E)2 .

4. We normalize amplitudes of the waveform such that :

k |s[k]| = 1 . (7.2.1)
The signal's amplitude will be estimated when fitting the model to the data.

With this in mind, we build two sets : the training set T train to train the network and a test set T test to assess the generalizability of the network. The detailed number of elements per set can be found in Appendix D.3.2.

Interpolatory Auto-Encoder model IAE Principle

The IAE is an AE which architecture differs from the usual "bottleneck" one described in Fig. 7.2. Its principle is illustrated in Fig. 7.4 : we consider the training set T train described in the previous section from which we label a certain number N a of elements as "anchor points" (APs) (e i on the figure). Now, we aim at learning a non-linear transform Φ such that the transform of any point of the training set s, Φ(s), can be expressed as a linear interpolation of the transformed APs {Φ(e i )} i . De facto, it is like looking for an alternative domain in which the data has a linear representation. In spirit, Φ is a transform that can "unfold" the manifold M described by the elements of the training set. The selected APs work as a representation dictionary in the encoded space, also called "latent space".

In parallel, we also learn a decoder Ψ that can transform back the interpolated data into a waveform. Therefore, any waveform s ∈ T train is encoded by its barycentric coordinates {λ i } i , i.e. its pseudo-coordinates in the dictionary of representation {Φ(e i )} i . Moreover, it is possible to produce a signal which morphology is similar to that of a waveform by applying Ψ on any vector from Span({Φ(e i )} i ).

Contrary to classical auto-encoders, IAE's architecture does not have to be bottlenecked. The bottleneck is associated with a projection on a low-dimensional vector subspace ; here, the interpolation plays this part. The low-dimensionality then only comes from the number N a of anchor points that were initially selected as the dimensionality reduction happens during the interpolation process happening in the latent space.

Architecture & Learning

The global learning process is represented in Fig. 7.5 : the elements x of the training set T train are encoded with Φ and then interpolated on the transformed APs {Φ(e i )} i using the interpolator I. Then the interpolated result is decoded using Ψ, and the obtained output Ψ I Φ(x) is compared to the corresponding input x. The three main blocks of an IAE that are the encoder Φ, the decoder Ψ and the interpolator I -combined with the choice of APs -entirely define the model.

We set the encoder Φ and decoder Ψ as classical dense neural networks with the same input and output dimensions combined with the respective activation functions mish [14] and sft [15]. Both Φ and Ψ are made of two hidden layers with the same input/output sizes. We choose a network that does not enforce dimension reduction as the interpolator I is responsible for this step.

For any interpolator I, the learning process is described by Eq. (7.2.2) :

Φ, Ψ = Argmin Φ,Ψ µ x∈T train x -Ψ I Φ(x) 2 2
Reconstruction error

+ x∈T train Φ(x) -I Φ(x) 2 2
Interpolation error .

(7.2.2)

We look for Φ and Ψ as the applications that minimize the global error. The latter is split in two terms : the first one represents the reconstruction error in the direct domain, i.e. the quadratic error between the input x and the estimation given by the network Ψ(I[Φ(x)]). The second one is the interpolation error, i.e. the error made when interpolating Φ(x) on the transformed APs in the latent domain. µ is called the regularization parameter and balances the effects of the two types of errors on the general cost function. The parameters used in our code can be found in App.D.3.2. Two elements will play a crucial role in the quality of the final model : the APs and the interpolator. They have to be mutually accounted for because they jointly impact the final result. The interpolator impacts the quality of approximation in the latent domain, whereas the APs impact the quality of both the interpolation and signal recovery. Combining this fact with the constraint of low dimensionality, it is clear that a viable selection strategy is required. This selection process should answer the following questions :

1. How many APs are needed ?

2. What waveforms should be selected as APs (i.e. finding a selection criterion) ?

3. Is there an optimal choice for APs ?

We can give a partial answer to these questions. About the required number of APs, based on the idea that the transform Φ is unfolding the manifold and that all the data should have a linear representation after that, we can say that the required number of APs is directly correlated with the dimension of the manifold. Here, the dimension of the manifold is 7 (corresponding to all the physical parameters except for the amplitude). Ideally, considering the interpolation scheme that we use, we would need 8 APs. However, this is only an optimistic lower bound, as the shape of the manifold is very complex and its features hard to learn.

Regarding the AP choice, the strategy consisting in randomly selecting a small number of elements in the training set often results in a poor representation of the training set. This suggests that there is indeed a "best choice" for the AP set. We propose a detailed analysis of the stakes posed by AP selection, as well as a way to perform this selection in the next section.

Interpolation & anchor points selection

This is the interpolator that defines how the data is approximated in the latent domain. Moreover the interpolator also turns the IAE to a generative model : we can easily approximate any GB waveform only from the APs. Fig. 7.4 shows how the IAE works with a linear interpolator of type :

I(x) = Na i=1 λ i (x)Φ(e i ) . (7.2.3)
The main issue now revolves around computing the barycentric coordinates λ i , which is a priori no easy task as this is highly correlated to the number and the choice of APs. In this section, we investigate how to automatically select the APs. Such a process would automatically yield the number of APs needed to have an accurate model. It would also solve the question of the optimal set of APs. Subsequently, we will provide a process to select a small-sized set of suboptimal APs.

Direct projection

The λ i were first computed to minimize the quadratic error of approximation :

λ i (x) i = Argmin {λ i } i s.t. i |λ i |=1 Φ(x) - i λ i Φ(e i ) 2 2 .
(7.2.4)

We enforce a supplementary constraint of type i |λ i | = 1. This removes a degree of freedom and adds an implicit constraint during the learning process, making it more robust.

A direct projection on the APs gives satisfying results when only few APs are considered. Yet, without any prior knowledge on how to select them, an easy choice -even if it goes against the model's low-dimensionality -can be to consider a large set of APs chosen randomly among the training set elements. In this case having a large number of APs becomes problematic because Eq. (7.2.4) becomes hard to solve as the APs (and a priori their transforms) tend to be highly correlated. This creates large numerical errors when solving the least square problem.

Regularized projection

We can address the direct projection issue by regularizing Eq.(7.2.4). As stated before, without any prior knowledge on how to select APs, it is tempting to consider a large set of them so that the model fits the data well. Still, we do not want to compromise on the low-dimensionality constraint as it is one of the model strong points. Since we are working to reduce as much as possible the number of APs, a first step toward that goal is to minimize for every input x the number of APs needed to give a good approximation of x. In short, we will ask for each Φ(x) to be represented with few elements of {Φ(e i )} i : this is exactly the sparsity framework we introduced in Chapter 5. Thus, the norm 1 comes forward as a natural regularization for this problem.

Replacing Eq. (7.2.4), the new regularized projector write as Eq. (7.2.5) : The projection term is the same as in Eq. (7.2.4). The new right hand side term corresponds to a weighted norm 1 for barycentric coefficients {λ i } i , which weights {γ i } i are also learned during the process. Eq. (7.2.5) is still a convex problem. However, there is no analytical solution to it : we approximate the solution using the first iterations of ISTA3 [16]. From this change, two observations can be done : first, this is not a problem anymore to consider a (reasonably) large set of APs. Fig. 7.6 shows that using a set of 50

λ i (x) i = Argmin {λ i } i s.t. i |λ i |=1 Φ(x) - i λ i Φ(e i )
APs randomly chosen in the training set, the numerical problem that was observed using Eq.(7.2.4) is corrected, and the global error on the training set is greatly reduced.

This change of behavior originates from how the { λ i } i are computed : solving Eq.(7.2.4) requires to inverse an ill-conditioned matrix4 , whereas Eq.(7.2.5) does not require to inverse any matrix and is well-posed thanks to the regularization. Second, as we constrain the elements of the training set to be represented with as little APs as possible, some APs end up being used a lot more than the others, as if they were carrying more information. Some others are little to never used to represent the data, as if they were too redundant with the "main" APs.

Automatic anchor points selection

Using this new interpolator, we enforce a stage of information compression. For our training set T train described in Section 7.2.1, information compression is noticeably efficient when starting from a rather high number of APs (20-25 -this is high compared to the expected dimension, meaning that there is a lot of information redundancy) and using the interpolator Eq. (7.2.3, 7.2.5). This last observation gives a "natural" approach to conduct an automatic AP selection : we can alternate between a learning stage and an elimination stage. The learning stage consists in iteratively solving Eq. (7.

Pre-processing stage

The GB signal clearly presents several features : short wavelength variations -with variations of one order of magnitude -and long wavelength variations -the global structure of the signal in Fourier domain. Experiments showed that large scale variations are harder to learn by IAE. A simple way to deal with this is to apply a scale-aware preprocessing, for which wavelets are well adapted.

Wavelet transforms aim at giving a multi-scale decomposition of a given signal [17].

Here, we will designate the scale at level j by L j from the coarsest (j = 0) to the finest (j = J max ). From all the scales it is possible to fully recover a signal using the right inverse wavelet transform. Many fast algorithms are developed to this aim. Numerical tests showed that the results were better when the signal was decomposed over the different scales, and a separate model is learned for each scale L j . One explanation may be that the change of representation -now as wavelet coefficients -as well as the separation of scales allow to split a complex signal in "less-complex" sub-signals that are easier to learn from. There is a sort of de-correlation of some of the effects making the initial signal complex.

Then the models are combined to identify the global signal.

In our case, using the python module pywt [18] with the wavelet transform associated with sym2 (symlets of order 2) and maximal scale of 2 (so including 3 scales L 0 , L 1 , L 2 ), we learn 3 separate models that are then recombined to create an estimate of the signal. The global process in presented in Fig. 7.8 : starting from the raw data s, one has to apply the wavelet transform and then normalize each scale L j to perform the learning. In order to reconstruct the global signal, one retrieves the output of each scale model, and look for their optimal linear combination : the one that represent the best the initial signal when the inverse wavelet transform is applied. The global signal recovery thus relies on the wavelet transform linearity. We apply the inverse wavelet transform independently to each scale, and then fit the amplitudes corresponding to every scale so that we minimize the quadratic error between the initial full signal s and its reconstruction 5 . This scale amplitude fit is a low-dimensional least square : the solution can be computed really fast and does not present any issue. The performance results of the next section are computed after recombining the different scales.

Let us highlight that the global model is still low-dimensional since we constrain all the scales to be as low-dimensional as possible through sparsity (as in Eq.(7.2.5)) and we impose the APs to be identical for all the scales.

During the benchmark performance, we will designate by "multi-scale " the IAE model based on this preprocessing ; the model based on the full data (without this pre-processing) will be labelled as "full". The multi-scale model has 3 scales (from j = 0 to J max = 2) and uses 9 APs that are the same for all the scales. The full model uses 28 APs.

Efforts will be made to show the pros and cons of each model. We will discuss the pros and cons of each model in the following sections.

Use, applicability & performances

We will show how the IAE model is used in practice and what are the performances that can be expected from it. We will then demonstrate the performances of these two models in various situations : waveform approximation, detection, denoising, but also source separation. When it calls for, we will compare the results of the IAE models with the ones obtained using the unstructured block sparsity (combined with BlockTree) described in Section 5.5.

The performance indicators remain the same as those introduced in Section 5.5.1 :

the quality factor Q dB (Eq.(5.5.4)), the false positive rate (FPR, defined in Eq.(5.5.5)), the false negative rate (FNR, defined in Eq.(5.5.6)), and the SN R (defined in Eq.(5.5.2), effective computation in App.D.4).

Using IAE to approximate waveforms Use & interpretation of IAE

In this section, we assume that the IAE model {Φ, I, Ψ} has already been built and learned following Section 7.2.2. We demonstrate how to use this model on a practical example.

Any GB waveform s can be approximated by finding the amplitude A BSP and the barycentric coordinates λ BSP i i encoding at best the information contained in s :

A BSP , λ BSP i (s) i = Argmin A,{λ i } i s -A • Ψ i λ i Φ(e i ) 2 2 = Argmin A,{λ i } i ∥s -A • Ψ ({λ i } i )∥ 2 2 .
(BSP)

This process is called barycentric span projection. Since the network was learned on normalized data, we have to jointly estimate the amplitude A BSP of the signal.

What does Barycentric Span Projection do ? It performs a projection of a signal s on the span of the APs in the latent domain.

De facto, Eq.(BSP) is solved using JAX [19]

Waveform approximation

APs play the same role as dictionary atoms in Section 5.1.1 but contrary to our previous model, the associated transform is not linear. Choosing the APs (and limiting their number) is like selecting the atoms representing best the signal (i.e. resulting in minimal recovery error). In a similar way, if we were considering a linear model with the same restrictions (low-dimensionality, best representation possible) then the best set of "anchor points" would be given by PCA components associated with its L largest principal values [2].

The IAE model used here has 9 APs and includes a preprocessing step for the data that is detailed in the next sections. In order to have a reference for performance assessment, we will compare IAE to the waveform approximation given by projecting the training set and the test set on the L = 9 first vectors given by the PCA of the training set. Since PCA is the "best linear approximation based on L components" 6 , any linear approximation based on the same number of components is likely to give an approximation of lesser quality : this we will show by projecting directly the training set and the test set on 9 elements randomly chosen in the training set. This last approximation will be referred to as "APs" 7 . Fig. 7.9 shows a typical example of waveform approximation by PCA (using linear projection) and IAE (using BSP). As expected, we observe that IAE could seize details that PCA was unable to because of the low number of PCA components that we considered. This is exactly the phenomenon that we explained with Fig. 7.1 : a low-dimensional linear model is not adapted to the approximation of a manifold with a complex geometry. This phenomenon can be quantified at larger scales. Fig. 7.10 shows the distribution of the quality estimator Q dB , defined in Eq. (5.5.4), of the recoveries of the training set's elements and the test set's elements. The performances of IAE, PCA and projection on elements of the training set ("APs") are represented. As stated before, the "APs" ap-6. In terms of approximation quadratic error. 7. It is like approximating the data by a linear combination of the APs in the direct domain using a linear projection. proximation gives lower quality estimations than PCA as it is a linear model that was not specifically designed for this data set. PCA already shows improvements in comparison : there is a heavy distribution tail of elements having a better recovery. However, PCA is much less efficient than IAE : there is an order of magnitude between the mean approximation error made by PCA and the one made by IAE for the same number of components. Moreover, the distribution of the quality estimator is more homogeneous for IAE than PCA approximation quality. This means that when constraining the number of PCA components to be as low as possible, PCA cannot give a good approximation of the waveforms. We observed that considering more PCA components reduces the approximation error. This is expected because we increased the model dimension. But there would have been less information compression in that case, which can prove to be problematic when we consider noisy data : in a context of trying to estimate a signal from noisy measurements, low-dimensionality helps limiting over-fitting.

Waveform denoising

Linear model versus non-linear model

Our initial objective was to characterize a GB from noisy measurements. In this context, it is central to have a model that is robust against noise. We tackle this robustness by considering inputs of form :

d = h 0 (SN R) • h + n , (7.3.1)
where h is an element of the test set, and h 0 is an amplitude computed such that the signal h 0 •h has a targeted SN R with regard to noise n. This process was already detailed in Section 5.5.2. We try to estimate the underlying signal for various SN R. We study comparatively the quality of the estimated signal made by a low dimensional IAE model (based on multi-scale decomposition and 9 APs, like in Section 7.3.1) and the two linear models already mentioned in Section 7.3.1. We compare the performances of IAE, PCA and APs 8 in this context for a test set with three SNRs : SNR ∈ {20, 40, 100} ; a typical input is shown in Fig. 7.11. Fig. 7.12 presents the distribution of the quality factor Q dB for waveform estimations given by the three models in the proposed cases. First, we notice that PCA and APs provide signal estimates which quality do not improve as the SNR increases. In fact, the quality of estimation remains strictly the same : this is the proof that the error is due not to the noise but to the fact that the model represents poorly the sought signals. On the contrary, the IAE has a performance that improves as the SNR increases. Moreover, it performs much better than both APs and PCA : in the "worst case scenario" (SNR = 20), the IAE-estimated signal has in average a quality 3 times better than the PCA or APs one. The impact of this difference of recovery quality is illustrated in Fig. 7.11 : for a noisy input at SNR = 20, we can see the recoveries given by IAE and PCA compared to the true signal. While IAE remains really close to the real signal, this is not the case for PCA.

The main conclusion here is that non-linear models can efficiently reproduce a complex signal. We surely could use more PCA components to approximate the signal. However, PCA vectors are hardly interpretable in terms of parameter estimation. On the contrary, the organisation of the latent space is highly correlated to the signal shape, and thus directly to the physical parameters (see Section 7.4.3 for more information). This could be a way to perform a fast parameter estimation on the signal -this work is still ongoing, and is presented in Section 7.4.3.

Sparse modeling versus IAE modeling

One of the goals of changing the waveform model from a sparse representation to a non-linear model was to improve the signal representation. This improvement is related to two levels of information : first, representing the signal itself. Second, being able to characterize the system from the model. Here, we tackle the first level. The second level is probed in an exploratory work presented in Section 7.4.3. We consider the same test set which elements are transformed as :

d = h 0 (SN R) • h + n , (7.3.2)
such that all elements have a fixed SNR. We study comparatively the answers of sparse representation versus IAE models -both full and multi-scale-in terms of the recovery quality Q dB for both channels A and E.

For the sparse representation, we selected the combined approach "atom reweighting for block sparsity" described in Section 5.4.4 (labelled simply as "sparse" here) .

We have to highlight that IAE models were learned over the frequency range [2, 3] mHz, which is more restricted than for the sparse representation benchmark (the main frequency was in [1, 10] mHz). Hence the sparse representation benchmark realized here having slightly different results than the one presented in Section 5.5.5. algorithm (labelled as "Sparse", blue triangles), for classical IAE (labelled as "IAE Full", green squares) and for multi-scale IAE (labelled as "IAE MS", orange circles). Fig. 7.13 shows the evolution of recovery quality with the input SN R. The quality of the sparsity-based approach increases linearly with SN R for these types of inputs. In comparison, IAE models are characterized by two behaviors. The IAE full model also has a recovery quality that improves linearly with the input SN R, faster than for the sparse model. For input SN R at 110, there is already a 5 dB difference of quality in the recovery. Looking at this statement, we can conclude that our objective of building a model that represents waveforms better has been fulfilled. On the other hand, the multi-scale IAE performs as well as the full IAE at low SN R but falters at higher SN R. This is due to the fact that we reach the model's limit : because of its low-dimensionality, fine details of the signal escape an IAE model with a low number of APs. We foresee a tradeoff between the low-dimensionality and the accuracy of the IAE representation. This remark concerns only the IAE architecture and not the principle of the detection algorithm. A bias clearly appears at high SN R and dominates the approximation error.

Detection algorithm performance IAE behavior in absence of signal

The IAE can provide a signal estimate for any type of input, as noisy as it can be. This is still true in the extreme case where the input consists in noise only. This is illustrated in Fig. 7.14 : for an input consisting in noise only, IAE is still able to estimate a signal with maximal likelihood of presence. However, its amplitude is close to noise amplitude : the detected signal is merely noise over-fitting. This gives us the right quantity to consider to tell apart noise and signal : the SNR of the estimated signal. We will call it "estimated SNR" and denote it by SNR est in the following study. Going back to an input containing signal like in Eq.7.3.2, it is now interesting to analyze the behavior of SNR est depending on the SNR of input signals. In Fig. 7.15 we plotted the distributions of SNR est , the response of IAE for inputs with various SNRs (SN R = 10, SN R = 5 and SN R = 0 -this corresponds to noise-only inputs). In order to get the target SN R on the test set, we proceed similarly as in Sec. 7.3.2.

Let us first look at the IAE response to noisy test sets with SN R = 10 and SN R = 5. Even if all input signals have the same SN R, the repartition of the estimated SN R SN R est on the output corresponds to a normal distribution with a mean and a standard deviation depending on the input SN R. This spread is due to the fact that the identification is not perfect because of the noise. The lower the SN R, the more biased is the mean of the estimated SN R with regard to the input SN R. This behavior can be explained : when the SN R is low, the signal amplitude is close to the noise level.

In this context, it is harder for IAE to separate the signal from the noise, and the risk of over-fitting noise becomes more important. Noise over-fitting artificially increases the recovered signal SN R (compared to the input SN R), this explains why SN R est is usually higher than the input SN R when signals are close to the noise level.

The aforementioned case where the signal is totally absent (i.e. SN R = 0 for the input) corresponds to the blue distribution on Fig. 7.15. We observe in Fig. 7.16 that the distribution of SN R est can still be assimilated to a Gaussian distribution SN R est ∼ N (µ est , σ 2 est ) which mean µ est and standard deviation σ est can be estimated :

µ est = 5.992 , σ est = 0.5356 , (7.3.3)
for the multi-scale model -we found a similar distribution in the full model case (µ est = 6.067, σ est = 0.5423). We can use this distribution to set a detectability threshold above which we are practically certain that there is a signal. The multi-scale IAE and full IAE models yield close distributions with nearly equal selection thresholds. This characteristic no-signal behavior can be used as a founding stone for a detection process. It is based on a hypothesis test similar to the ones we conducted before in Chapter 5 and it would confront two hypothesis : H 0 : There is no signal.

H 1 : There is a signal. The choice is made according to a p-value statistical test : we set a probability threshold corresponding to the expected false positive rate p FPR when trying to distinguish a signal from noise. Then the corresponding limit SN R lim is defined such that : It can be interpreted as the limit of detectability upon which we accept or reject H 1 :

P (SN R est > SN R lim ) = p FPR . (7.3.4)
Accept H 0 , Reject H 1 ⇐⇒ SN R est ≤ SN R lim , (7.3.5) Reject H 0 , Accept H 1 ⇐⇒ SN R est > SN R lim . (7.3.6)

False negative rate

The false negative rate (FNR) quantifies the detection capacity of the model for the proposed hypothesis test. This rate will depend on the input SN R. To assess the FNR, we apply IAE to the noisy test set with SN R varying between 0 and 30, and we apply the hypothesis test proposed above to the IAE output with a threshold probability p FPR = 0.01. Then we compute the proportion of undetected signals. For the chosen p FPR , the limit SN R is :

SN R lim = 7.238 .

(we found SN R lim = 7.33 for the full model). This detection threshold is represented in Fig. 7.15 : whatever falls on the left hand side will be assimilated to noise, everything falling on the right hand side will be detected as signal. With this value, we compute the proportion of undetected signals.

GB sources separation

Going toward more realistic data -such as the galaxy signal proposed during the LISA Data Challenges 1 and 2 [20], we wish to have a model that is adapted to the case where the measurement s is the sum of several overlapping signals. When several sources are lying in the same frequency range of interest, their identification becomes more complicated as they get closer in frequency, because their signals mutually contaminate each other by overlapping. This phenomenon is represented in Fig. 7.19 : when only the sum of two signals is measured, separating the individual signals becomes challenging. Even estimating the precise number of sources is difficult.

To face this problem, one can adopt a simple pursuit strategy, i.e. estimating at each step which signal seems to correlate the most with the residuals of the previous steps. This strategy can work if the signals do not overlap too much or if overlapping signals present very different characteristics (for instance, if one bright source is overlapping with a low SNR source, it might be possible to first identify the bright source and after subtracting the signal to identify the other one without any loss). However, there is generally no guarantee that the signals will be properly estimated without any leakage or over-fitting others overlapping signals. The presence of other signals is likely to bias the identification of one source.

A more appropriate strategy is to conduct a joint search of the signals : this way, the contamination that happens when signals are overlapping is directly taken into account during the search. Such a strategy can be conducted using an IAE-based algorithm, in a fashion similar to the barycentric span projection algorithm that we described earlier.

Indeed, the search Eq. (BSP) can be generalized for the case where K > 1 sources are lying in the same frequency range of interest (i.e. the frequency range serving as input for IAE). To the k-th signal, we can associate its barycentric coordinates {λ (k) i } i , its amplitude A (k) and a shift τ k -for now and on, signals are not centered anymore. Then we can identify the characteristics of each signal solving the joint problem :

A(s), Λ(s), T = Argmin A = {A (k) } Λ = {λ (k) i } i,k T = {τ k } k s - K k=1 D τ k A (k) • Ψ {λ (k) i } i 2 2 , (multi-BSP)
where Λ(s) gathers the barycentric coordinates of all the sources identified in s, A(s) is the collection of amplitudes, T is the vector of shifts for each source and D τ is the delay operator applying the shift τ to a source. More details can be found about shift estimation in App.D.6. De facto, Eq.(multi-BSP) is solved using JAX [19] optimizing tool.

In the next paragraph, we demonstrate the performances of source separation with IAE for an observation composed of a superposition of two sources and, elaborating on this example, we investigate a process to identify the number of sources that are present.

Signal estimation

We consider 1, 000 observations of type :

d = h 1 0 (SN R 1 ) • h 1 + h 2 0 (SN R 2 ) • h 2 + n , (7.4.1)
where h 1 , h 2 are two signals with amplitudes h 1 0 , h 2 0 computed such that they reach a chosen SN R with regard to noise n. Their main frequencies are chosen in [2,3] mHz such that the two signals are overlapping. For this benchmark, the SN R is chosen ac-cording to a uniform law : SN R ∼ U([10, 30]) 9 . We investigate here the relevance of using the multi-BSP instead of a classical BSP. The right hand side plot displays the residuals once the IAE estimate has been subtracted. After applying the detection test proposed in Section 7.3.3, they do not present any power excess that could be interpreted as signal leftover or as an undetected signal. Performing a joint estimate of the signals proves to be better than subtracting the signals one by one in a "pursuit" fashion. Indeed, if two GB signals overlap, then estimating the signals one by one can importantly mislead the results. Fig. 7.21 (middle plot) compares the results given by IAE, either combined with the simple BSP or the multi-BSP. For the simple BSP, we display the residual signal after identifying one single source (which corresponds to the first step of a pursuit mechanism) to show that the presence of two signals bias the recovery process. Indeed, the simple BSP blindly identifies a signal that lies exactly between the sum of signals, leaving significant residuals. We can directly compare the quality of recovered signals when applying simple BSP or multi-BSP. The distribution of the quality factor Eq.(5.5.4) is displayed in Fig. 7.22 for the simple BSP and multi BSP. Using the multi-scale IAE model, we see that multi-BSP contributes to improve the quality of the recovered signal in average by a factor 2. On the other hand, the full IAE model estimates signals of equivalent qualities when using simple BSP and multi-BSP. It reaches the same quality distribution as when using the multi-BSP of the multiscale model. Here, we see a clear difference of behavior between the two models : the full model cannot distinguish whether there are one or two signals (since both hypothesis result in average in the same approximation quality) whereas the multi-scale model can. Indeed, in presence of two signals using the multi-BSP produces an approximation of better quality than when using the simple BSP. Based on that, we can build a criterion to estimate the number of sources that are present in a frequency rangei.e. for model selection.

multi-scale model

Full model

Prospective work

There are many points deserving further exploration. We present here several preliminary results that will pave the ground of our future work along the lines of representation learning.

Criterion to estimate the number of sources

This work is based on the results of Section 7.4.2 on the quality of recovered signal for the full model and multi-scale model. The simple-BSP IAE identifies a single source which borrows features from the two true signals and leaves residuals with a manifest power excess at two separated locations reminiscent of the main frequencies of the corresponding GBs. This is a case of source confusion. In fact, this power excess provides us a criterion to identify the number of sources : if hypothesis H k : exactly k sources are present, is true, then for any i < k, we will see a power excess showing that the model cannot fit the sources in presence. In order to check if such a power excess can be a valid selection criterion when applying IAE, we applied the simple BSP to the residuals obtained by subtracting the estimated signal (for both simple and multi BSP) to the observed superposition. Then, we applied the detection described in Section 7.3.3 to check if the residuals are compatible with noise.

On Fig. 7.23, we plotted the distribution of the estimated SN R over the residuals obtained after applying simple BSP and multi-BSP, for a full IAE model and a multi-scale IAE model. In the case of IAE multi-scale model problem, we notice that the residuals obtained after applying the multi-BSP model are globally compatible with noise, whereas this is scarcely the case for simple BSP. Thanks to the low dimensionality of the model, we can clearly discriminate hypothesis H 1 and H 2 using the usual detectability threshold. On the other hand, the full IAE model does not allow to discriminate the two hypothesis since all the residuals -for both simple and multiple BSP -are compatible with noise. Because the dimension of the model is too high, it cannot constrain enough the shape of the sought signal.

multi-scale model

Full model As multi-BSP has a better recovery quality than simple BSP, the residuals better compare to noise.

IAE & fast parameter estimation

The last step towards GB characterization is the identification of the individual sources based on the measurements. We present here preliminary results about partial parameter estimation for GB systems using IAE.

An IAE model learns from signals morphology ; and the morphology is directly (although non-linearly) related to the physical parameters used to generate the signals. This is why we expect that the latent space organization is somehow correlated to the physical parameters, or at least shows sensitivity to some of them. Here we typically think about the parameter β : as mentioned in Section 5.4.3, cos(β) is related to the signal's width in Fourier domain. To a lesser extent, λ also has an impact on the signal's morphology. The low size of the training set does not provide a fine enough discretization of the joint parameter space for (β, λ) to have a more precise approximation. The second element is related to the number of APs and information redundancy : two elements with close shapes can be encoded using two different APs that are highly correlated. This is a priori not the case since we tried to select APs as much as possible, but 26 APs (for the full model) is still too many in comparison with the underlying number of parameters. This can mislead the closest neighbor search.

To palliate this issue, we could either increase the number of elements in the training set, use several neighbors instead of only 1, or change the estimator. For instance, we could learn a multi-layer perceptron (MLP) directly based on the latent domain to perform a regression in the parameter space. We could even imagine that this MLP would be learned simultaneously with the IAE encoder and decoder.

If we were able to estimate quickly intrinsic parameters β, λ, f 0 , we would have a very strong prior to perform parameter estimation : App. A detailed how the sampling can be accelerated when separating intrinsic and extrinsic parameters.

Conclusion

We introduced a way to create a non-parametric model for GB signals which, to the best of our knowledge, has so far not been considered in the field of GW analysis. By non-parametric, we mean that the model does not rely explicitly on the GB system physical parameters Even if the parameters are implicitly contained in the choice of the training set, we do not rely on them to estimate the signal. Instead, we rely on the use of "pseudo-parameters" (A, {λ i } i ).

This model can represent a wide range of signal morphologies and can be used either to compress information related to GB waveforms or to perform the denoising of observed signals. We paved the way for a new pipeline that could be dedicated to GB detection and partial characterization by extending the spectrum of potential applications to detection, partial parameter identification and source separation. On the whole, we expect IAE models to prove useful to LISA data analysis -not only because it is very simple to apply to any type of data, but also because of the wide range of applications it opens.

Prospects : Many improvements can be thought of to enrich this model. First of all, processing jointly all data channels is the first step towards a model that would even better constrain the extraction of sought signals. We can also think about changing the architecture of the encoder Φ and decoder Ψ, for instance by using convolutional layers. These convolutional layers could replace our wavelet-transform preprocessing for a preprocessing that would be even more adapted to the data we want to learn from. Then, there is the question of the frequency range on which the model is learned. For now, it is quite restricted to a range of about 1 mHz. However, the final goal should be to have a model that can represent all types of GB signals.

Finally, there is still a lot more work to do on the potential applications. Many more could come to mind, they would need to be investigated thoroughly. In particular, we think that IAE models bears some resemblance with the reduced basis models [5; 21; 22; 23; 24] introduced in LVC. In our opinion this connection deserves some further investigation.

Chapter 8 Conclusion Major results & prospects

We summed up here the "take home" messages concerning the results of this thesis.

Sparse Signal Representation :

We developed a model based on a very simple hypothesis that enables to get an instantaneous estimate of the underlying GB signals. It works independently of the number of sources. This type of model can be used in cases where we do not know precisely the waveform of the signal we seek, but nevertheless know some of its characteristics that could be used to identify it regardless. There is also a potential application with glitches fast detection and removal. On a wider scale, sparse modeling could be used to operate source separation.

Inpainting LISA data : We proposed two algorithms that can mitigate the impact of gaps on data analysis. In particular, the noise distribution is recovered, and signals are detected properly. These algorithms recreate the missing signals and missing noise, and the processing time is independent of the number of detected sources -all sources are managed simultaneously without having to identify them individually. The algorithms proved to be efficient with as much as 28% of missing data, and was shown to correct posteriors of MCMC identification on an example. Of course, we would like to extend this algorithm to the context where several types of sources are present -in particular transients in time domain. The generalized inpainting algorithm is expected to keep a similar layout, alternating between signal and noise estimation. We also wish to compare the inpainting technique with other existing methods dealing with data gaps.

Learning-based representation :

We introduced a model based on the interpolatory auto-encoder that ensures the low-dimensionality of the underlying representation. Thanks to this property, the denoising process is efficient and can lead to high detection rates. Moreover, the training set is kept to an acceptable size. Many more applications can be thought of, such as template bank compression or source separation. There are many other points to investigate : first, an improvement is expected when processing jointly all information channels. Then the latent space also showed organizational properties that could lead to a fast parameter estimation ; this should be investigated thoroughly. Finally, we could of course develop similar models for other source types or even for glitches.

Final words

There is so much left to do and so little time ! Working on LISA data analysis is a real challenge that twisted my brain for the last three years. Nonetheless, I am glad I was given the opportunity to work on such complex and yet thrilling problems.

Annexe A Parameter Estimation

In this section, we present the asymmetric role intrinsic and extrinsic parameters play in data analysis, and how it can be exploited when performing Bayesian data analysis. We follow the argument of [1].

A.1 Origin of the separation between intrinsic and extrinsic parameters

The distinction made in Tab.4.1 between extrinsic and intrinsic parameters originates directly from the form of the expected waveform for GBs. Indeed, any GB GW signal h I (θ) can be written under the form [2; 1] : Table A.1 -Extrinsic parameters needed to compute h j I (θ int ) in Eq.(A.1.1). Source : [3] Eq.(A.1.1) presents a great interest for parameter estimation, as there is a separation between the impact of the two types of parameters. The idea is now to make use of the asymmetry between intrinsic and extrinsic parameters by first estimating the intrinsic parameters and then deducing the corresponding extrinsic parameters.

h I θ = (θ int , θ ext ) =
Basically, if the log-posterior writes :

log L p (θ = (θ int , θ ext )) = - 1 2 d -

C.2 Linearity of low frequency inpainting algorithm

Let us denote :

F = 1 √ N T exp -2iπkn N T 0≤n≤N T -1 k∈[-K,K] ,
(C.2.1) the matrix containing the coefficients needed to perform the discrete Fourier transform. F is invertible such that :

F -1 = F † . (C.2.2)
For any measurement V ∈ R N T x2 we have :

v[k] = Fv[n] . (C.2.3)
Using the Fourier operator F, we can express the noise correlation matrix in time domain :

C I = F † SF . (C.2.4)
We denote by C the joint time correlation matrix for all the channels. Let v g [n] be a gapped measurement with mask M . We want to study the dependence of : Cette thèse porte sur la caractérisation des systèmes binaires galactiques (GBs) par ondes gravitationnelles dans le cadre de la mission LISA (Laser Space Interferometer Antenna) [1]. LISA est un détecteur spatial d'ondes gravitationnelles dont le lancement est prévu en 2034. Constitué de trois satellites formant des bras interférométriques de 2.5 millions de km de longueur, il sera capable de détecter des sources d'ondes gravitationnelles de fréquences beaucoup plus faibles que les observatoires terrestres. Parmi elles, les scientifiques estiment que les binaires galactiques seront des sources omniprésentes, car notre galaxie en compte plusieurs millions. Toutefois, seulement quelques dizaines de milliers seraient potentiellement identifiables, car ces sources sont de faible amplitude. De plus, LISA ne présente qu'un nombre limité de canaux d'observation et mesurera la somme de ces signaux. Dans ce contexte, les caractériser individuellement est un véritable défi car d'autres sources d'ondes gravitationnelles sont, elles aussi, présentes. De plus, LISA sera sujet à un bruit instrumental complexe et présente de nombreux artéfacts (interruptions de prises de données, glitches) [2].

f CI (v g ) = Argmin q M q[n]=vg[n] 1 2 ⟨q[k], q[k]⟩ S = Argmin q M q[n]=v[n]

D.2 Notation index & definitions

L'analyse des données LISA doit tenir compte des caractéristiques du bruit et des artéfacts et prouver qu'elle est robuste face à eux. C'est dans ce cadre que je propose d'introduire des méthodes non paramétriques d'analyse de données ; ces méthodes, si elles ne permettent pas l'identification directe des systèmes, sont très complémentaires avec l'approche Bayésienne qui est la plus représentée aujourd'hui au sein du consortium. Elles proposent un cadre permettant d'intégrer facilement la présence de bruit et d'artéfacts dans les données.

E.2 Modélisation parcimonieuse des signaux issus de binaires galactiques

Modélisation parcimonieuse

La parcimonie [3; 4; 5] est un cadre qui convient particulièrement bien aux signaux de GB. En effet, les GB émettent des signaux lisses, stationnaires et presque sinusoïdaux. Sur la base de cette observation, la meilleure solution est de procéder à leur étude directement dans le domaine de Fourier où toute l'information se rassemble principalement sur quelques atomes de Fourier. Par conséquent, nous pouvons dire qu'ils sont pour approcher le signal impacte la qualité de l'approximation Q dB . L'évolution de la qualité en fonction du nombre de coefficients utilisés pour l'approximation évolue en deux temps. Au début, la pente de progression est forte : chaque coefficient ajouté contribue fortement à l'amélioration de la qualité du signal reconstruit. Mais au-delà d'un certain nombre de coefficients (la limite se situe autour de 40 coefficients dans cet exemple), la pente se casse et l'amélioration apportée par l'ajout d'un coefficient supplémentaire est beaucoup plus faible. Néanmoins, en présence de bruit, la courbe de récupération est dégradée et la tâche d'estimer le signal malgré le bruit devient un véritable défi.

Application & Résultats

Nous proposons de construire une approximation parcimonieuse d'un jeu de mesures bruité afin de séparer les signaux GB du bruit. Cette séparation s'appuie sur un algorithme de type LASSO [6] pour lequel nous proposons plusieurs versions permettant de prendre en compte plusieurs canaux d'information ainsi qu'une modélisation basique de la morphologie des signaux recherchés. Ce dernier point s'appuie sur la mise en place d'une parcimonie dite "par blocs" [4; 3] permettant de traiter conjointement des clusters d'information. Ces blocs sont adaptés aux mesures via la mise en place d'un arbre (algorithme dit "BlockTree") et contribuent fortement à améliorer la qualité de reconstruction du signal ainsi que les capacités de détection. Nous effectuons par la suite une étude de performances détaillée afin de mieux caractériser le comportement du modèle. Ce travail a conduit à une publication [7].

Nous avons appliqué notre modèle au cas proposé lors du LISA Data Challenge 1 [8] sur l'ensemble des binaires de vérification (challenge 1-3). La mesure était composée de 10 signaux GB entâchés de bruit instrumental complexe. En le domaine temporel ont une phase qui coincide très bien avec celle du signal recherché, mais l'amplitude est moins bien estimée. Ce dernier point pourrait être amélioré en prenant mieux en compte les spécificités des signaux GB.

E.3 Corriger les données manquantes par inpainting

LISA subira des interruptions temporelles de prises de données de fréquences et durées variables [2]. Nous montrons que ces interruptions peuvent impacter la qualité d'identification des GB car elle se fait directement dans le domaine de Fourier pour des raisons de coûts calculatoires.

Prendre en compte ces interruptions dans l'analyse pourrait se faire simplement si l'estimation des paramètres se faisait sur les données temporelles. Mais le coût calculatoire serait alors beaucoup plus élevé. Afin de pouvoir continuer à réaliser l'analyse directement dans le domaine de Fourier même en présence d'interruptions de prises de données, nous proposons une méthode consistant à estimer les données manquantes -à la fois le signal et le bruit -via deux algorithmes d'inpainting. Dans les deux cas, l'estimation jointe du bruit est nécessaire en raison des caractéristiques du bruit instrumental.

Le premier algorithme, appelé Inpainting Classique, est entièrement déterministe et ne peut recouvrir que des bruits basses fréquences dans les périodes d'interruptions. Le deuxième, appelé Inpainting Modifié, et une variante du premier permettant de recouvrir les bonnes distributions de bruit au niveau des trous. La Fig.E.3 représente un bruit typique estimé par le premier algorithme dans un trou. Malgré une estimation très lisse du bruit, la transformée de Fourier du signal ainsi corrigé admet bien un spectre de bruit qui correspond à celui attendu, modulo une perte de puissance (due à l'absence de bruits haute fréquence dans les trous). L'algorithme d'Inpainting Modifié permet de Dans les deux cas, le signal est efficacement recouvert en utilisant le modèle parcimonieux mis au point auparavant. Les algorithmes ont été éprouvés pour des pertes allant jusqu'à 28% des données totales observées. Les résultats de l'étude détaillée des comportements des algorithmes ont été publiés [9]. La On touche à la limite du modèle de représentation parcimonieuse lorsqu'on commence à s'intéresser aux signaux GB individuellement. En effet, l'hypothèse de parcimonie ne permet pas de séparer des signaux de fréquences proches, ni de caractériser de manière précise la morphologie du signal détecté. Nous avons donc développé un nouveau modèle de faible dimensionnalité, toujours non-paramétrique1 et capable de représenter de manière précise la diversité morphologique des signaux issus de GB.

Nous avons choisi d'utiliser un modèle non linéaire en raison de la non-linéarité qui existe entre les paramètres physiques d'un système et le signal correspondant observé. Dans ce contexte, les Auto-Encodeurs sont les modèles non-linéaires construits par apprentissage qui donnent les meilleurs résultats [10]. Parmi tous les modèles existants, Une partie du travail réalisé pour la construction du modèle est liée au choix "optimal" de ces points d'ancrages. Une autre partie se concentre sur les applications possibles, à savoir : générer un modèle approximatif de signaux GB de faible dimension, débruiter des observations, détecter la présence de signaux GB au sein d'une mesure bruitée, mais aussi séparer des sources GB proches et réaliser une estimation rapide partielle des paramètres physiques correspondant au signal. La Fig.E.6 illustre la possibilité d'approcher des signaux GB par des modèles de faible dimension ; on compare avec un modèle linéaire de même dimensionnalité, à savoir une reconstruction basée sur les vecteurs principaux d'une PCA2 apprise sur l'ensemble d'entrainement. Si l'on contraint le modèle à avoir une très faible dimensionnalité, alors le modèle linéaire approche très mal les signaux car ils sont trop complexes. En revanche, le modèle appris arrive à représenter les données de manière précise. De plus, il est possible d'extraire des informations plus précises sur la source qui a émit le signal, comme typiquement sa position dans le ciel.

E.5 Conclusion

Nous résumons ici les messages à retenir concernant les résultats de cette thèse.

Représentation parcimonieuse des signaux : Nous avons développé un modèle basé sur une hypothèse très simple qui permet d'obtenir une estimation instantanée des signaux GB sous-jacents. Il fonctionne indépendamment du nombre de sources. Ce type de modèle peut être utilisé dans les cas où nous ne connaissons pas précisément la forme d'onde du signal que nous recherchons, mais nous connaissons certaines de ses caractéristiques qui pourraient être utilisées pour l'identifier malgré tout. Il existe également une application potentielle avec la détection et la suppression rapide des glitches. A plus grande échelle, la modélisation parcimonieuse pourrait être utilisée pour opérer la séparation des sources.

Corriger les données manquantes par inpainting : Nous avons proposé deux algorithmes qui peuvent atténuer l'impact des lacunes des données sur leur analyse. En particulier, la distribution du bruit est récupérée, et les signaux sont détectés correctement. Ces algorithmes recréent le signal manquant et le bruit manquant, et le temps de traitement est indépendant du nombre de sources détectées : toutes les sources sont gérées simultanément sans avoir à les identifier individuellement. Les algorithmes se sont avérés efficaces avec jusqu'à 28% de données manquantes, et il a été démontré qu'ils corrigeaient les distributions a posteriori de l'identification MCMC sur un exemple. Bien entendu, nous aimerions étendre cet algorithme au contexte où plusieurs types de sources sont présents -en particulier les transitoires dans le domaine temporel. L'algorithme d'inpainting généralisé devrait conserver une disposition similaire, alternant entre l'estimation du signal et du bruit. Nous souhaitons également comparer la méthode d'inpainting avec celle proposée dans [13].

Représentation apprise des signaux GB : Nous avons introduit un modèle basé sur un auto-encodeur interpolatoire qui garantit la faible dimensionnalité de la représentation sous-jacente. Grâce à cette propriété, le processus de débruitage est efficace et peut conduire à des taux de détection élevés. De plus, l'ensemble d'apprentissage est maintenu à une taille acceptable. De nombreuses autres applications peuvent être envisagées, comme la compression de banques de modèles ou la séparation de sources. Il existe de nombreux autres points à étudier : tout d'abord, une amélioration est attendue lors du traitement conjoint de tous les canaux d'information. Ensuite, l'espace latent a également montré des propriétés organisationnelles qui pourraient conduire à une estimation rapide des paramètres ; ceci devrait être étudié de manière approfondie. Enfin, nous pourrions bien sûr développer des modèles similaires pour d'autres types de sources ou même pour les glitches. Grâce à une modélisation simple des signaux recherchés, nous montrons qu'il est possible de détecter les signaux en présence de bruit et de les estimer. Nous expliquons ensuite comment ce modèle peut être utilisé pour atténuer efficacement l'impact des données manquantes sur l'analyse. Enfin, nous étudierons ce qu'un nouveau modèle appris peut apporter en termes de caractérisation du signal.

Title : Characterization of galactic binaries by gravitational waves Keywords : Gravitational waves, Galactic binaries, Sparse modeling, Missing data, Inpainting, Learned representation, Auto-encoder Abstract : The forthcoming space-based gravitational wave observatory LISA will open a new window for the measurement of gravitational waves, making it possible to observe emitting systems hardly visible with the current Earth-based observatories.

Among these sources, the detection of galactic binaries promises an unprecedented wealth of information about these systems, but also raises several challenges in signal processing. In particular the large number of expected sources and the presence of both complex instrumental noise and artifacts tainting the data call for the development of robust methods.

Through simple modeling of the sought signals, we show that it is possible to detect them accurately in presence of instrumental noise and to recover the signals. We then explain how this model can be used to efficiently mitigate the impact of missing data on the analysis. Finally, we investigate what a new learning-based model can bring in terms of signal characterization.

T 2 QmgF

 2 Partial derivative with regard to variable x β I ij (t) Second mass momentum µ Mass density distribution x(t), y(t), z(t) System coordinates with time R, M, Ω System radial distance to observer, system mass and system orbital frequency L(t) = L * + δL Distance between two test masses, observing a variation δL y 1 (t), y 2 (t) Interferometric measurements L 1 , L 2 Michelson interferometer arm lengths n 1 (t), n 2 (t) Noises other than laser noise for each arm h 1 (t), h 2 (t) Gravitational wave strain measured for each arm c Speed of light in vacuum, equals to 299.792.458 m.s -1 X, Y, Z Pseudo-Michelson interferometric measurements, reconstructed by TDI. Data in time domain, data in frequency domain d, d I Measurements for all channels / for channel I h, h I Signal for all channels / for channel I h, h I Signal estimator for all channels / for channel I n, n I Noise for all channels / for channel I S, S IJ Cross power spectral density, Cross power spectral density between channels I and J in frequency domain E[X] Mathematical expectation of random variable X θ,Θ GB Waveform parameters , Set of admissible waveform parameters θ = (θ int , θ ext ) Separation between extrinsic and intrinsic parameters h(θ), h I (θ) Signal for all channels / for channel I corresponding to parameter θ f 0 , ḟ , β, λ, A, ι, ψ, ϕ 0 Parameters for a GB waveform : frequency, frequency derivative, ecliptic latitude, ecliptic longitude, amplitude, inclination, polarization, initial phase p Dictionary representing the signal ∥α∥ 0 = j 1[α j = 0] Norm 0 of h in dictionary T (counts the number of non-zero components of a signal h when approximated by T) ∥α∥ p = j |α j | p 1/p Norm p of α ∥α∥ 1 = j |α j | Norm 1 of α t k = e 2πikn/N Fourier atom for DFT at frequency f k ⟨v 1 , v 2 ⟩ S Inner product of v 1 and v 2 reweighted by covariance matrix S A S I (v I ), A S (v) Signal amplitude per Fourier atom for channel I / for all channels γ[k] Threshold associated to k-th Fourier atom χ Probability distribution of a χ 2 of order Q ρ Rejection rate x 0 p-value associated to rejection rate ρ γ m [k] Threshold associated to k-th Fourier atom, computed at iteration m ϵ γ Stopping criteria for threshold computation κ Parameter used to compute γ m which values impacts the reweighting process B = B i i Fourier domain decomposition over disjoint frequency ranges B i A B S (v) Signal amplitude per neighborhood of Fourier atoms B ∈ B γ B [B], γ m B [B] Threshold associated to Fourier neighborhood B, threshold associated to Fourier neighborhood B, computed at iteration m ϕ(t) Phase of GB signal with time ρ tree Rejection rate associated to BlockTree algorithm ρ reweighting Rejection rate associated to reweighting process E 2 (h, h) Quadratic relative error made when estimating h by h Q dB (h, h) Quality factor (also known as Normalised Mean Square Error) Missing data estimator J (v, u) Cost function associated with classical inpainting algorithm d p+1 Missing data iterative solution, converging to d m h p+1 Signal iterative solution, converging to h n Noise estimator u Dummy variable for missing data. u ∈ Ker(M ) q Dummy variable for noise n p Gapped estimated noise at iteration p n p+1 Noise iterative solution, converging to n d p Updated full data at iteration p f CI (v g ) Classical inpainting operator (inpaints gapped data v g ) D N Difference between gapped and inpainted data n samp Random noise sample drawn according to the noise expected distribution n p+1 gap Gapped noise at iteration p augmented with gapped sample M n samp L gap , T gap Gap length and period f c Cut-off frequency of operator f CI η I Whitened noise D KL [k] Kullback-Leibler divergence at frequency k Matrix of Fourier coefficients C I Noise correlation matrix in time domain L v I (q I , Λ) Lagrangian function associated to the problem prox τ G Proximal operator of function G D, D Primal domain, dual domain prox D G (u) Proximal operator for pre-conditioned algorithm wavelet transform SN R est Estimated SNR on the recovered waveform µ est , σ est Mean and standard deviation of estimated SN R distribution produced by IAE for noise-only inputs H 0 , H 1 Hypothesis tested with p-value test SN R lim Detection limit based on estimated SN R A(s), Λ(s), T Multi-BSP amplitudes, barycentric coordinates and delays H k Hypothesis associated to the presence of k sources Λ train Set of barycentric coordinates for all training set elements {λ test i } i Barycentric coordinates of an element of the test set Chapter 1

Figure 2 . 1 -

 21 Figure 2.1 -Left : LIGO interferometers detection principle. This is based on Michelson interferometer.(source :[13]) Right : Test mass motion in spacetime in absence of GWs corresponding to a flat spacetime (orange, labelled as "Flat") and when a GW is passing (blue, labelled as "GW")

Figure 3 . 1 -

 31 Figure 3.1 -GW sources sorted according to their expected emission frequencies. Different types of interferometers are sensitive to different frequency ranges, and thus will observe different phenomena -with potential overlaps. The "Merging white dwarfs in our Galaxy" represent most part of GB signals that will be observed by LISA. (Source : https://www.esa.int/ESA_Multimedia/Images/2021/09/The_ spectrum_of_gravitational_waves)

Figure 3 . 2 -

 32 Figure 3.2 -Left : The LISA constellation follows the Earth on its orbit around the Sun. It forms a 60 degrees angle with the ecliptic plane. (Source : [1]) Right : One of the three spacecrafts of the LISA constellation. Two laser beams are pointing toward the other spacecrafts. (Source :[8])

Figure 3 . 3 -

 33 Figure 3.3 -Simplified diagram of the two optical benches aboard a spacecraft and the interferometric measurements s 1 , s ′ 1 , τ 1 , τ ′ 1 that are done. Each bench has its own laser (p 1 , p ′ 1 ), making compulsory to compare the inputs and outputs on the two benches. The optical benches movement ⃗ ∆ 1 , ⃗ ∆ ′

Figure 3 . 4 -

 34 Figure 3.4 -GW are transverse waves admitting two types of polarization : h + and h × , forming an angle of π/4 Rad in the plane orthogonal to the direction of propagation. Were LISA constellation in that plane, the inter-satellite distance would be affected by the GW in a continuous deformation depending on the wave polarization. The distance variation measured on each arm will be different. (Source : [2])

Figure 3 . 5 -

 35 Figure 3.5 -Michelson interferometer with two different armlengths L 1 , L 2 and two interferometric measurements y 1 , y 2 . The laser is stabilized around frequency ν 0 but has small fluctuations C(t). In this case, we can find a combination of the delayed measurements depending on the armlengths that cancels the laser fluctuations. (Source : [9])

Figure 3 . 6 -

 36 Figure 3.6 -LISA noise origin. Frequency dependence and origin of the noise in Fourier domain. (Source : [13])

Figure 3 . 7 -

 37 Figure 3.7 -LISA sensitivity strain curve with Fourier frequency. Over this sensitivity curve are represented the main GW sources that should be detected by LISA. (source : [1])

4. 1 . 1

 11 Single GB source

Figure 4 . 1 -

 41 Figure 4.1 -Frame changes between the Solar System Barycenter frame (left) and the source frame (right). This involves a rotation of angle ψ illustrated by the picture in the middle. (Source :[1])

Figure 4 . 2 -

 42 Figure 4.2 -GB waveform with central frequency f 0 = 3 mHz. Left column : modulus of the signal in Fourier domain ; Right column : signal in time domain. Top : full signal ; Bottom : zooms.In time domain, we observe that the "high frequency" signal, which is the one that we seek, admits a very low frequency envelop that is due to LISA constellation movement around the Sun. In Fourier domain, this Doppler effect causes the signal to spread over neighboring Fourier frequencies.

Figure 4 . 3 -

 43 Figure 4.3 -In Fourier domain, simulation of realistic data containing the GB signals emitted within a galaxy (red), MBHB signals (green) and instrumental noise (corresponding noise PSD in black). The total signal is plotted in blue. The galaxy is made of millions of GB systems, all emitting GWs in a restricted frequency band. Some signals are brighter than the others, which will facilitate their identification. Others will be assimilated to a continuous background noise, the "confusion noise". (Source :[4], challenge 2a)

( 4 . 3 . 6 )

 436 Equivalently, we can directly write an estimator for the parameters associated to the signal h just by introducing a slight change in Eq. (4.3.6) :

Figure 4 . 4 -

 44 Figure 4.4 -Plot realised using [7]. Sampling realised using [8; 9]. Behavior of the log-posterior L p (θ)

Figure 5 . 1 -

 51 Figure 5.1 -Left : A random picture of my mom's cat. (credit : A. Blelly ) Right : Sorted normalized amplitudes of the Bi-orthogonal wavelet transform atoms at chosen scale. They decrease like a power law of the coefficient index.

Figure 5 . 2 -

 52 Figure 5.2 -Left : For a fixed threshold γ (dashed red line), we approximate a waveform (in blue) based on the largest 34 (upper row) and 80 (lower row) components. The corresponding approximations are plotted in orange.(To increase the number of coefficients used to approximate the signal, we increased the input signal's amplitude while keeping the threshold at constant value). Middle : The number of components used for the approximation will mainly impact the estimation of the signal time envelop. The approximation gets better with increasing number of atoms considered for the approximation. Right : Zooming over a short time period, both estimates (with 34 and 80 components) are in phase with the

Figure 5 . 3 -

 53 Figure 5.3 -Quality Q dB of GB signal approximation with number of Fourier atoms h 0 used to approximate the signal. Two regimes are observed : the first 40 atoms are key to obtain a trustworthy approximation. Considering more atoms then still improves the approximation, but with a lesser impact.

Figure 5 . 4 -

 54 Figure 5.4 -Soft thresholding applied to a sinusoidal signal. The true signal is in light orange, the noisy signal is in blue and the corresponding sparse solution is in dashed green. The constant threshold γ is represented in dashed red. Left : Signal in time domain ; Right : Amplitude in Fourier domain. Soft thresholding consists in setting to 0 any coefficient below the threshold γ ; the coefficients above threshold are kept up to a coefficient depending on the threshold. The amplitude bias is directly related to the choice of threshold γ, and has consequences on the time series amplitude.

Figure 5 . 5 -

 55 Figure 5.5 -GW signals on channel E as a function of time : target noiseless signal h E [n] (green) and sparse recovery h E [n] (red). Upper row : signal estimation without reweighting process. The threshold determination causes an underestimation of the signal. Lower row : signal estimation with reweighting. The bias is mostly corrected. (Figure extracted from Ref.[10])

Figure 5 . 6 -

 56 Figure 5.6 -Joint signal amplitude A S (d)[k] as a function of Fourier atom frequency (solid blue) and compared to the threshold γ (dashed orange).Upper row : When the signal is present, there is a power excess over all Fourier atoms with frequency close to f 0 , forming a peak. However, within the peak, some

Fig. 5 .

 5 Fig.5.7 portrays this diversity through two examples where the signals are obtained using exactly the same physical parameters except for the ecliptic latitude β. On the upper figure, choosing β = 0 means that the source is in the ecliptic plane : this setting maximizes the Doppler broadening of the frequency peak. On the other hand, on the lower figure choosing β = π/2 corresponds to a source that is at local noon with regards to the ecliptic plane. This setting minimizes the Doppler broadening of the frequency peak. Both peaks are presented with the parameters indicated in Appendix B.2.

Figure 5 . 7 -

 57 Figure 5.7 -For an observation duration long enough, the peak width is proportional to | cos(β)|, where β is the ecliptic latitude for the observed GB. Thus, the largest peak is obtained for β = 0 (top plot), whereas the thinnest peak is obtained for β = π/2 (bottom plot). (Figure extracted from Ref.[10])

Figure 5 . 8 -

 58 Figure 5.8 -Illustration of the BlockTree algorithm. We use the notations of Algorithm 3. Evolution of the domain decomposition in a simple case. The red blocks represent inactive blocks, i.e. those which signal power do not exceed the block threshold. The green blocks represent active blocks, i.e. those which signal power exceeds the block threshold. The top row represents the initial block decomposition, where all blocks have the same size n B . There are a signal peak in the middle and a FP on the right. This FP will be discarded during the first iteration (middle row). The bottom row represents the final block decomposition. The two blocks on the left side are not merged because their sizes are too different (their ratio is larger than R comp ). (Figure extracted from Ref.[10])

Figure 5 . 9 -

 59 Figure5.9 -Average Q dB as a function of the chosen rejection rate ρ for a single binary detection for a single channel resolution ("mono" in blue, corresponding to problem (L 1 , I)) and joint channels resolution ("multi", in orange, corresponding to problem (Joint L 1 )). The averaging is done over the results obtained for a set of 1000 waveforms with parameters randomly sampled and SN R = 50.

Figure 5 . 10 -

 510 Figure 5.10 -Evolution of the FP rate as a function of different initial thresholds for different resolution classes : unstructured sparsity (blue), uniform blocks resolution with a block size |B| = 10 (orange), adapted blocks BlockTree resolution based on the previous uniform decomposition (green). The probability threshold is chosen so that under H 0 , P(χ 2 > x 0 ) = ρ. Left : For structured block sparsity. Right : For unstructured block sparsity. The expected value is given by the first bisector, which is perfectly matched by the two first methods. The BlockTree enables to have an effective FP rate that is much lower than the expected value due to the averaging process. Moreover, unstructured block sparsity has lower FP rate than structured block sparsity.

Figure 5 . 11 -

 511 Figure 5.11 -Evolution of average Q dB with SN R for the 3 main methods presented here : resolutionatom-by-atom ("Unstructured", blue triangles), with uniform domain decomposition ("Uniform blocks", orange circles) and with the BlockTree domain decomposition ("BlockTree ", green squares). For the two block decompositions, the reweighting is operated by block (left) and atom-by-atom (right), for different sets of parameters. The "unstructured" solution is identical on both sides. Up : Unstructured sparsity : ρ = 10 -4 ; Block decompositions : ρ unif orm = ρ tree = 10 -3 . For atom reweighting, ρ reweighting = 10 -1 . Down : Unstructured sparsity : ρ = 10 -6 ; Block decompositions : ρ unif orm = 10 -4 ,ρ tree = 10 -3 . For atom reweighting, ρ reweighting = 10 -1 . The "unstructured" solution has a high rate of both FN and FP.

Figure 5 . 12 -

 512 Figure 5.12 -Average evolution of FN rate with SN R. for unstructured sparsity (Joint L 1 ) (blue,"Unstructured") and BlockTree algorithm (L B 1,2 ) (orange,"Blocktree"). BlockTree has a better detection rate than Unstructured sparsity.

Figure 5 . 13 -

 513 Figure 5.13 -Average quality Q dB of signal recovered using BlockTree for rejection rate ρ tree = 10 -4 . Input signals have a SN R of 50. There is an operating plateau between |B| = 4 and |B| = 20. The area filled in light blue shows how 95% of the values are distributed around their average (dark blue line).

Figure 5 . 14 -

 514 Figure 5.14 -LDC1-3 : Verification galactic white-dwarf binaries : 10 verification binaries have to be found in this signal which noise was generated with the LISACode noise generator. The target signal h is the sum of all GB signals.

Figure 5 . 15 -

 515 Figure 5.15 -In blue : standard deviation σ M AD computed with the MAD estimator and a sampling window length of 5000 points. In red, polynomial fit of order 3 of σ M AD , referred to as σ pol . (Figure extracted from Ref.[10])

  imaginary parts of signals d I , I ∈ {A, E}. The MAD estimate of the standard deviation for ℜ(d E [k]) in presented in Fig.5.15.

Figure 5 . 16 -

 516 Figure 5.16 -Impact of the PSD correction on the noise distribution. The random variable |A S (d)[k]| 2 for a frequency k is expected to follow a law of χ 2 4 as explained in Section 5.3.3. On the left, the signal is normalized with the theoretical PSD PSD th , whereas it is normalized with the effective PSD PSD ef f on the right, for f < 0.027Hz. We computed the noise distribution directly on the signal as the amount of signal is negligible compared to the amount of noise. (Figure extracted from Ref.[10])

Figure 5 . 17 -

 517 Figure 5.17 -Results for the LDC challenge with ρ = 10 -6 and a minimal block size |B| min = 10. For this rejection rate, no FP was found. The global solutions are displayed in the upper row, and a zoom was done on the lower row. Left : solutions in Fourier domain, with a zoom on the 10 th peak (from the left). Right : solutions in time domain, with a zoom on the solutions. The phases match quite well, the agreement between amplitudes is less satisfactory.

Figure 5 . 18 -

 518 Figure 5.18 -In green is represented the distribution of whitened residuals |A S (r)[k]| 2 in Fourier domain. In red is plotted the probability density of the law χ 2 4 . The residuals manifestly follow the expected law. (Figure extracted from Ref.[10])

Figure 6 . 1 -

 61 Figure 6.1 -Fourier transform of simulated noise for an ungapped signal (in orange) and a gapped signal (in blue) for daily short unplanned gaps and maintenance gaps (see simulation details in Sec. 6.3.1). The presence of gaps in time domain leads to a modification of the expected noise distribution in Fourier domain. This modification impacts the mHz frequency band of interest. (Figure extracted from Ref.[4])

Figure 6 . 2 -

 62 Figure 6.2 -Simulated signal received from a GB in Fourier domain : in orange the real, ungapped signal, and in blue the gapped signal (small and planned gaps). Up : Full data. Down : Zoom around the theoretical GB signal. The gaps in time domain lead to a deformation of the signal in Fourier domain : the gapped signal is 'less sparse' than the ungapped signal. The more information we lose, the greater the deformation of the signal.

Figure 6 . 3 -

 63 Figure 6.3 -Fourier transform of LDC1-3 data, composed of 10 verification binaries, for the real ungapped signal (orange) and the signal that we gapped with daily short unplanned gaps and maintenance gaps (blue). Up : Full data. Down : Zoom around the frequency range where the sought signals lie. Signals that could have been identified without gaps are now completely drowned in noise.

Figure 6 . 4 -

 64 Figure 6.4 -Gaps apodization : instead of considering a crenel mask (blue), it is often considered a better idea to smooth the crenel, i.e. to use the orange mask instead, using for instance a cosine function.

Figure 6 . 5 -

 65 Figure 6.5 -Impact of apodization on noise distribution. Upper row : For crenel gaps ; Lower row : For smooth gaps. When using smooth gaps, the noise distribution seems to have the right PSD on higher frequencies (f > 2 mHz).

Figure 6 . 6 -

 66 Figure 6.6 -Impact of apodization on waveform : an ungapped GB signal (orange) is superposed on a crenel-gapped signal (up, blue) and on a smoothlygapped signal (down, green). The gapped signals are more spread than the ungapped signal. Around the emission frequency, smooth or crenel gaps have similar impact on the signal deformation : apodization does not correct waveforms.

( 6 . 1 . 1 )

 611 Let us introduce the "mask" time function m :m[n] = 0 if the data is missing at time nδt ,

.1. 4 )

 4 Their Fourier transforms will be respectively noted d g [k] and d m [k], as stated in Section 4.2.1. Thus, we have the relation between variables : d = d g + d m . (6.1.5)

  ) and M u[k] will refer to the elements of the Fourier transform of M u[n]. As M is a simple projector in time domain, we will use the following algebra notations : Ker(M ) : kernel of M (data vectors u such that M u = 0), Ran(M ) : range of M (image of data vectors in time domain by M ), I d : identity operator. By definition, we have d g ∈ Ran(M ) and d m ∈ Ker(M ) = Ran(I d -M ).

Figure 6 . 8 -

 68 Figure 6.8 -Impact of gaps on the Fourier transform of a discrete sine with finite horizon depending on the amount of lost data and its distribution. The mask m applied to the monochromatic signal is displayed on the left in time domain. On the right, the corresponding ungapped and gapped data are represented in Fourier domain. Upper row : single small gap equivalent to a 5% data loss. Middle row : single big gap equivalent to a 25% data loss. Bottom row : randomly placed small gaps equivalent to a 25% data loss. (Figure extracted from Ref.[4])

Figure 6 . 9 -

 69 Figure 6.9 -Upper row : Inpainted noise in time domain in a gap for classical inpainting algorithm (C.I.).We only recover noise frequencies for which the correlation length is superior to the gap size. This is not a shortcoming, as Fig.6.1 shows that the PSD deformation is more important for low frequencies than high frequencies. Lower row : Comparison between gapped signal and inpainted signal in frequency domain : the noise spectrum was corrected from the imprint of the gaps. (Figure extracted from Ref.[4])

  Figure 6.9 -Upper row : Inpainted noise in time domain in a gap for classical inpainting algorithm (C.I.).We only recover noise frequencies for which the correlation length is superior to the gap size. This is not a shortcoming, as Fig.6.1 shows that the PSD deformation is more important for low frequencies than high frequencies. Lower row : Comparison between gapped signal and inpainted signal in frequency domain : the noise spectrum was corrected from the imprint of the gaps. (Figure extracted from Ref.[4])

Figure 6 . 10 -

 610 Figure 6.10 -Upper row : Inpainted noise in time domain in a gap for modified inpainting algorithm (M.I.). The red line shows the position of the gap. Lower row : Comparison between gapped signal and inpainted signal in frequency domain : the noise spectrum was corrected from the imprint of the gaps. (Figure extracted from Ref.[4])

Figure 6 . 11 -

 611 Figure 6.11 -Distribution over 50 samples n samp of real and imaginary parts of S -1/2 f CI (M (n+n samp ))n sample -n , with n the input noise (before gaps) and f CI (M (n+n samp ))-n samp the recovered inpainted noise for M.I. algorithm and sample n samp . Left : Small gaps. Middle : Medium gaps. Right : Large gaps. The algorithm acts like a low-pass filter of order 1. The sample that we use in order to fill the gaps only impacts frequencies higher than the cut-off frequency f c (represented as a dashed black line), which is directly linked to the gap size. The cut-off frequency depends on the gap duration L gap (in seconds) as : f c = 1/(2L gap ). (Figure extracted from Ref.[4])

Figure 6 . 12 -

 612 Figure 6.12 -Evolution of Kullback-Leibler divergence versus frequency for reference ungapped signal (blue), gapped signal (orange), signal inpainted with C.I. (green) and signal inpainted with M.I. (red) for different types of gaps. Left : Small gaps. Middle : Medium gaps. Right : Large gaps. The power loss becomes more important as the gap size increases, which explains the incompatibility between reference noise statistic and C.I. noise statistic. Gaps impact noise distribution in the mHz band in a non-negligible way. C.I. is effective for small-and medium-sized gaps, but less on large gaps as it cannot make up for the power loss. (Figure extracted from Ref.[4])

Figure 6 . 13 -

 613 Figure 6.13 -FP rate with rejection rate for real ungapped signal (blue), gapped signal (red), signal inpainted with C.I. (orange) and signal inpainted with M.I. (green) for different types of gaps. Left : Small gaps. Middle : Medium gaps. Right : Large gaps. The FP rate for algorithms C.I. and M.I. is identical to the one obtained for an ungapped signal.

  112

Figure 6 . 14 -

 614 Figure 6.14 -Median Q dB with input SN R for different types of gaps. Left : Small gaps. Middle : Medium gaps.Right : Large gaps. This was realized over 10 samples for each point. We compare the recovery quality in absence of gaps (green, labelled as "REF") and in presence of gaps using classical inpainting (blue, labelled as "C.I.") and modified inpainting (orange, labelled as "M.I."). The envelops of corresponding colors show the repartition of all samples. In presence of small or medium gaps, the final Q dB is not

Figure 6 . 15 -

 615 Figure 6.15 -Median Q dB with input SN R for different types of gaps. Left : Small gaps. Middle : Medium gaps. Right : Large gaps. The final Q dB is not affected by small and medium-sized gaps. However, larger gaps deteriorate the quality of extracted signal. C.I. performs better than M.I. as it virtually increases the SN R by not making up for the power loss in the gaps. On the central plot, the outlier observed for M.I. is due to a FP detection.

Figure 7 . 2 -

 72 Figure 7.2 -Scheme of a basic auto-encoder : the data s is encoded to a low-dimensional latent space using the encoder Φ.The low-dimensionality is ensured by the architecture which brings large inputs to small-sized encrypted data c. This specific network design is widely known as the "bottleneck" architecture. The encoded information c can then be decoded using the decoder Ψ to recover an approximation s of the input.

Figure 7 . 3 -

 73 Figure 7.3 -Training set waveform example. Upper row : typical element contained in the training set. Lower row : Logarithm of its module. The waveform module can take values with several orders of magnitude.

Figure 7 . 4 -

 74 Figure 7.4 -IAE principle : it is impossible to describe any signal s in the direct domain as a linear combination of the anchor points {e i } because of data complexity. However, we can learn a non-linear transform Φ such that, in the transformed domain, Φ(s) is now a linear combination of the transformed anchor points {Φ(e i )} i .

Figure 7 . 5 -

 75 Figure 7.5 -Global learning process : elements x of the training set T train are encoded with Φ and then interpolated on the transformed anchor points {Φ(e i )} i using interpolator I. Then the interpolated result is decoded using Ψ, and the obtained output is compared to the corresponding input x.

Figure 7 . 6 -

 76 Figure 7.6 -We trained a IAE network based on 50 APs on the training set. The evolution of the global error Eq.(7.2.2) of the two interpolators Eq.(7.2.4) ("Direct projection",blue) and Eq.(7.2.5) ("Sparse", orange) are plotted against the number of learning iterations. The interpolator based on direct projection has some numerical issues which prevent it from reaching a low reconstruction error. On the other hand, the regularized interpolator (Sparse) converges fast to a low error (between 2 and 3 orders of magnitude below the one obtained by direct projection).

Figure 7 . 7 -

 77 Figure 7.7 -Distribution of normalized module of barycentric coordinates on the training set for several anchor points. The text above each figure gives the proportion of elements that are 0-valued. Left : Anchor point widely used by training set elements : most of the time, the corresponding barycentric coordinate has a high, non-zero value. Middle : Anchor point often used by training set elements. In 20% of cases, due to sparsity constraint the anchor point barycentric coordinate is 0. However, in 80% of cases, this value is non-zero, making it essential to describe well the data. Right : Anchor point never used by training set elements : the corresponding barycentric coordinate is always 0. We can remove this anchor point as it is not used to represent the data.

2 . 2 )

 22 with a fixed number of APs. The elimination stage consists in assessing the usage rate of every APs on the whole training set at fixed Φ and Ψ and eliminate the least used ones. Fig.7.7 shows how we process with this selection : plotting the normalized module of barycentric coordinates for each AP on the training set, we assess its extent of use. Then we rule out the ones that are the less used, i.e. the ones which take the most often values that are close to 0. The APs thus selected are shown in App.D.5. They present a wide morphological diversity.

Figure 7 . 8 -

 78 Figure 7.8 -Wavelet transform based preprocessing.We decompose the signal over several scales and we learn one model per scale. APs are the same for the three scales. We recombine the estimates over the different scales and estimate the amplitude of each scale so that the residuals are as small as possible.

Figure 7 . 9 -

 79 Figure7.9 -Example of waveform approximation of a real signal (dotted blue) by PCA (orange) and IAE (green). The upper row shows one of the signals present in the test set as such, and the lower row represents the logarithm of its module. In this case, IAE performs much better than PCA as PCA could not capture all of the variations of the waveform.

Figure 7 . 10 -

 710 Figure 7.10 -Distribution of the quality estimator Q dB when the waveforms of the training set (left) and test set (right) are approximated using IAE (green), PCA (blue) and projection on elements of the training set (orange). All methods generalized well (errors are consistent between training set and test set), but IAE gives waveform approximations 10 times better (in terms of recovery quality) than PCA for the chosen number of components.

Figure 7 . 11 -

 711 Figure 7.11 -For a noisy input (grey) corresponding to the dotted blue signal at SN R = 50, example of waveform denoising by PCA (orange) and IAE (green). Upper row : comparison in linear scale. Bottom row : module in log scale. For the same model dimension, the IAE recovers all the variations of the waveform whereas PCA can only capture part of the waveform morphology features.

8 .Figure 7 . 12 -

 8712 Figure 7.12 -Distribution of quality estimator Q dB for various input SN R for IAE (green), PCA (blue) and APs (orange). Left : SNR = 20, Middle : SNR = 40, Right : SNR = 100. IAE performance improves with the SN R at variance with the two other models.

Figure 7 . 13 -

 713 Figure 7.13 -Quality factor of the recovered signal for sparse signal modeling combined with BlockTree

Figure 7 . 14 -

 714 Figure 7.14 -Comparison between an input consisting in noise only (light grey) and the signal estimate that IAE performs in this case (green) in linear scale (upper row) (scale on which IAE works) and log scale (lower row) of the module. IAE estimates a signal which amplitude is close to noise level.

Figure 7 . 15 -

 715 Figure 7.15 -Estimated SNR on IAE output for an input consisting in noise only (blue), in signal with SNR = 5 (orange) and in signal with SNR = 10 (green). The dashed red line shows the limit of detectability set like in Eq.(7.3.7). Left : multi-scale model ; Right : Full model. Both models give coherent results when assessing the detectability threshold.

Figure 7 . 16 -

 716 Figure 7.16 -For an input made of noise only, the IAE-estimated SNR matches with a Gaussian distribution.

Fig. 7 .

 7 Fig.7.17 displays the evolution of the estimated FNR with the input SN R compared to the performances of the unstructured sparse representation combined with the BlockTree algorithm realised in Section 5.5.6. For SN R > 7, all signals are detected. When the SN R decreases to 0, the detection becomes harder as the signals estimated SN R are more likely to be compatible with noise. Yet we see that IAE is able to detect signals with very low SN R -down to SN R = 2, even if their recovery has a poor quality. Moreover, we notice a huge improvement when comparing the performances of the new algorithm to the detection performed when using sparse modeling. Indeed, signals are detected by IAE at a much lower SN R than using sparse modeling. The corresponding false negative rate decreases much faster with SN R : at SN R = 7, nearly all signals are detected by IAE, whereas less than 10% are detected by sparse modeling.

Figure 7 . 18 -

 718 Figure 7.18 -Left : Module of the whitened signal with Fourier frequency index. There is 1 GB signal that we want to detect. Right : Estimated SN R of signal found vs sliding window starting index. The SN R is maximal when the signal is centered in the input window.

Figure 7 .

 7 Figure 7.19 -Left : Example of overlapping GB signals in linear scale (upper row) and the log of its module (lower row). Right : In practice, we only observe the sum of the two signals (upper row in linear scale and lower row in log scale). This makes identification complicated, as signals are contaminating each others. Even estimating the exact number of sources can become challenging.

Figure 7 . 20 -

 720 Figure 7.20 -Result of source separation using IAE for a sum of overlapping signals such that SN R 1 = 44 and SN R 2 = 32.

Fig. 7 .

 7 Fig.7.20 (left) shows the type of input we are dealing with : the two signals can be so close that they nearly totally overlap, and have different SNRs. On the middle plot, the noisy input is compared to the noiseless sum of signals and the recovery done by IAE. The IAE estimate is close to the superposition of signals, even in the presence of noise.The right hand side plot displays the residuals once the IAE estimate has been subtracted. After applying the detection test proposed in Section 7.3.3, they do not present any power excess that could be interpreted as signal leftover or as an undetected signal. Performing a joint estimate of the signals proves to be better than subtracting the signals one by one in a "pursuit" fashion. Indeed, if two GB signals overlap, then estimating the signals one by one can importantly mislead the results. Fig.7.21 (middle plot) compares the results given by IAE, either combined with the simple BSP or the multi-BSP. For the simple BSP, we display the residual signal after identifying one single source (which corresponds to the first step of a pursuit mechanism) to show that the presence of two signals bias the recovery process. Indeed, the simple BSP blindly identifies a signal that lies exactly between the sum of signals, leaving significant residuals.

Figure 7 . 21 -

 721 Figure 7.21 -Result of source separation using IAE simple BSP (red) and multi-BSP (green) for a sum of overlapping signals such that SN R 1 = 25 and SN R 2 = 51. Performing a pursuit (i.e. trying to estimate the signals in presence one after another) would result on an erroneous signal estimation.

Figure 7 . 23 -

 723 Figure 7.23 -Applying simple BSP on the residuals after subtracting the signal found by either simple BSP (blue) or multi-BSP (orange), we estimate the SN R of the result to check its compatibility with noise.

Figure 7 .

 7 Figure 7.25 -Fast parameter estimation through nearest neighbor. Quadratic error when using a fast parameter estimation based on closest neighbor. Upper line : without noise. Lower line : With high noise level. Most times, the parameters are properly estimated. Noise degrades the quality of the estimation (the distribution tails are heavier), but the results remains encouraging.

4 j=1a 4 j=1a

 44 j (θ ext ) h I (θ ext = θ 0 ext,j , θ int )h j I (θ int ) = j (θ ext )h j I (θ int ) , (A.1.1)and the θ 0 ext,j are constants defined as in TableA.1.

4 j=1a

 4 j (θ ext )h j I (θ int ) † S -1 d -4 j=1 a j (θ ext )h j I (θ int ) +log [p((θ int , θ ext ))] .(A.1.2)

1

 1 ], q[n]⟩ C . (C.2.5)with the input v g . The lower line is the same expression as the upper line, but in time domain. The resulting cost function is separable with respect to each channel. For a channel I, the Lagrangian of this problem writes solely in time domain as :L v I (q I , Λ) = 1 2 q T I F * S -1 I Fq I + ⟨Λ, v I -M q I ⟩ , (C.2.6)where ⟨, ⟩ denotes the classical hermitian inner product. The solution of the problem can be written using the Lagrangian :f CI (v) = Argmin q I Argmax Λ L v I (q I , Λ) ∂q I (q I , Λ) = F * S -1 I Fq I -M Λ = 0 , ∂L v I ∂Λ (q I , Λ) = v I -M q I = 0 . (C.2.8)They result in :q I = (F * S -1 I F) -1 M Λ , M q I = v I .(C.2.9) Open source code IAE[1] code is open source and can be found online at https://github.com/jbobin/ IAE on version 3 of the GPL (GPLv3).

Norm 1 :Norm 2 :D. 3 D. 3 . 1

 12331 for a signal h in discrete Fourier domain, denoting h = {h[k]} k : for a signal h in discrete Fourier domain, denoting h = {h[k]} k : h[k] * is the complex conjugate of h[k]. Simulation parameter settings Parameter space Our training, validation and test sets are composed of waveforms which parameters are sampled under the statistic distributions described in Tab.D.1. The waveforms were produced using the LDC code [2], with observation parameters set as follow : Total Observation duration : T obs = 31457280 s (for one year of data) Time step : ∆t = 15 s D.3.2 IAE model parameters The sets used for the learning and the tests are composed of : Training set : 5, 000 waveforms Validation set : 1, 000 waveforms

Figure E. 1 -

 1 Figure E.1 -Qualité d'approximation Q dB d'un signal GB en fonction du nombre d'atomes de Fourier h 0 utilisés pour approcher le signal. Deux régimes sont observables : les 40 premiers atomes contribuent fortement à améliorer la qualité de la solution. Au-delà, tout atome supplémentaire considéré ne contribue que faiblement à l'amélioration de la qualité d'approximation.

  Fig.E.2 nous présentons l'estimation de signal que produit notre modèle. Toutes les sources ont été identifiées et aucune "fausse détection" (faux positif) n'a été observée. Les signaux reconstruits dans

Figure E. 2 -

 2 Figure E.2 -Application au LDC 1-3 : extraction des signaux GB à partir de mesures bruitées simulées par modélisation parcimonieuse. Le signal réel est représenté en bleu et le signal estimé par notre modèle est représenté en orange. Les solutions globales sont représentées en ligne supérieure, et des zooms sont représentés sur la ligne inférieure. Gauche : solutions dans le domaine de Fourier, avec un zoom sur le 10ème pic en partant de la gauche. Droite : solutions dans le domaine temporel, avec un zoom sur une période temporelle d'environ 45 minutes. Le modèle a extrait un signal dont la phase correspond exactement à celle du signal réel. En revanche, l'amplitude est moins bien estimée.

Figure E. 3 -

 3 Figure E.3 -Haut : Estimation de bruit réalisée par un de nos algorithmes d'inpainting (C.I.) au niveau d'une interruption temporelle de prises de données. Cet algorithme peut uniquement recouvrir les fréquences de bruit dont la longueur de corrélation est supérieure à la taille du trou. Bas : Comparaison entre l'amplitude du bruit pour un signal troué (bleu) et un signal inpainté (orange). Dans le domaine de Fourier, le spectre du bruit a bien été corrigé de l'empreinte des trous.

Figure E. 4 -

 4 Figure E.4 -Comparaison des distributions a posteriori obtenues par échantillonnage des transformées de Fourier d'un signal complet (vert), d'un signal troué (bleu) totalisant 28% de pertes et du même signal troué corrigé par inpainting (orange). L'estimateur utilisé est le maximum de la postérieure, et doit se rapprocher le plus possible de la valeur vraie (trait rouge). Alors que l'estimation sur des données trouées donne lieu à un biais dans les paramètres estimés, l'inpainting permet de corriger ce biais.

E. 4 Figure E. 5 -

 45 Figure E.5 -Principe de l'IAE : dans le domaine direct, il n'est pas possible de décrire n'importe quel signal s comme une combinaison linéaire des points d'ancrages {e i } car les données sont trop complexes. On peut toutefois apprendre une transformation non-linéaire Φ telle que, pour n'importe quel signal GB s, Φ(s) puisse être exprimé comme une combinaison linéaire des transformations des points d'ancrages {Φ(e i )} i .

Figure E. 6 -

 6 Figure E.6 -Exemple d'approximation d'un signal GB (en pointillés bleus) par le modèle IAE (en vert) et par une PCA à 9 éléments (orange) en échelle linéaire (haut). La ligne basse représente le logarithme du module des signaux réel et approchés. L'IAE donne une meilleure approximation que la PCA car la PCA n'arrive pas bien à capturer toutes les variations du signal.

Titre:

  Caractérisation des systèmes binaires galactiques par ondes gravitationnelles Mots clés : Ondes gravitationnelles, Binaires galactiques, Modélisation parcimonieuse, Données manquantes, Inpainting, Représentations apprises, Auto-encodeur Résumé : Le futur observatoire spatial d'ondes gravitationnelles LISA ouvrira une nouvelle fenêtre pour la mesure des ondes gravitationnelles, permettant d'observer des sources difficilement visibles avec les observatoires terrestres actuels. Parmi ces sources, la détection des binaires galactiques promet une richesse d'informations sans précédent, mais soulève également plusieurs défis d'analyse de données. En particulier, le grand nombre de sources attendues et la présence à la fois de bruit et d'artefacts entachant les données nécessitent le développement de méthodes d'analyse robustes.

Table of Contents IAE

		Interpolatory Auto Encoder
	AP	Anchor Points
	BSP	Barycentric Span Projection
	ISTA	Iterative Shrinkage-Thresholding Algorithm

Notations The two (or more) bodies problem in gravitation : from Newton to gravitational waves

  Gravitational constant, equals to 6.67430 × 10 -11 m 3 .kg -1 .s -2

	-→ F A→B (t)	Newton's gravitational interaction between A and B
	G
	m A , m B	Masses of bodies A and B
	-→ AB(t)	Distance between bodies A and B
	Φ A	Newtonian gravitational potential generated by A
	-→ x , -→ x A (t)	Coordinates of a random vector, coordinates of body A
	m a , m I	Acceleration mass, inertial mass
	G αβ	Einstein's tensor
	T αβ , T αβ	Stress energy tensor
	g αβ , g αβ

Table 4 .1 -Parameters

 4 needed to fully describe a GB without eccentricity. The parameters can be split in two categories : intrinsic parameters and extrinsic parameters, corresponding to their roles in parameter estimation.

  . The most popular are of type Metropolis-Hastings and generically possess the following two stages : Initialization : Set a number of walkers that will probe the search space. To each of them corresponds a Markov Chain. Initialize them randomly over the plausible parameters set. Set a probability distribution g(θ ′ |θ) that suggests a sample θ ′ for the next iteration given that the current sample is θ. Most of the time, g is symmetric : g(θ ′ |θ) = g(θ|θ ′ ).

	Evolution : For each Markov chain :
	■ From the current sample θ t at iteration t, draw a new sample θ ′ according to
	probability g
	■ Compute the acceptance rate α = min	p(h(θ ′ )|d) p(h(θ t )|d)	, 1
	■ Accept or reject : draw u ∼ U([0, 1]), compare α to u :
	-If u ≤ α : accept the new sample :

Table 4 .2 -Priors

 4 p(θ) used for the sampling Fig.4.4.estimate as it was only meant to be used as a simple way to compare our work with what is currently done in the community.

Table 5 . 2

 52

	-3 )

-Peak to peak recovery quality factor Q dB . Peak are numbered from lowest frequency (1) to highest frequency

Table 6 . 1 -

 61 Description of the three types of gaps used for the study. Small and large gaps are both unplanned, whereas medium sized gaps represent the planned interruptions (or maintenance gaps). We simulate one gap of length L gap once every T gap .

	Type	Duration L gap	Period T gap
	Small	10 min	Every 24 hours
	Medium	7 hours	Every 2 weeks
	Large	3 days	Every 12 days

Table 6 . 2

 62

	Peak	SNR	REF	C.I.	M.I.
	1	26.6 9.513	7.566	8.806
	2	22.5 9.081	7.067	7.936
	3	12.0	4.85	3.086	1.799
	4	9.56 1.918	0.0	0.0
	5	65.3 15.773 12.047 13.604
	6	52.7 13.199 10.972 11.456
	7	48.7 13.661 11.180 11.467
	8	19.8 4.927	3.897	2.199
	9	12.3 2.696	0.834	0.743
	10	72.5 13.719 10.677 10.575
	Global		13.108 12.177 10.235

-Peak to peak Q dB for extracted signal from ungapped signal (REF), with classical inpainting (C.I.) and modified inpainting (M.I.). Peaks are ordered from lowest to highest frequency. For this experiment the rejection rate was set to ρ = 10 -6 .

  optimization tools. For the sake of simplicity, we use the notation Ψ ({λ i } i ) instead of Ψ i λ i Φ(e i ) . We do not have a massive training set ; yet having a continuous representation of the signals like the one IAE can generate is beneficial since we simply interpolate the missing waveforms from the one that we already know in the latent domain. Choosing the anchor points among the elements of the training set helps making the model robust despite the low size of T train .

multi-scale model Full model

	0.16					SNR lim Noise only	0.16					SNR lim Noise only
	0.14					SNR = SNR = 10 5	0.14					SNR = SNR = 10 5
	0.12						0.12
	0.10						0.10
	0.08						0.08
	0.06						0.06
	0.04						0.04
	0.02						0.02
	0.00	4	6	8 Estimated SNR 10	12	14	0.00	4	6	8 Estimated SNR 10	12	14

This convention will only be used in this section.

Other types of linearizations are possible, like Post Newtonian developments.

Each new TDI generation is based on hypothesis that are less restrictive than the previous one. More information can be found in[9].

In this simplified case we do not try to estimate the number of sources that are present. Thus we can discard all the terms that are related to the model since we are not performing model selection.

The sin originates from the definition of β as the co-latitude (see Fig.4.1).

The unit production cost is not necessarily high. But MCMC type methods make hundreds of thousands of calls to the waveform, hence the significant computational cost.

Since S is diagonal (refer to noise modeling described in Section 4.2.2), it is easy to compute S 1/2 by simply taking the square root of the diagonal elements of S.

We can solve the equation independently for each Fourier atom.

This is similar to the introduction of a bias through Eq. (5.3.5) and Eq. (5.3.6) as made manifest by the common structure of these thresholding equations.

This value was set after computing the SN R of LDC 1-3 verification binaries. They are given as an indicator in Section 5.6.2

Here we consider the 2•10 6 measurement points in Fourier domain, independently of the frequency range of interest.

The "enough" is directly quantified by the compressibility of x in T, i.e. the number of non-zero coefficients needed to approximate x accurately.

This colored noise is equivalent to having noise time correlation lengths depending on the noise frequency. As the noise is not white, we have to account for its behavior as much as possible.

Private communication from N. Korsakova.

The code works on a single core. It takes about 2 hours to process a 2 years data set with 28% of data loss ; thus experiments requiring repetitive solving of the problem having high computing cost. It could probably be re-written to be more efficient, but I did not have the time to do so.

In further developments, this should be investigated properly. It is all the more possible that investigations tend to show that latent space organization is correlated to the values of physical parameters.

2. It can seem like a harsh approximation. However, we did not notice any additional bias by applying the model learned on real part of A to all other channels.

We use an improved version of ISTA that only performs the 10 first iterations ; the corresponding gradient steps are learned as model parameters to increase the convergence speed.

The matrix is ill-conditioned because the APs are correlated. The more APs are considered, the greater the correlation, the more ill-conditioned the matrix is. This is the origin of numerical instabilities.

We need to retrieve amplitude because we learned a normalized model for each scale, loosing the amplitude information at input.

We chose this interval as it was plausible considering the SN R found for sources of LDC1-3.

In the sense that we consider that the KL divergence is constant over a few consecutive frequency bins

Par cela, nous entendons que la morphologie du signal est modélisée indépendamment des paramètres physiques qui ont servi à générer le signal.

Principal Component Analysis

Acknowledgments

Conferences

Quality factor (dB) : it measures the quality of approximating a signal h by an estimator h. It is measured in decibel based on the reconstruction relative error :

Q dB (h, h) = -10 log 10 E 2 (h, h) .

(5.5.4)

The greater it is, the better the approximation. Q dB gives an idea of the global quality of the reconstructed signal. However, it does not provide any information about the FP and FN rates. This quantity is also known as "Normalised Mean Square Error" in classical signal processing field.

FP rate :

The number and the amplitude of FPs can greatly impact the signal recovery.

In absence of a meaningful signal, the FP rate is defined as the average number of atoms detected as active whereas they should not. For a given rejection rate ρ, we can estimate the FP rate as :

where #FP denotes the number of FPs for the rejection rate ρ.

FN rate :

We can define the FN rate in a similar fashion as the FP rate : this is the average number of atoms detected as inactive whereas they are active (because they carry signal). In practice, we define a FN rate as follow : for a given rejection rate ρ for which the FP rate is low, the FN rate will essentially depend on the input signal SN R (i.e. its power compared to the noise power). It will be obtained as an average over N N oise different noise realisations. For each input signal at chosen SN R and each noise realisation, one will assess whether the signal was detected by the algorithm. Then, the FN rate will be defined as the proportion of experiments (noise realisations) for which the signal was not detected : (5.5.6) where ϵ i = 1 if signal with chosen SNR is undetected when added to i-th noise realisation, 0 otherwise.

Building a representative set

In order to assess the overall performance of our algorithm, it is important to consider the wide variety of waveforms that the algorithm could encounter. To represent properly this diversity, we will benchmark our algorithms on a set of 1000 waveforms sampled through the parameter space. The prior used for this sampling are detailed in Table 5.1. The overall results of the benchmark will be obtained by averaging the results over the randomly chosen waveforms.

The amplitude of the signal can be computed as a function of the wanted SN R : if we produce the signal with an initial amplitude h 0 = 1, we then compute the final amplitude as a function of the SN R given by Eq.(5.5.2) : for any multiplicative factor h 0 and initial signal v, we have :

let us consider the update of d p+1 m :

(6.2.3)

We could remove the sparsity term ∥•∥ 1 as it only depends on h p . Introducing q the noise variable defined by the following change of variables :

one can recast Eq. (6.2.3) as an equation on noise. We define the estimated noise outside the gaps (i.e. observed noise) at iteration p :

This eventually leads to the following equivalent problem :

where n p converges to the noise estimator n for p → +∞.

What is n p+1 ? In Eq.(6.2.6), we are looking for the noise q that matches the observation n p g = d g -M h p outside the gaps and which PSD would match S in Fourier domain. The solution n p+1 is the noise realization that is equal to n p g outside the gaps and that minimizes the residuals (reweighted by the PSD S) in Fourier domain.

Solving this problem is challenging because the noise PSD S is known in Fourier domain while the equality constraint n p g = M q is defined in time domain. This problem does not admit a closed-form expression. To evaluate numerically its minimizer, we use the Chambolle and Pock primal-dual algorithm [16]. This algorithm has two main advantages : i) it remains computationally simple, and ii) it can further be preconditioned to speed up convergence, which is particularly convenient since the noise PSD is ill-conditioned. We refer to Appendix C. 4 for more details about this algorithm and its implementation.

Updating the signal

Let us introduce the full updated data :

which is the data whose gaps have been filled in at the p-th iteration. If d p+1 m is solution of Eq.6.2.3 and n p+1 is solution of Eq.6.2.6, then they are related through the change of variable Eq.6.2.4 by :

f CI (v g ) is the inpainted noise solution corresponding to observed gapped data v g obtained after solving Eq.(6.2.6). In particular, it is equal to v g out of the gaps. In the case of classical inpainting Alg.(1), we had :

Let us draw an ungapped noise random sample n samp that follows the expected noise distribution (see Section 4.2.2, Eq.(4.2.7)). The corresponding gapped measurement is v g = M n samp . Let D N be the difference :

What is D N ? By construction, f CI (M n samp ) has the same value as n samp outside the gaps. Thus D N is a signal term whose value is 0 outside the gaps. Additionally, since f CI (M n samp ) only keeps the low frequency information of n samp (see Section 6.2.2), D N only contains information at high frequency since the low-frequency content has been removed. D N is a high frequency term, defined only where data is missing, and that is compatible with any inpainted noise f CI (M n).

By "compatible" we mean that it could have been the high frequency part of the noise that we could not recover using classical inpainting Eq.(6.2.6). Consequently, it can be virtually added to any gapped measurement that has been inpainted with C.I. without altering the noise PSD.

The resulting modified noise inpainting writes :

where Eq.(6.2.14) exhibits a linearity property with regard to its input, as demonstrated in Appendix C.2.

Modified inpainting (M.I.) algorithm :

The modified algorithm builds upon the C.I. algorithm by adding to the noise solution a stochastic component, like in Eq. (6.2.17). We further prove in Appendix C.3 that the overall process consists in minimizing a cost function that shares similarities with the standard inpainting defined in Eq. (6.2.1), but to which we added a correcting term. It can be solved using the exact same minimization scheme based on the BCD architecture up to a modification of the input and output. More precisely, the resulting algorithm is detailed in Alg. 2.

To further evaluate the impact of the new algorithm on the inpainted noise, we make the exact same test as in Section 6.2.2. We consider simulated noise realizations for LISA (not containing any GW signal), whose statistics is described by the theoretical PSD of the LDC1-3 (see Section 4.2.2). Gapped data are generated with small unplanned gaps and planned gaps (see Section 6.3.1 for more details about gaps generation). Then, inpainting is applied to the resulting gapped noise according to Eq. (6.2.17). Fig. 6.10 shows the solution which has been inpainted with the proposed algorithm. This time the real signal and the inpainted one do not visually show any differences. The power loss was corrected using the added stochastic term. Performances of both algorithms are assessed in details in Section 6.3.

Classical Inpainting

Modified Inpainting Figure 6.16 -Comparison of the posteriors obtained for a full, ungapped signal (green), a crenel-gapped signal (blue) and an inpainted signal (orange). For the gapped and full signal, the posteriors are identical to the one presented in Fig. 6.7. The data loss rate is about 28%, combining the three types of gaps. Left : Using classical inpainting. Right : Using modified inpainting. In both cases, the posteriors of inpainted data are really close to the full data posteriors. Using inpainting corrects the impact of the mask on the signal, and thus helps improving the parameter estimation. Using modified inpainting makes the main frequency f 0 estimation a bit harder, as it also adds noise around that frequency. For the sky location, both inpainting methods provide similar results. of noise distribution, detection capacity, accuracy of extracted signal or even improvement for parameter estimation, to situations with ungapped data. More precisely the C.I. algorithm yields an efficient signal reconstruction but outputs a less satisfactory noise distribution. On the contrary the M.I. algorithm successfully recovers the expected noise distribution but is challenged by the extraction of low SNR signals. We showed on an example that both could help un-biasing the parameter estimation if part of the data is missing. The benefits of the M.I. algorithm mainly lies in the fact that it can recover more precisely the noise distribution, especially in the gaps : whereas the C.I. algorithm provides a smooth solution that correct the overall noise distribution in Fourier domain but seems "incorrect" locally in time domain (see Fig. 6.9), the M.I. algorithm successfully recovers a solution that corrects the noise distribution both globally and locally (see Fig. 6.10). This feature will play a central role if one has to perform a data analysis directly in time domain after inpainting the gaps.

Prospects :

The general framework described here can be used both as a detection mean and a pre-processing step in the LISA pipeline. As it can be adapted to a wide range of gravitational events providing that the source admits a sparse representation on a specific domain, the present study paves the way for further investigations and extensions of this type of methods.

Chapter 7

Learning-based representations for GW signals

Introduction

We already emphasized in Chapter 5 (see also [1]) that non-parametric algorithms based on a simplified modeling of the sought data could already enable not only a precise detection that is robust to noise, but could also help to efficiently mitigate the impact of missing data (Chapter 6) on the analysis. One of the outcome of our algorithm was an estimate of the total sum of all GB signals detected. This could be convenient in the case where the main objective is to detect all GB signals and to subtract them from the measurements. However this approach does not allow to characterize the GBs individually ; we have yet to tackle this issue.

Designing a new non-parametric model

GB characterization requires to study GB signals individually. The model developed in Chapter 5 could approximate GB signals properly in absence of noise, but had more difficulties identifying the signal properly when the SNR became low. This can be explained by the lack of constraint the model was imposing on the sought signal. The only constraint was "sparsity in the Fourier dictionary", which is not enough to define a GB signal. Moreover, even if sparsity could compress the signal over few Fourier atoms, the number of Fourier coefficients required to have an accurate approximation was still quite high compared to the number of physical parameters needed to generate the signal. This lack of constraints made the search fast, but impaired the possibility to separate close signals and to characterize them.

Therefore, we wish to build a new model that would represent accurately the GB signals -even in presence of noise. Considering what we said before, the model should be as low-dimensional as possible, and yet must properly represent a wide range of GB signals. All this should also be parameters-agnostic, i.e. it should only rely on the signal's observed morphology.

Linear Models

The physical parameters θ of a GB are bounded for both geometrical reasons (angular parameters β, λ, ι, ψ, ϕ 0 -their bounds are defined in Table 4.1 ) and physical reasons (f 0 , ḟ0 , h 0 are constrained by the composition of the system). We can thus define Θ GB as the set of admissible parameters θ for a binary system. Now, if we consider h the signature of a GB signal observed on the TDI channels, the set of all possible signatures is defined by :

and describes a manifold which dimension is that of the space of parameters Θ GB , i.e. at most 8. In this context, using a linear model is not optimal. Indeed, we have already studied the implementation of such a model in Chapter 5 ; it was based on Fourier atoms and sparsity ensured the low-dimensionality. But it did not constrain enough the sought signals and therefore limited our ability to characterize the GBs individually.

A similar -yet different -approach would be to learn a model that would be locally linear. Typically, this would involve learning a sparse adapted dictionary (more adapted than the dictionary of Fourier atoms) ; but for the learning process to be efficient, it must rely on a large bank of templates to learn from, which is not always possible. Moreover, this would result in a large dictionary that is difficult to handle.

We could also investigate the applicability of an "optimal" linear model of the PCA type [2]. Fig. 7.1 (right) sketches the reason why such a model is unlikely to be able to represent the manifold well : the geometry of M is too complex to be approximated uniformly by a low-dimensional linear model.

We can find similar approaches in LVC. [3] proposes an approach to detect poorly modeled GW transients of short durations. The proposed linear model is based on PCA vectors that were learned from a data base. On a different note, [4; 5] used reduced order modeling of a bank of GW templates based on a linear model to accelerate detection and parameter estimation in LVC. These methods have been combined with neural networks to improve their performances [6; 7; 8]. Yet, in order to represent accurately .17 -False negative rate (FNR) as a function of SNR for IAE models (in green, labelled as "IAE") compared to the one obtained for BlockTree sparse modeling (in orange, labelled as "Blocktree"). Left :

multi-scale IAE model ; Right : Full IAE model. IAE models perform a more precise detection than the BlockTree algorithm as they are able to detect a signal with a much lower SN R.

Towards a pipeline for GB characterization

The study conducted in Section 7.3 focuses on the ability of IAE models to represent individual GB signals, centered in the analysis window -whether with or without noise. In fact, the data we expect are much more complex. Fig. 4.3 is a simulated example of a GW signal produced by all GBs in the galaxy : it is not known how many sources are present and detectable and where they are located ; moreover, the signals may have significant overlap. Under these conditions, characterizing GB signals becomes much more complex. It is necessary to develop a framework answering several questions : ■ How to localize GBs ? ■ How to separate overlapping GB signals ? ■ How to determine the number of GBs emitting in a given frequency range ? ■ How to characterize an individual signal ?

These questions are key elements to build a robust GB characterization pipeline. We propose here several answers based on IAE models.

Detection

The detection test lead in Section 7.3 focused on inputs similar to what the training set contains, i.e. centered signals to which noise is added. In a more realistic situation, even if we focus on a limited frequency range, we do not know a priori where individual GB signal stand. Previously, we showed that the estimated SN R could play the role of a selection criterion for detection. Here we want to show that it is still adapted to detect the position of a signal.

We consider the following situation : a full data set containing only one GB signal roughly located between f min and f max . IAE is applied to windows of consecutive frequency data like the one delimited in orange in Fig. 7.18 (left). This window slides from f min to f max ; for each input the IAE estimates the underlying signal and its correspon-Let us consider the training set T train with corresponding known parameters θ ∈ Θ train . We also know their barycentric coordinates in the latent space :

We would like to estimate the parameters θ test of a waveform that does not belong to the set but which parameters are still compatible with the distribution over which the training set parameters were sampled. For that purpose, we propose to make use of the latent space to recover (at least partially) parameters β and λ.

We perform a simple "nearest neighbor" estimate on the barycentric coordinates : using (BSP), we can estimate the waveform's barycentric coordinates {λ test i } i . Then we compute the distance between {λ test i } i and all of the elements of Λ train :

Ideally, the element of Λ train minimizing the distance with {λ test i } i should have parameters that are close with the sought parameters. Fig. 7.24 shows for an element of the test set the distance between {λ test i } i and the elements of Λ train as a function of parameter cos(β) and λ. The "closest neighbor" minimizing the distance has parameters that are close to the true parameters. On the overall test set, this is often verified. However, since the training set has a limited size, an estimator only based on the closest neighbor can only have limited performances. Fig. 7.25 shows the distribution of the quadratic error on the overall test set without and with noise. Without noise, there is a accurate estimation of β with quadratic error essentially inferior to 0.01. For sin(λ), the approximation is less efficient : the distribution tails are heavier, and the standard deviation is bigger. When noise is added, the estimation is further degraded with heavier distribution tails. Yet the approximation seems to be close to the real value in number of cases.

The error is mainly due to 2 elements : first and already mentioned is that the approximation is based on few elements, i.e. the elements available in the training set.

Appendices

Then we will work in the "reduced" log-posterior, restricted to intrinsic parameters :

Solving the max with regard to parameters {α j } j is simply solving the least squares minimizing the residuals between the data d and its approximation on the set of functions {h j I (int)} j .

We can define the estimator of intrinsic parameters as :

Now, if one has found the maximum of the log-posterior θ int , one can build an estimator for the extrinsic parameters based on the {α j } j . Let us introduce the amplitudes of the two GW polarizations :

Then the {α j } j are related to the extrinsic parameters as :

A.2 Solving the extrinsic parameters

Conversely, we can compute the extrinsic parameters knowing the {α j } j , by setting :

Then, we have :

and :

This can lead to consequent optimization when doing parameter estimation.

Annexe B Sparse Signal Representation B.1 Open Source code

The code is open source and can be found online at https://github.com/GW-IRFU/ gw-irfu on version 3 of the GPL (GPLv3).

B.2 Galactic binary parameters

The example of Fig. 5.7 was plot choosing the GB physical parameters contained in Tab.B.1.

Parameter Value

Frequency

# Try to group elementary blocks by 2 ;

B.4 PSD correction parameters

The polynomial correction that we used for LDC1-3 for frequencies f < 0.027 Hz is given by : Table C.1 -Waveform parameters. Signal used to show the impact of gaps on parameter estimation. This is one of the verification binaries used in LDC data sets.

Thus :

The solution of the Lagrangian problem thus linearly depends on v I , i.e. for any measurements v 1 , v 2 gapped with the same mask M , we have :

(C.2.11)

C.3 Algorithms : proofs of convergence

As we built Alg. 1 as a block coordinate descent (BCD) procedure (the variables ( h, d m )

minimize the cost function (6.2.1)), the argument of [1; 2] applies and establishes the convergence of Alg. 1. In order to prove that Alg. 2 also converges, we similarly recast it as the minimization of a cost function through a BCD procedure. Consider a fixed noise sample :

corresponding to the expected noise distribution in frequency domain, and the following algorithm :

initialized with h 0 = 0. Instead of looking for the solution as a decomposition signal/noise ( h, N), we seek the decomposition in terms of signal/missing data ( h, d m ). To this aim, we change variables similarly to Eqs. (6.2.3-6.2.4) that link the noise variable q and the missing data variable u. This relation writes :

By definition :

and the constraint on the noise solution reads :

Combining both with the change of variables formula provides :

because d g ∈ Ran(M ). Thus, the constraint on the missing data variable u writes :

Considering the alternate gapped data :

and its corresponding Fourier transform, the equation on the missing data (for classical inpainting) reads :

and thanks to the change of variables the following equality holds :

(C.3.10)

Subtracting n samp from both sides yields :

which can be plugged in the equation for h p+1 :

(C.3.12) before expanding the quadratic norm :

(C.3.13)

As we are optimizing with regard to the variable v, we can remove all the terms that are independent of it. Finally h p+1 is defined as :

The cost function :

is block-convex and its minimum can be reached through :

with initialization h 0 = 0. We recognize the form of a BCD algorithm, hence justifying the convergence of the system (C.3.2) to the global minimum of the cost function J M.I. .

C.4 Solving the noise with the Chambolle-Pock algorithm

To keep notations compact, we will respectively denote D = R N T ×2 and D = C (2N f +1)×2

the time and frequency domains.

C.4.1 Primal-dual formulation

Chambolle and Pock developed the primal-dual algorithm to solve problems with the general form :

where K is a matrix, F and G are convex functions also satisfying some extra assumptions [3] not reminded here.

Denoting by K * the conjugate operator of K and by F * , G * the conjugate applications of F, G (see [4] for a definition of conjugate applications), it was shown in [3] that using the following primal-dual algorithm : where L denotes the norm of the matrix K. The proximal function prox αf is defined for a function f and a real α > 0 by :

In the present context of gapped data in LISA, we have to minimize a function defined in the Fourier domain D (where the PSD is diagonal) subjected to the constraint u ∈ Ker(M ) expressed in the time domain D. The inpainted noise within the gaps is solution of Eq. (6.2.6) reminded here :

where n p+1 gap = d g -M h p+1 and the characteristic function 1 satisfies :

This problem is amenable to a resolution with the Chambolle-Pock algorithm with the following identification :

After some algebra involving Moreau's identity, we obtain :

(C.4.10)

C.4.2 Preconditioned formulation

The parameters τ and σ are constrained by Eq. (C.4.3) involving the norm of the matrix K, i.e. the norm of the inverse square root of the noise PSD. Since this PSD is very ill-conditioned, we will have to select very small τ, λ to satisfy Eq. (C.4.3) which leads to a slow convergence. This calls for a preconditioning of the primal-dual formulation along the lines of [5]. The main idea consists in changing the inner products equipping the time D and frequency D domains. This results in a mere modification of the proximal operators used in the iteration (C.4.2) which leaves its computational complexity basically unchanged.

We thus define new prox operators using two symmetric definite positive matrices T and Λ :

(C.4.12)

The preconditioned primal-dual algorithm writes :

with the same initialization as before at arbitrary x 0 ∈ D, y 0 ∈ D and x 0 = x 0 . This algorithm converges if the norm of the matrix Λ 1/2 KT 1/2 is (strictly) smaller than 1.

In the case of LISA data with a known PSD Σ, we chose the following :

Thanks to this choice, the prox operators can still be computed in closed form.

C.5 Effective computation of KL-Divergence

For a frequency f k , let us define a neighborhood of size n k :

Then defining :

an approximation 1 of Kullback-Leibler divergence is given by :

where V[Ω(f k )] and E[Ω(f k )] respectively denote the variance and the expectation over the elements of Ω(f k ).

For the plots, we chose a window width of ∆f = 0.1 mHz with an overlap of frequencies between two consecutive estimations of the divergence.

C.6 Inpainting behavior on a specific GB

The study was conducted choosing the parameters reported in Table C. These parameters are those needed to create a GB signal by the LDC code [6].

C.7 Algorithms Parameters

Signal extraction : We refer to [7] for details about the implementation of the signal extraction algorithm for ungapped data and the tuning of its parameters. Table D.1 -The algorithms will be tested on a set of 1, 000 waveforms which parameters have been chosen according to the the laws reported here.

Test set : 1, 000 waveforms

The network parameters are chosen as : Hence :

(D.4.7)

Now, when considering multiple noise-decorrelated channels, we define the total SN R as :

(D.4.8)

D.5 Anchor Points selection

Some of the APs selected by the process described in Section 7.

D.6 Estimating the shift

De facto, solving Eq.(multi-BSP) is not trivial -in particular because it requires to estimate the shift. To make the resolution easier, we consider that the shift can be obtained directly from comparing the observation as the delay maximizing the correlation between A (k) •Ψ {λ (k) i } i and the residual s-k-1 j=1 D τ j A (j) • Ψ {λ (j) i } i . This correlation is computed according to the method investigated in [4]. The overall minimization is done using JAX [5] minimizers that perform well in such cases.