Marco Vignuzzi

Pierre Roques

Dorothée Missé

Dominique Pontier

Barbara G Maxime

Marie Pierre

Marie Wilhelm

Marlène, Maryline Margot Lucie

De Nova-Ocampo

A Novel System for the Launch of Alphavirus RNA Synthesis Reveals a Role for the Imd Pathway in Arthropod Antiviral Response

Mes premiers remerciements vont à Dimitri, pour m'avoir donné l'opportunité tout d'abord de faire cette thèse, mais également d'avoir eu la chance de découvrir la Chine, l'Institut Pasteur de Shanghai et de pouvoir développer de magnifiques collaborations. Le parcours n'a pas été un long fleuve tranquil, mais on voit enfin le bout !! Merci à Catherine d'avoir accepté d'être ma co-directrice après le départ de Dimitri à Shanghai. Merci pour ta présence, tes nombreux conseils constructifs, nos réunions, ton soutien. Merci aux membres de mon jury, d'avoir accepté de participer à ma soutenance, à Marco Vignuzzi, Pierre Roques et Jin Zhong d'avoir accepté d'être rapporteurs, et à Dorothée Missé, et Dominique Pontier d'avoir examiné mon travail. Au cours de cette thèse, j'ai eu la chance d'être accueillie au sein de plusieurs laboratoires et équipes, ce qui a été riche en rencontres. Je remercie tout d'abord Yvan Moenne Loccoz de m'avoir permis d'intégrer le laboratoire d'Ecologie Microbienne. Merci à Patrick et Claire pour l'accueil dans l'équipe 7. A Yoann, on a commencé notre thèse en même temps, Van, pour toutes les longues commandes, Stéphanie, pour nos fous rires pendant les captures de moustiques, Guillaume, pour m'avoir montré comment trouver des larves. Une pensée toute particulière pour Flo, partie trop tôt, tu restes dans ma mémoire comme une personne aimante, soucieuse du bien-être de chacun. Merci également à tous les membres de l'unité Eco Mic, qui m'ont toujours accueilli avec le sourire. A l'unité de Gerland (Infections Virales et Pathologie Comparée), j'ai passé 5 ans mémorables, merci à tous pour votre soutien, nos nombreuses discussions et conseils. A Fabienne, pour diriger l'unité à merveille, depuis le départ (trop tôt) de Christophe. Aux membres des AlphaGirls, Carine, Marie, Céline et Catherine, merci pour votre bonne humeur, votre soutien, votre aide aussi, pour les qPCR et les titrations en grand nombre. Merci aussi à la team P3 : Maxime, Barbara G, Barbara V, Marie Pierre, merci d'avoir fait en sorte que je travaille en sécurité et dans de bonnes conditions dans ma deuxième maison : mon très cher et adoré P3, sans qui je n'aurai rien fait. Merci aux étudiants lyonnais pour nos discussions scientifiques ou pas : Marie, Wilhelm, Lucie, Marlène, Maryline et Margaux (Tea time), Nico et Nico de la Tour. Je n'oublie pas les anciens, Sarah, Franck, Claire et Najate. Pardon à ceux que

I also thank the Collaborators of the cytometry: Sébastien and Thibault, for every hour spent in the LSR, our discussions, and rich in precious advices.

My friends, who were of a big support during all these years. Charlotte, Apo, Virginie, Anna, Stephan, Val, Lucas, and the others. In spite of the years which pass and our often-fascinating lives, we are always be there for each other.

Thank you also to my family, My parents, my brother who always supported me, in spite of this "intrusive" thesis. My grandparents, I hope to make you proud. My cousin, a small thought for you, who just begins your thesis, I wish you the best for you, so many enjoyments and adventures, the good ones as the least good (it is also a part of the thesis). My family in-law also, which always had a thought for me, for your encouragements, Beatrice, Jacques, Emilie, Laurent, Mémie, Loen, Tom and the others.

The best for the end … My biggest support, Bastien, thank you for having been there during these years, to and against me. For your patience during my long weeks spent in Shanghai, or my Parisian internship, even if I know that you prefer keeping me at home. You have an incredible patience. Thank you for your support, your presence; even during the moments of doubts. This thesis is for you, who always believed in me and who pushed me to give the best. Also activated after the DENV infection, the RNA interference pathway is counteract by the action of DENV NS4, acting as an RNAi suppressor [START_REF] Kakumani | Role of RNA Interference (RNAi) in Dengue Virus Replication and Identification of NS4B as an RNAi Suppressor[END_REF]. A common pathway could be an enhancer for DEN viruses. Mutations on NS4 to affect the RNA interference regulation could reduce differences of permissiveness and viral production observed in mono and dual-infection.


Stéphane Breuils

Vincent Nozick

Akihiro Sugimoto

Eckhard Hitzer

Quadric Conformal Geometric Algebra of R 9,6

Keywords: Mathematics Subject Classification (2010). Primary 99Z99; Secondary 00A00 Quadrics, Geometric Algebra, Conformal Geometric Algebra, Clifford Algebra

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Geometric Algebra provides useful and, more importantly, intuitively understandable tools to represent, construct and manipulate geometric objects. Intensively explored by physicists, Geometric Algebra has been applied in quantum mechanics and electromagnetism [START_REF] Doran | Geometric algebra for physicists[END_REF]. Geometric Algebra has also found some interesting applications in data manipulation for Geographic Information Systems (GIS) [START_REF] Luo | A Hierarchical Representation and Computation Scheme of Arbitrary-dimensional Geometrical Primitives Based on CGA[END_REF]. More recently, it turns out that Geometric Algebra can be applied even in computer graphics, either to basic geometric primitive manipulations [START_REF] Vince | Geometric algebra for computer graphics[END_REF] or to more complex illumination processes as in [START_REF] Papaefthymiou | Real-time rendering under distant illumination with conformal geometric algebra[END_REF] where spherical harmonics are substituted by Geometric Algebra entities.

This paper presents a Geometric Algebra framework to handle quadric surfaces which can be applied to detect collision in computer graphics, and to calibrate omnidirectional cameras, usually embedding a mirror with a quadric surface, in computer vision. Handling quadric surfaces in Geometric Algebra requires a high dimensional space to work as seen in subsequent sections. Nevertheless, no low-dimensional Geometric Algebra framework is yet introduced that handles general orientation quadric surfaces and their construction using contact points.

High dimensional Geometric Algebras

Following Conformal Geometric Algebra (CGA) and its well-defined properties [START_REF] Dorst | Geometric algebra for computer science, an object-oriented approach to geometry[END_REF], the Geometric Algebra community recently started to explore new frameworks that work in higher dimensional spaces. The motivation for this direction is to increase the dimension of the relevant Euclidean space (R n with n > 3) and/or to investigate more complex geometric objects.

The Geometric Algebra of R p,q is denoted by G p,q where p is the number of basis vectors squared to +1 and q is that of basis vectors squared to -1. Then, the CGA of R 3 is denoted by G 4,1 . Extending from dimension 5 to 6 leads to G 3,3 defining either 3D projective geometry (see Dorst [START_REF] Dorst | 3d oriented projective geometry through versors of R 3,3[END_REF]) or line geometry (see Klawitter [START_REF] Klawitter | A Clifford algebraic approach to line geometry[END_REF]). Conics in R 2 are represented by the conic space of Perwass [START_REF] Perwass | Geometric algebra with applications in engineering[END_REF] with G 5,3 . Conics in R 2 are also defined by the Double Conformal Geometric Algebra (DCGA) with G 6,2 introduced by Easter and Hitzer [START_REF] Easter | Double conformal geometric algebra[END_REF]. DCGA is extended to handle cubic curves (and some other even higher order curves) in the Triple Conformal Geometric Algebra with G 9,3 [START_REF] Easter | Triple conformal geometric algebra for cubic plane curves[END_REF] and in the Double Conformal Space-Time Algebra with G 4,8 [START_REF] Benjamin | Double conformal space-time algebra[END_REF]. We note that the dimension of the algebras generated by any n-dimensional vector spaces (n = p + q) grows exponentially as they have 2 n basis elements. Although most multivectors are extremely sparse, very few implementations exist that can handle high dimensional Geometric Algebras. This problem is discussed further in Section 7.2.

Geometric Algebra and quadric surfaces

A framework to handle quadric surfaces was introduced by Zamora [START_REF] Zamora-Esquivel | G 6,3 geometric algebra; description and implementation[END_REF] for the first time. Though this framework constructs a quadric surface from control points, it supports only axis-aligned quadric surfaces.

There exist two main Geometric Algebra frameworks to manipulate general quadric surfaces.

On one hand, DCGA with G 8,2 , defined by Easter and Hitzer [START_REF] Easter | Double conformal geometric algebra[END_REF], constructs quadric and more general surfaces from their implicit equation coefficients specified by the user. A quadric (torus, Dupin-or Darboux cyclide) is represented by a bivector containing 15 coefficients that are required to construct the implicit equation of the surface. This framework preserves many properties of CGA and thus supports not only object transformations using versors but also differential operators. However, it is incapable of handling the intersection between two general quadrics and, to our best knowledge, cannot construct general quadric surfaces from control points.

On the other hand, quadric surfaces are also represented in a framework of G 4,4 as first introduced by Parkin [START_REF] Spencer | A model for quadric surfaces using geometric algebra[END_REF] and developed further by Du et al. [START_REF] Du | Modeling 3D Geometry in the Clifford Algebra R 4,4[END_REF]. Since this framework is based on a duplication of the projective geometry of R 3 , it is referred to as Double Perspective Geometric Algebra (DPGA) hereafter. DPGA represents quadric surfaces by bivector entities. The quadric expression, however, comes from a so-called "sandwiching" duplication of the product. DPGA handles quadric intersection and conics. It also handles versors transformations. However, to our best knowledge, it cannot construct general quadric surfaces from control points. This incapability seems true because, for example, wedging 9 control points together in this space results in 0 due to its vector space dimension.

Contributions

Our proposed framework, referred to as Quadric Conformal Geometric Algebra (QCGA) hereafter, is a new type of CGA, specifically dedicated to quadric surfaces. Through generalizing the conic construction in R 2 by Perwass [START_REF] Perwass | Geometric algebra with applications in engineering[END_REF], QCGA is capable of constructing quadric surfaces using either control points or implicit equations. Moreover, QCGA can compute the intersection of quadric surfaces, the surface tangent, and normal vectors for a quadric surface point.

Notation

We use the following notation throughout the paper. Lower-case bold letters denote basis blades and multivectors (multivector a). Italic lower-case letters refer to multivector components (a 1 , x, y 2 , • • • ). For example, a i is the i th coordinate of the multivector a. Constant scalars are denoted using lowercase default text font (constant radius r). The superscript star used in x * represents the dualization of the multivector x. Finally, subscript on x refers to the Euclidean vector associated with the point x of QCGA.

Note that in geometric algebra, the inner product, contractions and outer product have priority over the full geometric product. For instance, a ∧ bI = (a ∧ b)I.

QCGA definition

This section introduces QCGA. We specify its basis vectors and give the definition of a point.

QCGA basis and metric

The QCGA G 9,6 is defined over a 15-dimensional vector space. The base vectors of the space R 9,6 are basically divided into three groups: {e 1 , e 2 , e 3 } (corresponding to the Euclidean vectors in R 3 ), {e o1 , e o2 , e o3 , e o4 , e o5 , e o6 }, and {e ∞1 , e ∞2 , e ∞3 , e ∞4 , e ∞5 , e ∞6 }. The inner products between them are as defined in Table 1.

For some computation constraints, a diagonal metric matrix may be required. The orthonormal vector basis of R

• • • • • • • • • • • • e 2 0 1 0 • • • • • • • • • • • • e 3 0 0 1 • • • • • • • • • • • • e o1 • • • 0 -1 • • • • • • • • • • e ∞1 • • • -1 0 • • • • • • • • • • e o2 • • • • • 0 -1 • • • • • • • • e ∞2 • • • • • -1 0 • • • • • • • • e o3 • • • • • • • 0 -1 • • • • • • e ∞3 • • • • • • • -1 0 • • • • • • e o4 • • • • • • • • • 0 -1 • • • • e ∞4 • • • • • • • • • -1 0 • • • • e o5 • • • • • • • • • • • 0 -1 • • e ∞5 • • • • • • • • • • • -1 0 • • e o6 • • • • • • • • • • • • • 0 -1 e ∞6 • • • • • • • • • • • • • -1 0
squares to +1 along with six other basis vectors {e -1 , e -2 , e -1 , e -4 , e -5 , e -6 } each of which squares to -1 corresponds to a diagonal metric matrix. The transformation from the original basis to this new basis (with diagonal metric) can be defined as follows:

e ∞i = e +i + e -i , e oi = 1 2 (e -i -e +i ), i ∈ {1, • • • , 6}. (2.1)
For clarity, we also define the 6-blades

I ∞ = e ∞1 ∧ e ∞2 ∧ e ∞3 ∧ e ∞4 ∧ e ∞5 ∧ e ∞6 , I o = e o1 ∧ e o2 ∧ e o3 ∧ e o4 ∧ e o5 ∧ e o6 , (2.2)
the 5-blades

I ∞ = (e ∞1 -e ∞2 ) ∧ (e ∞2 -e ∞3 ) ∧ e ∞4 ∧ e ∞5 ∧ e ∞6 , I o = (e o1 -e o2 ) ∧ (e o2 -e o3 ) ∧ e o4 ∧ e o5 ∧ e o6 , (2.3)
the pseudo-scalar of

R 3 I = e 1 ∧ e 2 ∧ e 3 , (2.4)
and the pseudo-scalar

I = I ∧ I ∞ ∧ I o . (2.5)
The inverse of the pseudo-scalar results in

I -1 = -I. (2.6)
The dual of a multivector indicates division by the pseudo-scalar, e.g., a * = -aI, a = a * I. From eq. (1.19) in [START_REF] Hitzer | Carrier method for the general evaluation and control of pose, molecular conformation, tracking, and the like[END_REF], we have the useful duality between outer and inner products of non-scalar blades a and b in Geometric Algebra:

(a ∧ b) * = a • b * , a ∧ (b * ) = (a • b) * , a ∧ (bI) = (a • b)I, (2.7)
which indicates that

a ∧ b = 0 ⇔ a • b * = 0, a • b = 0 ⇔ a ∧ b * = 0. (2.8)

Point in QCGA

The point x of QCGA corresponding to the Euclidean point x = xe 1 +ye 2 + ze 3 ∈ R 3 is defined as

x = x + 1 2 (
x 2 e ∞1 +y 2 e ∞2 +z 2 e ∞3 )+xye ∞4 +xze ∞5 +yze ∞6 +e o1 +e o2 +e o3 .

(2.9) Note that the null vectors e o4 , e o5 , e o6 are not present in the definition of the point. This is merely to keep the convenient properties of CGA points, namely, the inner product between two points is identical with the squared distance between them. Let x 1 and x 2 be two points, their inner product is from which together with Table 1, it follows that

x 1 • x 2 = x 1 e 1 + y 1 e 2 + z 1 e 3 + 1 2 x 2 1 e ∞1 + 1 2 y 2 1 e ∞2 + 1 2 z 2 1 e ∞3 + x 1 y 1 e ∞4 + x 1 z 1 e ∞5 +
x 1 • x 2 = x 1 x 2 + y 1 y 2 + z 1 z 2 - 1 2 x 2 1 - 1 2 x 2 2 - 1 2 y 2 1 - 1 2 y 2 2 - 1 2 z 2 1 - 1 2 z 2 2 = - 1 2 x 1 -x 2 2 .
(2.11)

We see that the inner product is equivalent to minus half the squared Euclidean distance between x 1 and x 2 .

QCGA objects

QCGA is an extension of CGA, thus the objects defined in CGA are also defined in QCGA. The following sections explore the plane, the line, and the sphere to show their definitions in QCGA, and similarity between these objects in CGA and their counterparts in QCGA.

3.1. Plane 3.1.1. Primal plane. As in CGA, a plane π in QCGA is computed using the wedge of three linearly independent points x 1 , x 2 , and x 3 on the plane:

π = x 1 ∧ x 2 ∧ x 3 ∧ I ∞ ∧ I o . (3.1)
The multivector π corresponds to the primal form of a plane in QCGA, with grade 14, composed of six components. The e o2o3 , e o1o3 , e o1o2 components have the same coefficient and can thus be factorized, resulting in a form defined with only four coefficients x n , y n , z n and h: Using distributivity and anticommutativity of the outer product, we obtain

π = x
x ∧ π = xx n + yy n + zz n - 1 3 h(1 + 1 + 1) I = xx n + yy n + zz n -h I = (x • n -h) I, (3.4)
which corresponds to the Hessian form of the plane with Euclidean normal n = x n e 1 + y n e 2 + z n e 3 and with orthogonal distance h from the origin. Proof. Consequence of (2.8).

Because of (2.11), a plane can also be obtained as the bisection plane of two points x 1 and x 2 in a similar way as in CGA.

Proposition 3.3. The dual plane π * = x 1 -x 2 is the dual orthogonal bisecting plane between the points x 1 and x 2 .

Proof. From Proposition 3.2, every point x on π * satisfies x • π * = 0,

x • (x 1 -x 2 ) = x • x 1 -x • x 2 = 0. (3.6)
As seen in (2.11), the inner product between two points results in the squared Euclidean distance between the two points. We thus have

x • (x 1 -x 2 ) = 0 ⇔ x -x 1 2 = x -x 2 2 . (3.7)
This corresponds to the equation of the orthogonal bisecting dual plane between x 1 and x 2 .

3.2. Line 3.2.1. Primal line. A primal line l is a 13-vector constructed from two linearly independent points x 1 and x 2 as follows:

l = x 1 ∧ x 2 ∧ I ∞ ∧ I o . (3.8)
The outer product between the 6-vector I ∞ and the two points x 1 and x 2 removes all their e ∞i components (i ∈ {1, • • • , 6}). Accordingly, they can be reduced to x 1 = (e o1 + e o2 + e o3 + x 1 ) and x 2 = (e o1 + e o2 + e o3 + x 2 ) respectively. For clarity, (3.8) is simplified "in advance" as

l = x 1 ∧ (e o1 + e o2 + e o3 + x 2 ) ∧ I ∞ ∧ (e o1 -e o2 ) ∧ (e o2 -e o3 ) ∧ e o4o5o6 = x 1 ∧ (x 2 ∧ (e o1 -e o2 ) ∧ (e o2 -e o3 ) + 3e o1 ∧ e o2 ∧ e o3 ) ∧ I ∞ ∧ e o4o5o6 = 3e o1 ∧ e o2 ∧ e o3 ∧ (x 2 -x 1 ) + x 1 ∧ x 2 ∧ (e o1 -e o2 ) ∧ (e o2 -e o3 ) ∧ I ∞ ∧ e o4o5o6 .
(3.9)

Setting u = x 2 -x 1 and v = x 1 ∧ x 2 gives l = 3e o1 ∧ e o2 ∧ e o3 ∧ u + v ∧ (e o1o2 -e o1o3 + e o2o3 ) ∧ I ∞ ∧ e o4o5o6 = -3 u I ∞ ∧ I o + v I ∞ ∧ I o . (3.10)
Note that u and v correspond to the 6 Plücker coefficients of a line in R 3 . More precisely, u is the support vector of the line and v is its moment.

Proposition 3.4. A point x with Euclidean coordinates x lies on the line l iff x ∧ l = 0.

Proof.

x ∧ l = (x + e o1 + e o2 + e o3 ) ∧ (-

3 u I ∞ ∧ I o + v I ∞ ∧ I o ) = -3x ∧ u I ∞ ∧ I o + x ∧ v I ∞ ∧ I o + v I ∞ ∧ (e o1 + e o2 + e o3 ) ∧ I o = -3(x ∧ u -v ) I ∞ ∧ I o + x ∧ v I ∞ ∧ I o . (3.11)
The 6-blade I ∞ ∧ I o and the 5-blade I ∞ ∧ I o are linearly independent. Therefore, x ∧ l = 0 yields

x ∧ l = 0 ⇔ x ∧ u = v , x ∧ v = 0. (3.12)
As x , u and v are Euclidean entities, (3.12) corresponds to the Plücker equations of a line [START_REF] Kanatani | Understanding geometric algebra: Hamilton, Grassmann, and Clifford for Computer Vision and Graphics[END_REF].

Dual line.

Dualizing the entity l consists in computing with duals: Proof. Consequence of (2.8).

l * = (-3 u I ∞ ∧ I o + v I ∞ ∧ I o )(-I) = 3 u I + (e ∞3 + e ∞2 + e ∞1 ) ∧ v I . ( 3
Note that a dual line l * can also be constructed from the intersection of two dual planes as follows:

l * = π * 1 ∧ π * 2 .
(3.14)

3.3. Sphere 3.3.1. Primal sphere. We define a sphere s using four points as the 14-blade

s = x 1 ∧ x 2 ∧ x 3 ∧ x 4 ∧ I ∞ ∧ I o . (3.15)
The outer product of the points with I ∞ removes all e ∞4 , e ∞5 , e ∞6 components of these points, i.e., the cross terms (xy, xz, and yz). The same remark holds for I o and e o4 , e o5 , e o6 . For clarity, we omit these terms below. We thus have

s = x 1 ∧x 2 ∧x 3 ∧ 1
s =x 1 ∧ x 2 ∧ (x 3 ∧ x 4 I o ∧ I ∞ + 3(x 3 -x 4 )I o ∧ I ∞ + 1 2 x 4 2 x 3 I ∞ ∧ I o - 1 2 x 3 2 x 4 I ∞ ∧ I o (3.17) + 3 2 ( x 4 2 -x 3 2 )I o ∧ I ∞ .
Again we remark that the resulting entity has striking similarities with a point pair of CGA. More precisely, let c be the Euclidean midpoint between the two entities x 3 and x 4 , d be the unit vector from x 3 to x 4 , and r be half of the Euclidean distance between the two points in exactly the same way as Hitzer et al in [START_REF] Hitzer | Carrier method for the general evaluation and control of pose, molecular conformation, tracking, and the like[END_REF], namely

2r = |x 3 -x 4 | , d = x 3 -x 4 2r , c = x 3 + x 4 2 . (3.18)
Then, (3.17) can be rewritten by

s =x 1 ∧ x 2 ∧ 2r d ∧ c I o ∧ I ∞ (3.19) + 3d I o ∧ I ∞ + 1 2 (c 2 + r 2 )d -2c c • d I ∞ ∧ I o .
The bottom part corresponds to a point pair, as defined in [START_REF] Hitzer | Carrier method for the general evaluation and control of pose, molecular conformation, tracking, and the like[END_REF], that belongs to the round object family. Applying the same development to the two points x 1 and x 2 again results in round objects:

s = - 1 6 x c 2 -r 2 I ∧ I ∞ ∧ I o + e 123 ∧ I ∞ ∧ I o + (x c I ) ∧ I ∞ ∧ I o . (3.20)
Note that x c corresponds to the center point of the sphere and r to its radius. It can be further simplified into

s = x c - 1 6 r 2 (e ∞1 + e ∞2 + e ∞3 ) I, (3.21)
which is dualized to

s * = x c - 1 6 r 2 (e ∞1 + e ∞2 + e ∞3 ), (3.22)
where x c corresponds to x c without the cross terms xy, xz, yz. Since a QCGA point has no e o4 , e o5 , e o6 components, building a sphere with these cross terms is also valid. However, inserting these cross terms (that actually do not appear in the primal form) raises some issues in computing intersections with other objects.

Proposition 3.6. A point x lies on the sphere s iff x ∧ s = 0.

Proof. Since the components e ∞4 , e ∞5 and e ∞6 of x are removed by the outer product with s of (3.17), we ignore them to obtain

x ∧ s = x ∧ (s * I) = x • s * I (3.23) = x + e o1 + e o2 + e o3 + 1 2 x 2 e ∞1 + 1 2 y 2 e ∞2 + 1 2 z 2 e ∞3 (3.24) • x c - 1 6 r 2 (e ∞1 + e ∞2 + e ∞3 ) I,
which can be rewritten by

x ∧ s = xx c + yy c + zz c - 1 2 x 2 c - 1 6 r 2 - 1 2 y 2 c - 1 6 r 2 - 1 2 z 2 c - 1 6 r 2 - 1 2 x 2 - 1 2 y 2 - 1 2 z 2 I = 0. (3.25)
This can take a more compact form defining a sphere Proof. Consequence of (2.8).

(x -x c ) 2 + (y -y c ) 2 + (z -z c ) 2 = r 2 . ( 3

Quadric surfaces

This section describes how QCGA handles quadric surfaces. All QCGA objects defined in Section 3 become thus part of a more general framework.

Primal quadric surfaces

The implicit formula of a quadric surface in R 3 is F (x, y, z) = ax 2 + by 2 + cz 2 + dxy + exz + fyz + gx + hy + iz + j = 0. (4.1)

A quadric surface is constructed by wedging 9 points together with 5 null basis vectors as follows

q = x 1 ∧ x 2 ∧ • • • ∧ x 9 ∧ I o . (4.2)
The multivector q corresponds to the primal form of a quadric surface with grade 14 and 12 components. Again 3 of these components have the same coefficient and can be combined together into the form defined by 10 coefficients a, b, . . . , j, as in q = e 123 2ae o1 + 2be o2 + 2ce o3 + de o4 + ee o5 + fe o6 where in the second equality we used the duality property. The expression for the dual quadric vector is therefore q * = -2ae o1 + 2be o2 + 2ce o3 + de o4 + ee o5 + fe o6

+ ge 1 + he 2 + ie 3 -j 3 (e ∞1 + e ∞2 + e ∞3 ). (4.4)

Proposition 4.1. A point x lies on the quadric surface q iff x ∧ q = 0.

Proof.

x ∧ q = x ∧ (q * I) = (x • q * )I = x • -2ae o1 + 2be o2 + 2ce o3 + de o4 + ee o5 + fe o6

+ ge 1 + he 2 + ie 3 -j 3 (e ∞1 + e ∞2 + e ∞3 ) I = ax 2 + by 2 + cz 2 + dxy + exz + fyz + gx + hy + iz + j I. (

This corresponds to the formula (4.1) representing a general quadric surface.

Dual quadric surfaces

The dualization of a primal quadric surface leads to the 1-vector dual quadric surface q * of (4.4). We have the following proposition whose proof is a consequence of (2.8).

Proposition 4.2. A point x lies on the dual quadric surface q * iff x • q * = 0.

Normals and tangents

This section presents the computation of the normal Euclidean vector n and the tangent plane π * of a point x (associated to the Euclidean point x = xe 1 + ye 2 + ye 3 ) on a dual quadric surface q * . The implicit formula of the dual quadric surface is considered as the following scalar field

F (x, y, z) = x • q * . ( 5.1)
The normal vector n of a point x is computed as the gradient of the implicit surface (scalar field) at x:

n = ∇F (x, y, z) = ∂F (x, y, z) ∂x e 1 + ∂F (x, y, z) ∂y e 2 + ∂F (x, y, z) ∂z e 3 . (5.2)
Since the partial derivative with respect to the x component is defined by

∂F (x, y, z) ∂x = lim h →0 F (x + h, y, z) -F (x, y, z) h , (5.3)
we have

∂F (x, y, z) ∂x = lim h →0 x 2 • q * -x • q * h = lim h →0 x 2 -x h • q * , (5.4)
where x 2 is the point obtained by translating x along the x-axis by the value h. Note that x 2 -x represents the dual orthogonal bisecting plane spanned by x 2 and x (see Proposition 3.3). Accordingly, we have

lim h →0 x 2 -x h = xe ∞1 + ye ∞4 + ze ∞5 + e 1 = (x • e 1 )e ∞1 + (x • e 2 )e ∞4 + (x • e 3 )e ∞5 + e 1 . (5.5)
This argument can also be applied to the partial derivative with respect to the y and z components. Therefore, we obtain

n = (x • e 1 )e ∞1 + (x • e 2 )e ∞4 + (x • e 3 )e ∞5 + e 1 • q * e 1 + (x • e 2 )e ∞2 + (x • e 1 )e ∞4 + (x • e 3 )e ∞6 + e 2 • q * e 2 + (x • e 3 )e ∞3 + (x • e 1 )e ∞5 + (x • e 2 )e ∞6 + e 3 • q * e 3 . (5.6)
On the other hand, the tangent plane at a surface point x can be computed from the Euclidean normal vector n and the point x. Since the plane orthogonal distance from the origin is -2(e o1 + e o3 + e o3 ) • x, the tangent plane π * is obtained as

π * = n + 1 3
e ∞1 + e ∞2 + e ∞3 -2(e o1 + e o3 + e o3 ) • x.

(5.7)

Intersections

Let us consider two geometric objects corresponding to dual quadrics1 a * and b * . Assuming that the two objects are linearly independent, i.e., a * and b * are linearly independent, we consider the outer product c * of these two objects c * = a * ∧ b * . (6.1) If a point x lies on c * , then

x • c * = x • (a * ∧ b * ) = 0. (6.2)
The inner product computation of (6.2) leads to

x • c * = (x • a * )b * -(x • b * )a * = 0. (6.3)
Our assumption of linear independence between a * and b * indicates that (6.3) holds if and only if x • a * = 0 and x • b * = 0, i.e. the point x lies on both quadrics. Thus, c * = a * ∧ b * represents the intersection of the linearly independent quadrics a * and b * , and a point x lies on this intersection if and only if x • c * = 0.

6.1. Quadric-Line intersection For example, in computer graphics, making a Geometric Algebra compatible with a raytracer requires only to be able to compute a surface normal and a line-object intersection. This section defines the line-quadric intersection.

Similarly to (6.1), the intersection x ± between a dual line l * and a dual quadric q * is computed by l * ∧ q * . Any point x lying on the line l defined by two points x 1 and x 2 can be represented by the parametric formula x = α(x 1 -x 2 ) + x 2 = αu + x 2 . Note that u could also be computed directly from the dual line l * (see (3.13)). Any point x 2 ∈ l can be used, especially the closest point of l from the origin, i.e. x 2 = v •u -1 . Accordingly, computing the intersection between the dual line l * and the dual quadric q * becomes equivalent to finding α such that x lies on the dual quadric, i.e., x • q * = 0, leading to a second degree equation in α, as shown in (4.1). In this situation, the problem is reduced to computing the roots of this equation. However, we have to consider four cases: the case where the line is tangent to the quadric, the case where the intersection is empty, the case where the line intersects the quadric into two points, and the case where one of the two points exists at infinity. To identify each case, we use the discriminant δ defined as:

δ = β 2 -4(x 2 • q * ) 6 i=1 (u • e oi )(q * • e ∞i ), (6.4)
where

β = 2u • (a(x 2 • e 1 )e 1 + b(x 2 • e 2 )e 2 + c(x 2 • e 3 )e 3 )+ d (u ∧ e 1 ) • (x 2 ∧ e 2 ) + (x 2 ∧ e 1 ) • (u ∧ e 2 ) + e (u ∧ e 1 ) • (x 2 ∧ e 3 ) + (x 2 ∧ e 1 ) • (u ∧ e 3 ) + f (u ∧ e 2 ) • (x 2 ∧ e 3 ) + (x 2 ∧ e 2 ) • (u ∧ e 3 ) + q * • u . (6.5)
If δ < 0, the line does not intersect the quadric (the solutions are complex). If δ = 0, the line and the quadric are tangent. If δ > 0 and 6 i=1 (u • e oi )(q * • e ∞i ) = 0, we have only one intersection point (linear equation). Otherwise, we have two different intersection points x ± computed by

x ± = u(-β ± √ δ)/ 2 6 i=1 (u • e oi )(q * • e ∞i ) + x 2 . ( 6
• • quadrics intersection • • • quadric plane intersection • • • versors • • Darboux cyclides • • •
in DPGA and DCGA. However, QCGA also faces some limitations that do not affect DPGA and DCGA, as summarized in Table 2. First, DPGA and DCGA are known to be capable of transforming all objects by versors [START_REF] Du | Modeling 3D Geometry in the Clifford Algebra R 4,4[END_REF][START_REF] Easter | Double conformal geometric algebra[END_REF] whereas it is not yet clear whether objects in QCGA can be transformed using versors. An extended version of CGA versors can be used to transform lines in QCGA (and probably all round and flat objects of CGA), but more investigation is needed. Second, the number of basis elements spanned by QCGA is 2 15 ( 32, 000) components for a full multivector. Although multivectors of QCGA are in reality almost always very sparse, this large number of elements may cause implementation issues (see Section 7.2). It also requires some numerical care in computation, especially during the wedge of 9 points. This is because some components are likely to be multiplied at the power of 9.

Implementations

There exist many different implementations of Geometric Algebra, however, very few can handle dimensions higher than 8 or 10. This is because higher dimensions bring a large number of elements of the multivector used, resulting in expensive computation. In many cases, the computation then becomes impossible in practice. QCGA has a 15 vector space dimension and hence requires some specific care during the computation.

We conducted our tests with an enhanced version of Breuils et al. [START_REF] Breuils | A geometric algebra implementation using binary tree[END_REF][START_REF] Breuils | A hybrid approach for computing products of high-dimensional geometric algebras[END_REF] which is based on a recursive framework. We remark that most of the products involved in our tests were the outer products between 14-vectors and 1-vectors, applying one of the less time consuming products of QCGA. Indeed, QCGA with vector space dimension of 15 has 2 15 elements and this number is 1,000 times as large as that of elements for CGA with vector space dimension of 5 (CGA with vector space dimension of 5 is needed for the equivalent operations with QCGA with dimension of 15). The computational time required for QCGA, however, did not need 1,000 times but only 70 times of that for CGA. This means that the computation of QCGA runs in reasonable time on the enhanced version of Breuils et al. [START_REF] Breuils | A geometric algebra implementation using binary tree[END_REF][START_REF] Breuils | A hybrid approach for computing products of high-dimensional geometric algebras[END_REF]. More detailed analysis in this direction is left for future work. Figure 1 depicts a few examples generated with our OpenGL renderer based on the outer product null-space voxels and our ray-tracer. From left to right: a dual hyperboloid built from its equation, an ellipsoid built from its control points (in yellow), the intersection between two cylinders, and a hyperboloid with an ellipsoid and planes (the last one was computed with our ray-tracer).

Conclusion

This paper presented a new Geometric Algebra framework, Quadric Conformal Geometric Algebra (QCGA), that handles the construction of quadric surfaces using the implicit equation and/or control points. QCGA naturally includes CGA objects and generalizes some dedicated constructions. The intersection between objects in QCGA is computed using only outer products. This paper also detailed the computation of the tangent plane and the normal vector at a point on a quadric surface. Although QCGA is defined in a high dimensional space, most of the computations run in relatively low dimensional subspaces of this framework. Therefore, QCGA can be used for numerical computations in applications such as computer graphics.

3. 1 . 2 .

 12 Dual plane. The dualization of the primal form of the plane is π * = n + 1 3 h(e ∞1 + e ∞2 + e ∞3 ). (3.5) Proposition 3.2. A point x with Euclidean coordinates x lies on the dual plane π * iff x • π * = 0.

. 13 ) 3 . 5 .

 1335 Proposition A point x lies on the dual line l * iff x • l * = 0.

.26) 3 . 3 . 2 .

 332 Dual sphere. The dualization of the primal sphere s gives: s * = x c -1 6 r 2 (e ∞1 + e ∞2 + e ∞3 ). (3.27) Proposition 3.7. A point x lies on the dual sphere s * iff x • s * = 0.

Figure 1 .

 1 Figure 1. Example of our construction of QCGA objects.From left to right: a dual hyperboloid built from its equation, an ellipsoid built from its control points (in yellow), the intersection between two cylinders, and a hyperboloid with an ellipsoid and planes (the last one was computed with our ray-tracer).

Table 1 .

 1 9,6 with the Euclidean basis {e 1 , e 2 , e 3 }, and 6 basis vectors {e +1 , e +2 , e +3 , e +4 , e +5 , e +6 } each of which Inner product between QCGA basis vectors.e 1 e 2 e 3 e o1 e ∞1 e o2 e ∞2 e o3 e ∞3 e o4 e ∞4 e o5 e ∞5 e o6 e ∞6

	e 1	1 0 0

  y 1 z 1 e ∞6 + e o1 + e o2 + e o3 • x 2 e 1 + y 2 e 2 + z 2 e 3 +

	1 2	x 2 2 e ∞1 +	1 2	y 2 2 e ∞2 +	1 2	z 2 2 e ∞3

+ x 2 y 2 e ∞4 + x 2 z 2 e ∞5 + y 2 z 2 e ∞6 + e o1 + e o2 + e o3 .

(2.10)

  Proposition 3.1. A point x with Euclidean coordinates x lies on the plane π iff x ∧ π = 0. ∞3 + xye ∞4 + xze ∞5 + yze ∞6 + e o1 + e o2 + e o3 ∧ (x n e 23 -y n e 13 + z n e 12 )I ∞ ∧ I o

	Proof.
	x ∧ π = xe 1 + ye 2 + ze 3 + z 2 e + 1 2 x 2 e ∞1 + 1 2 y 2 e ∞2 + 1 2 h 3 e 123 I ∞ ∧ (e o2o3 -e o1o3 + e o1o2 ) ∧ e o4o5o6 .	(3.3)

n e 23 -y n e 13 + z n e 12 I ∞ ∧ I o + h 3 e 123 I ∞ ∧ e o2o3 -e o1o3 + e o1o2 ∧ e o4o5o6 . (3.2)

  ∞3 )I ∞ ∧I o -3I ∞ ∧I o +x 4 I ∞ ∧I o . (3.16)Note the similarities with a CGA point x 4 + e o + 1 2 x 4 2 e ∞ . Then, the explicit outer product with x 3 gives:

	2	(x 2 4 e ∞1 +y 2 4 e ∞2 +z 2 4 e

  • I ∞ ∧ I o + ge 1 + he 2 + ie 3 e 123 I ∞ ∧ I o + j 3 e 123 I ∞ ∧ (e ∞1 + e ∞2 + e ∞3 ) • I o = -2ae o1 + 2be o2 + 2ce o3 + de o4 + ee o5 + fe o6 + ge 1 + he 2 + ie 3 -j 3 e ∞1 + e ∞2 + e ∞3 I = q * I, (4.3)

Table 2 .

 2 Comparison of properties between DPGA, DCGA, and QCGA. The symboli • stands for "capable", • for "incapable" and for "unknown".

	.6)

7.1. Limitations

The construction of quadric surfaces by the wedge of conformal points presented in Sections 3 and 4 is a distinguished property of QCGA that is missing

The term "quadric" (without being followed by surface) encompasses quadric surfaces and conic sections.


C Kervazo

J Bobin

C Chenot

Blind separation of a large number of sparse sources

Keywords: Blind source separation, sparse representations, block-coordinate optimization strategies, matrix factorization

Blind separation of a large number of sparse sources

Introduction

Problem statement

Blind source separation (BSS) is the major analysis tool to retrieve meaningful information from multichannel data. It has been particularly successful in a very wide range of signal processing applications ranging from astrophysics [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF] to spectroscopic data in medicine [START_REF] Biswal | Blind source separation of multiple signal sources of fMRI data sets using independent component analysis[END_REF] or nuclear physics [START_REF] Nuzillard | Application of blind source separation to 1-D and 2-D nuclear magnetic resonance spectroscopy[END_REF], to name only a few. In this framework, the observations {x i } i=1,...,m are modeled as a linear combination of n unknown elementary sources {s j } j=1,...,n :

x i = n j=1 a ij s j + z i . The coefficients a ij are measuring the contribution of the j-th source to the observation x i , while z i is modeling an additive noise as well as model imperfections. Each datum x i and source s j is supposed to have t entries. This problem can be readily recast in a matrix formulation:

X = AS + N (1)
where X is a matrix composed of the m row observations and t columns corresponding to the entries (or samples), the mixing matrix A is built from the {a ij } i=1,...,m,j=1,...,n coefficients and S is a n × t matrix containing the sources. Using this formulation, the goal of BSS is to estimate the unknown matrices A and S from the sole knowledge of X.

Blind source separation methods

It is well-known that BSS is an ill-posed inverse problem, which requires additional prior information on either A or S to be tackled [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. Making BSS a better-posed problem is performed by promoting some discriminant information or diversity among the sources. A first family of standard techniques, such as Independent Component Analysis (ICA), assumes that the sources are statistically independent [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF].

In this study, we will specifically focus on the family of algorithms dealing with the case of sparse BSS problems (i.e. where the sources are assumed to be sparse), which have attracted a lot of interest during the last two decades [START_REF] Zibulevsky | Blind source separation by sparse decomposition in a signal dictionary[END_REF][START_REF] Bronstein | Sparse ICA for blind separation of transmitted and reflected images[END_REF][START_REF] Li | Underdetermined blind source separation based on sparse representation[END_REF]. Sparse BSS has mainly been motivated by the success of sparse signal modeling for solving very large classes of inverse problems [START_REF] Starck | Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity[END_REF]. The Generalized Morphological Component Analysis (GMCA) algorithm [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF][START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF] builds upon the concept of morphological diversity to disentangle sources that are assumed to be sparsely distributed in a given dictionary. The morphological diversity property states that sources with different morphologies are unlikely to have similar large value coefficients. This is the case of sparse and independently distributed sources, with high probability. In the framework of Independent Component Analysis (ICA), Efficient FastICA (EFICA) [START_REF] Koldovsky | Efficient variant of algorithm Fas-tICA for independent component analysis attaining the Cramér-Rao lower bound[END_REF] is a FastICA-based algorithm that is especially adapted to retrieve sources with generalized Gaussian distributions, which includes sparse sources. In the seminal paper [START_REF] Zibulevsky | Blind source separation with relative Newton method[END_REF], the author also proposed a Newton-like method for ICA called Relative Newton Algorithm (RNA), which uses quasi-maximum likelihood estimation to estimate sparse sources. A final family of algorithms builds on the special case where it is known that A and S are furthermore non-negative, which is often the case on real world data [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF].

However, the performances of most of these methods decline when the number of sources n becomes large. As an illustration, Fig. 1 shows the evolution of the mixing matrix criterion (cf. sec. 3.1, [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF]) as a function of the number of sources for various BSS methods. This experiment illustrates that most methods do not perform correctly in the "large-scale"regime. In this case, the main source of deterioration is very likely related to the non-convex nature of BSS. Indeed, for a fixed number of samples t, an increasing number of sources n will make these algorithms more prone to be trapped in spurious local minima, which tends to hinder the applicability of BSS on practical issues with a large n. Consequently, the optimization strategy has a huge impact on the separation performances.

Contribution

In a large number of applications such as astronomical [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF] or biomedical signals [START_REF] Biswal | Blind source separation of multiple signal sources of fMRI data sets using independent component analysis[END_REF], designing BSS methods that are tailored to precisely retrieve a large number of sources is of paramount importance. For that purpose, the goal of this article is to introduce a novel algorithm dubbed bGMCA (block-Generalized Morphological Component Analysis) to specifically tackle sparse BSS problems when a large number of sources need to be estimated.

In this setting, which we will later call the large-scale regime, the algorithmic strategy has a huge impact on the separation quality since BSS requires solving highy challenging non-convex problems. For that purpose, the proposed bGMCA algorithm builds upon the sparse modeling of the sources, as well as an efficient minimization scheme based on block-coordinate descent. In contrast to state-of-the art methods [START_REF] Zibulevsky | Blind source separation with relative Newton method[END_REF][START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF][START_REF] Rapin | NMF with sparse regularizations in transformed domains[END_REF][START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF], we show that making profit of block-based minimization with intermediate block sizes allows the bGMCA to dramatically enhance the separation performances, particularly when the number of sources to be estimated becomes large. Comparisons with the state-of-the art methods have been carried out on various simulation scenarios. The last part of the article will show the flexibility of bGMCA, with an application to sparse and non-negative BSS in the context of spectroscopy.

Optimization problem and bGMCA

General problem

Sparse BSS [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF][START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF] aims to estimate the mixing matrix A and the sources S by minimizing a penalized least-squares of the form:

min A,S 1 2 X -AS 2 F + J (A) + G(S) (2)
The first term is a classical data fidelity term that measures the discrepancy between the data and the mixture model. The . F norm refers to the Frobenius norms, whose use stems from the assumption that the noise is Gaussian. The penalizations J and G enforce some desired properties on A and S (e.g. sparsity, non-negativity). In the following, we will consider that the proximal operators of J and G are defined, and that J and G are convex. However, the whole matrix factorization problem (2) is non-convex.

Consequently, the strategy of optimization has a critical impact on the separation performances, especially to avoid spurious local minimizers and to reduce the sensitivity to initialization. A common idea of several strategies (Block Coordinate Relaxation -BCR [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF], Proximal Alternating Linearized Minimization -PALM [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], Alternating Least Squares -ALS) is to benefit from the multi-convex structure of (2) by using blocks [START_REF] Xu | A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF] in which each sub-problem is convex. The minimization is then performed alternately with respect to one of the coordinate blocks while the other coordinates stay fixed, which entails solving a sequence of convex optimization problems. Most of the already existing methods can then be categorized in one of two families, depending on the block sizes:

-Hierarchical or deflation methods: these algorithms use a block of size 1. For instance, Hierarchical ALS (HALS) ( [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF] and references therein)

updates only one specific column of A and one specific row of S at each iteration. The main advantage of this family is that each subproblem is often much simpler as their minimizer generally admits a closed-form expression. Moreover, the matrices involved being small, the computation time is much lower. The drawback is however that the errors on some sources/mixing matrix columns propagate from one iteration to the other since they are updated independently.

-Full-size blocks: these algorithms use as blocks the whole matrices A and S (the block size is thus equal to n). For instance, GMCA [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF],

which is reminiscent of the projected Alternating Least Squares (pALS) algorithm, is part of this family. One problem compared to hierarchical or deflation methods is that the problem is more complex due to the simultaneous estimation of a high number of sources. Moreover, the computational cost increases quickly with the number of sources.

The gist of the proposed bGMCA algorithm is to adopt an alternative approach that uses intermediate block sizes. The underlying intuition is that using blocks of intermediate size can be recast as small-scale source separation problems, which are simpler to solve as testified by Fig. 1. As a byproduct, small-size subproblems are also less costly to tackle.

Block based optimization

In the following, bGMCA minimizes the problem in eq. ( 2) with blocks, which are indexed by a set of indices I of size r, 1 r n. In practice, the minimization is performed at each iteration on submatrices of A (keeping only the columns indexed by I ) and S (keeping only the rows indexed by I ).

Minimizing multi-convex problems

Block coordinate relaxation (BCR, [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]) is performed by minimizing [START_REF] Biswal | Blind source separation of multiple signal sources of fMRI data sets using independent component analysis[END_REF] according to a single block while the others remain fixed. In this setting,

Tseng [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF] proved the convergence of BCR to minimize non-smooth optimization problems of the form (2). Although we adopted this strategy to tackle sparse NMF problems in [START_REF] Rapin | NMF with sparse regularizations in transformed domains[END_REF], BCR requires an exact minimization for one block at each iteration, which generally leads to a high computational cost. We therefore opted for Proximal Alternating Linearized Minimization (PALM), which was introduced in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. It rather performs a single proximal gradient descent step for each coordinate at each iteration. Consequently, the PALM algorithm is generally much faster than BCR and its convergence to a stationary point of the multi-convex problem is guaranteed under mild conditions. In the framework of the proposed bGMCA algorithm, a PALMbased algorithm requires minimizing at each iteration eq. ( 2) over blocks of size 1 r n and alternating between the update of some submatrices of A and S (these submatrices will be noted A I and S I ). This reads at iteration (k) as:

1 -Update of a submatrix of S using a fixed A:

S (k) I = prox γG(.) A (k-1) T I A (k-1) I 2   S (k-1) I - γ A (k-1) T I A (k-1) I 2 A (k-1) T I (A (k-1) S (k-1) -X)   (3)
2 -Update of a submatrix of A using a fixed S:

A (k) I = prox δJ (.) S (k) I S (k) T I 2   A (k-1) I - δ S (k) I S (k) T I 2 (A (k-1) S (k) -X)S (k) T I   (4)
In eq. ( 3) and ( 4), the operator prox f is the proximal operator of f (cf. Appendix and [17] [18]). The scalars γ and δ are the gradient path lengths. The . 2 norm is the matrix norm induced by the 2 norm for vectors. More specifically, if x is a vector and .

2 is the 2 norm for vectors, the . 2 induced matrix norm is defined as:

M 2 = sup x =0 Mx 2 x 2 (5)

Block choice

Several strategies for selecting at each iteration the block indices I have been investigated: i) Sequential : at each iteration, r sources are selected sequentially in a cyclic way; ii) Random: at each iteration, r indices in [1, n] are randomly chosen following a uniform distribution and the corresponding sources updated; iii) Random sequential : this strategy combines the sequential and the random choices to ensure that all sources are updated an equal number of times. In the experiments, random strategies tended to provide better results. Indeed, compared to a sequential choice, randomness is likely to make the algorithm more robust with respect to spurious local minima.

Since the results between the random strategy and the random sequential one are similar, the first was eventually selected.

Examined cases and corresponding proximal operators

In several practical examples, an explicit expression can be computed for the proximal operators. In the next, the following penalizations have been considered:

1 -Penalizations G for the sources S:

-1 sparsity constraint in some transformed domain: The sparsity constraint on S is enforced with a 1 -norm penalization: G(S) = Λ S (SΦ T S ) 1 , where the matrix Λ S contains regularization parameters and denotes the Hadamard product. Φ S is a transform into a domain in which S can be sparsely represented. In the following, Φ S will be supposed to be orthogonal. The proximal operator for G in ( 3) is then explicit and corresponds to the softthresholding operator with threshold Λ S , which we shall denote S Λ S (.) (cf. Appendix). Using γ = 1 and assuming Φ S orthogonal, the update is then:

S (k) I = S Λ S   S (k-1) I Φ S T - 1 A (k-1) I A (k-1) T I 2 A (k-1) T I (A (k-1) S (k-1) -X)Φ S T   Φ S (6)
-Non-negativity in the direct domain and 1 sparsity constraint in some transformed domain: due to the non-negativity constraint, all coefficients in S must be non-negative in the direct domain in addition to the sparsity constraint in a transformed domain Φ S . It can be formulated as

G(S) = Λ S SΦ S T 1 + ι {∀j,k;S[j,k]≥0} (S)
where ι U is the indicator function of the set U . The difficulty is to enforce at the same time two constraints in two different domains, since the proximal operator of G is not explicit. It can either be roughly approximated by composing the proximal operators of the individual penalizations to produce a cheap update or computed accurately using the Generalized Forward-Backward splitting algorithm [START_REF] Raguet | A generalized forward-backward splitting[END_REF].

-Penalizations J for the mixing matrix A:

-Oblique constraint: to avoid obtaining degenerated A and S matrices ( A → ∞ and S → 0), the columns of A are constrained to be in the 2 ball, i.e. ∀j ∈ [1, n], A j 2 1. More specifically, J can be written as J (A) = ι {∀i; A i 2 2 ≤1} (A). Following this constraint, the proximal operator for J in eq. ( 4) is explicit and can be shown to be the projection Π . 2 (cf. Appendix) on the 2 unit ball of each column of the input. The update (4) of A I becomes:

A (k) I = Π . 2 1   A (k-1) I - 1 S (k) I S (k) T I 2 (A (k-1) S (k) -X)S (k) T I   (7)
-Non-negativity and oblique constraint: Adding the non-negativity constraint on A reads:

J (A) = ι ∀i; A i 2 2 ≤1 (A) + ι ∀i,j;A[i,j]≥0 (A).
The proximal operator can be shown to be the composition of the proximal operator corresponding to non-negativity followed by Π . 2 1 . The proximal operator corresponding to non-negativity is the projection Π K + (cf. Appendix) on the positive orthant K + .

The update is then:

A (k) I = Π . 2 1   Π K +   A (k-1) I - 1 S (k) I S (k) T I 2 (A (k-1) S (k) -X)S (k) T I     (8)

Minimization: introduction of a warm-up stage

While being provably convergent to a stationary point of (2), the above PALM-based algorithm suffers from a lack of robustness with regards to a bad initialization, which makes it more prone to be trapped in spurious local minima. Moreover, it is quite difficult to automatically tune the thresholds Λ so that it yields reasonable results. On the other hand, algorithms based on GMCA [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF][START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF] have been shown to be robust to initialization. Furthermore, in this framework, fixing the parameters Λ can be done in an automatic manner. However, GMCA-like algorithms are based on heuristics, which preclude provable convergence to a minimum of (2).

The proposed strategy consists in combining the best of both approaches to build a two-stage minimization procedure (cf. Algorithm 1): i) a warm-up stage building upon the GMCA algorithm to provide a fast and reliable first guess, and ii) a refinement stage based on the above PALM-based algorithm that provably yields a minimizer of (2). Moreover, the thresholds Λ in the refinement stage will be naturally derived from the first stage. Based on the GMCA algorithm [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF][START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF], the warm-up stage is summarized below: 0 -Initialize the algorithm with random A. For each iteration (k):

1 -The sources are first updated assuming a fixed A. A submatrix S I is however now updated instead of S. This is performed using a projected least square solution:

S (k) I = prox G(.) (A (k-1) † I R I ) (9)
where: R I is the residual term defined by R I = X -A (k )

I C S (k )
I C (with I C the indices of the sources outside the block), which is the part of X to be explained by the sources in the current block

I . A (k) † I is the pseudo-inverse of A (k)
I , the estimate of A I at iteration (k).

2 -The mixing sub-matrix A I is similarly updated with a fixed S:

A (k) I = prox J (.) (R I S (k) † I ) ( 10
)
The warm-up stage stops after a given number of iterations. Since the penalizations are the same as in the refinement stage, the proximal operators can be computed with the formulae described previously, depending on the implemented constraints. For S, eq. ( 6) can be used to enforce sparsity. To enforce non-negativity and sparsity in some transformed domain, the cheap update described in section 2.2.1 consisting in composing the proximal operators of the individual penalizations can be used. For A, equations ( 7) and ( 8) can be used depending on the implemented constraint.

Heuristics for the warm-up stage

In the spirit of GMCA, the bGMCA algorithm exploits heuristics to make the separation process more robust to initialization, which mainly consists in making use of a decreasing thresholding strategy. In brief, the entries of the threshold matrix Λ first start with large values and then decrease along the iterations towards final values that only depend on the noise level. This stategy has been shown to significantly improve the performances of the separation process [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF][START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF] as it provides: i) a better unmixing, ii) an increased robustness to noise, and iii) an increased robustness to spurious local minima.

In the bGMCA algorithm, this strategy is deployed by first identifying the coefficients of each source in I that are not statistically consistent with noise.

Assuming that each source is contaminated with a Gaussian noise with standard deviation σ, this is performed by retaining only the entries whose amplitude is larger than τ σ, where τ ∈ [START_REF] Biswal | Blind source separation of multiple signal sources of fMRI data sets using independent component analysis[END_REF][START_REF] Nuzillard | Application of blind source separation to 1-D and 2-D nuclear magnetic resonance spectroscopy[END_REF]. In practice, the noise standard deviation is estimated empirically using the Median Absolute Deviation (MAD)

estimator. For each source in I, the actual threshold at iteration k is fixed based on a given percentile of the available coefficients with the largest amplitudes. Decreasing the threshold at each iteration is then performed by linearly increasing the percentage of retained coefficients at each iteration:

Percentage = k iterations × 100.

Convergence

The bGMCA algorithm combines sequentially the above warm-up stage and the PALM-based refinement stage. Equipped with the decreasing thresholding strategy, it cannot be proved that the warm-up stage neither converges to a stationary point of eq. ( 2) nor converges at all. In practice, after consecutive iterates, the warm-up stage tends to stabilize. However, it plays a key role to provide a reasonable starting point, as well as threshold values Λ for the refinement procedure. In the refinement stage, the thresholds are computed from the matrices estimated in the warm-up and fixed for the whole refinement step. Based on the PALM algorithm, and with these fixed thresholds, the refinement stage converges to a stationary point of eq. ( 2).

The convergence is also guaranteed with the proposed block-based strategy, as long as the blocks are updated following an essentially cyclic rule [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF] or even if they are chosen randomly and updated one by one [START_REF] Patrascu | Efficient random coordinate descent algorithms for large-scale structured nonconvex optimization[END_REF].

Required number of iterations

Intuitively, the required number of iterations should be inversely proportional to r, since only r sources are updated at each iteration, requiring n/r times the number of iterations needed by an algorithm using the full matrices. As will be emphasized later on, the number of required iterations will be smaller than expected, which reduces the computation time.

In the refinement stage, the stopping criterion is based on the angular distance for each column of A, i.e. the angle between the current column and that of the previous iteration. Then, the mean over all the columns is taken:

∆ = j∈[1,n] A (k) j A (k-1) j 1 n (11)
The stopping criterion itself is then a threshold τ used to stop the algorithm when ∆ > τ . In addition, we also fixed a maximal number of iterations.

Numerical experiments on simulated data

In this part, we present our results on simulated data. The goal is to show and to explain on simple data how bGMCA works.

Experimental protocol

The simulated data were generated in the following way: The number of observations m is taken equal to the number of sources:

I = prox γG(.) A (k-1) T I A (k-1) I 2   S (k-1) I - γ A (k-1) T I A (k-1) I 2 A (k-1) T I (A (k-1) S (k-1) -X)   A (k) I = prox δJ (.) S (k) I S (k) T I 2   A (k-1) I - δ S (k) I S (k) T I 2 (A (k-1) S (k) -X)S (k) T I   ∆ = j∈[1,n] A (k) j A (k-1) j 1 n k = k + 1 end while return A, S 1
m = n.
In this first simulation, no noise is added. The algorithm was launched with 10, 000 iterations. It has to be emphasized that since neither A nor S are non-negative, the corresponding proximal operators we used did not enforce non-negativy. Thus, we used soft-thresholding for S and the oblique constraint for A according to section 2.2.1.

To measure the accuracy of the separation, we followed the definition in [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF] to use a global criterion on A:

C A = median(|PA † A * | -I d ),

Modeling block minimization

In this section, a simple model is introduced to describe the behavior of the bGMCA algorithm. As described in section 2.2, updating a given block

R I = X -A I C S I C = A * I S * I + E + N (12)
A way to further describe the structure of E is to decompose the S matrix in the true matrix plus an error: S I = S * I + I and S I C = S * I C + I C , where S is the estimated matrix, and is the error on S * . Assuming that the errors are small and neglecting the second-order terms, the residual R I can now be written as:

R I = X -A I C S I C = A * I S * I + A * I C S * I C -A I C S * I C -A I C I C + N (13)
This implies that:

E = (A * I C -A I C )S * I C -A I C I C (14)
Equation ( 14) highlights two terms. The first term can be qualified as interferences in that it comes from a leakage of the true sources that are outside the currently updated block. This term vanishes when A I C is perfectly estimated. The second term corresponds to interferences as well as artefacts.

It originates indeed from the error on the sources outside the block I. The artefacts are the errors on the sources induced by the soft thresholding corresponding to the 1 -norm.

Equation ( 14) also allows us to understand how the choice of a given block size r n will impact the separation process:

-Updating small-size blocks can be recast as a small-size source separation problem where the actual number of sources is equal to r. The residual of the sources that are not part of the block I then plays the role of extra noise. As testified by Fig. 1, updating small-size block problems should be easier to tackle.

-Small-size blocks should also yield larger errors E. It is intuitively due to the fact that many potentially badly estimated sources in I C are used for the estimation of A I and S I through the residual, deteriorating this estimation. It can be explained in more details using equation ( 14

Experiment

In this section, we investigate the behavior of the proposed block-based GMCA algorithm with respect to various parameters such as the block size, the number of sources, the conditioning of the mixing matrix and the sparsity level of the sources.

Study of the impact of r and n

In this subsection, bGMCA is evaluated for different numbers of sources n = 20, 50, 100. Each time the block sizes vary in the range 1 ≤ r ≤ n. In this experiment and to complete the description of section 3.1, the parameters for the matrices generation were: p = 0.1, t = 1, 000, C d = 1, m = n, with a Bernoulli-Gaussian distribution for the sources. These results are displayed in Fig. 2a. Interestingly, three different regimes characterize the behavior of the bGMCA algorithm:

-For intermediate and relatively large block sizes (typically r > 5 and r < n -5): we first observe that after an initial deterioration around r = 5 , the separation quality does not vary significantly for increasing block sizes. A degradation of several dB can then be observed for r close to n. In all this part of the curve, the error term E is composed of residuals of sparse sources, and thus E will be rather sparse when the block size is large. Based on the MAD, the thresholds are set according to dense and not to sparse noise. Consequently the automatic thresholding strategy of the bGMCA algorithm will not be sensitive to the estimation errors.

-A very prominent peak can be observed when the block size is of the order of 3. Interestingly, the maximum yields a mixing matrix criterion of about 10 -16 , which means that perfect separation is reached up to numerical errors. This value of 160 dB is at least 80 dB larger than in the standard case r = n, for which the values for the different n are all below 80 dB. In this regime, error propagation is composed of the mixture of a larger number of sparse sources, which eventually entails a densely distributed contribution that can be measured by the MAD-based thresholding procedure. Therefore, the threshold used to estimate the sources is able to filter out both the noise and the estimation errors. Moreover, r = 5 is quite small compared to n. Following the modeling introduced in section 3.2, small block sizes can be recast as a sequence of low-dimensional blind source separation problems, which are simpler to solve.

-For small block sizes (typically r < 4), the separation quality is deteriorated when the block size decreases, especially for large n values. In this regime, the level of estimation error E becomes large, which entails large values for the thresholds Λ. Consequently, the bias induced by the soft-thresholding operator increases, which eventually hampers the performance quality. Furthermore, for a fixed block size r, E increases with the number of sources n, making this phenomenon more pronounced for higher n values.

Condition number of the mixing matrix

In this section, we investigate the role played by the conditioning of the mixing matrix on the performances of the bGMCA algorithm. Fig. 2b displays the empirical results for several condition numbers C d of the A matrix.

There are n = 50 sources generated in the same way as in the previous experiment: with a Bernoulli-Gaussian distribution and p = 0.1, t = 1, 000.

One can observe that when C d increases, the peak present for r close to 5 tends to be flattened, which is probably due to higher projection errors. At some iteration k, the sources are estimated by projecting X -A I c S I c onto the subspace spanned by A I . In the orthogonal case, the projection error is low since A I c and A I are close to orthogonality at the solution. However, this error increases with the condition number C d .

Sparsity level p

In this section, the impact of the sparsity level of the sources is investigated. The sources are still following a Bernoulli-Gaussian distribution.

The parameters are: n = 50, t = 1, 000, C d = 1. As featured in Figure 3, the separation performances at the maximum value decrease slightly with larger p, while a slow shift of the transition between the small/large block size regimes towards larger block sizes operates. Furthermore, the results tend to deteriorate quickly for small block sizes (r < 4). Indeed, owing to the model of subsection 3.2, the contribution of S * I C and I C in the error term [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF] increases with p, this effect being even more important for small r (which could also explain the shift of the peak for p = 0.3, by a deterioration of the results at its beginning, r = 3). When p increases, the sources in S I become denser. Instead of being mainly sensitive to the noise and E, the MAD-based thresholding tends to be perturbated by S I , resulting in more artefacts, which eventually hampers the separation performances. This effect increases when the sparsity level of the sources decreases. Beyond improving the separation performances, the use of small block sizes decreases the computational cost of each iteration of the bGMCA algorithm. Since it is iterative, the final running time will depend on both the complexity of each iteration and of the number of iterations. In this part, we focus only on the warm-up stage, which is empirically the most computationally expensive stage. Each iteration of the warm-up stage can be decomposed into the following elementary steps: i) a residual term is computed with a complexity of O(mtr), where m is the number of observations, t the number of samples and r the block size; ii) the pseudo-inverse is performed with the singular value decomposition of a r × r matrix, which yield an overall complexity of O(r 3 + r 2 m + m 2 r); iii) the thresholding-strategy first requires the evaluation of the threshold values, which has a complexity of rt; iv) then the soft-thresholding step which has complexity O(rt); and v) updating A is finally performed using a conjugate gradient algorithm, whose complexity is known to depend on the number of non-zero entries in S and on the condition of this matrix C d (S). An upperbound for this complexity is thus O(rt C d (S)). The final estimate of the complexity of a single iteration is finally given by:

Complexity and computation time

r[mt + rm + m 2 + r 2 + t C d (S)] (15)
With C d (S) the conditioning number of S. Thus, both the r factor and the behavior in r 3 show that small r values will lower the computational budget of each iteration. We further assess the actual number of iterations required by the warm-up stage to yield a good initialization. To this end, the following experiment has been conducted:

1. First, the algorithm is launched with a large number of iterations (e.g. 10000) to give a good initialization for the A and S matrices. The corresponding value of C A is saved and called C * A .

2. Using the same initial conditions, the warm-up stage is re-launched and stops when the mixing matrix criterion reaches 1.05 × C * A (i.e. 5% of the "optimal"initialization for a given setting).

The number of iterations needed to reach the 5% accuracy is reported in Fig. 4. Intuitively, one would expect that when the block size decreases, the required number of iterations should increase by about n/r to keep the number of updates per source constant. This trend is displayed with the straight curve of Fig. 4. Interestingly, Fig. 4 shows that the actual number of iterations to reach the 5% accuracy criterion almost does not vary with r.

Consequently, on top of leading to computationally cheaper iterations, using small block sizes does not require more iterations for the warm-up stage to give a good initialization. Therefore, the use of blocks allows a huge decrease of the computational cost of the warm-up stage and thus of sparse BSS.

Experiment using realistic sources

Context

The goal of this part is to evaluate the behavior of bGMCA and show its efficiency in a more realistic setting. Our data come from a simulated LC -1 H NMR (Liquid Chromatography -1 H Nuclear Magnetic Resonance) experiment. The objective of such a experiment is to identify each of the chemicals compounds present in a fluid, as well as their concentrations. The LC -1 H NMR experiment enables a first physical imperfect separation during which the fluid goes through a chromatography column and its chemicals are separated according to their speeds (which themselves depend on their physical properties). Then, the spectrum of the output of the column is measured at a given time frequency. These measurements of the spectra at different times can be used to feed a bGMCA algorithm to refine the imperfect physical separation.

The fluids on which we worked could for instance correspond to drinks. The goal of bGMCA is then to identify the spectra of each compound (e.g. caffein, saccharose, menthone...) and the mixing coefficients (which are proportional to their concentrations) from the LC -1 H NMR data. BSS has already been successfully applied [START_REF] Toumi | Effective processing of pulse field gradient NMR of mixtures by blind source separation[END_REF] to similar problems but generally with lower number of sources n.

The sources (40 sources with each 10, 000 samples) are composed of elementary sparse non-negative theoretical spectra of chemical compounds taken from the SDBS database 1 , which are further convolved with a Laplacian having a width of 3 samples to simulate a given spectral resolution. Therefore, each convolved source becomes an approximately sparse non-negative row of S. The mixing matrix A of size (m,n) = (320,40) is composed of Gaussians (see Fig. 5), the objective being to have a matrix that could be consistent with the first imperfect physical separation.

Experiments

There are two main differences with the previous experiments of section 3: i) the sources are sparse in the undecimated wavelet domain Φ S , which is These results show that non-negativity yields a huge improvement for all block sizes r, which is expected since the problem is more constrained. This is probably due to the fact that all the small negative coefficients are set to 0, thus artificially allowing lower thresholds and therefore less artefacts. This is especially advantageous in the present context with very low noise2 (the Signal to Noise Ratio -SNR -has a value of 120 dB) where the thresholds do not need to be high to remove noise. Furthermore, the separation quality tends to be constant for r ≥ 10. In this particular setting, non-negativity helps curing the failure of sparse BSS when large blocks are used. However, using smaller block sizes still allows reducing the computation cost while preserving the separation quality. The bGMCA with non-negativity also compares favorably with respect to other tested standard BSS methods (cf. Section 1 for more details), yielding better results for all values of r. In particular, it is always better than HALS, which also uses non-negativity. As an illustration, a single original source is displayed in the right panel of Fig. 6 after its convolution with a Laplacian.

Its estimation using bGMCA with a non-negativity constraint is plotted in dashed line on the same graph, showing the high separation quality because of the nearly perfect overlap between the two curves. Both sources are drawn in the direct domain.

The robustness of the bGMCA algorithm with respect to additive Gaussian noise has further been tested. Fig. 7 reports the evolution of the mixing matrix criterion for varying values of the signal-to-noise ratio. It can be observed that bGMCA yields the best performances for all values of SNR.

Although it seems to particularly benefit from high SNR compared to HALS and EFICA, it still yields better results than the other algorithms for low SNR despite the small block size used (r = 10), which could have been particularly prone to error propagations.

Conclusion

While being central in numerous applications, tackling sparse BSS problems when the number of sources is large is highly challenging. In this article, we describe the block-GMCA algorithm, which is specifically tailored to solve sparse BSS in the large-scale regime. In this setting, the minimiza- All the numerical comparisons conducted show that bGMCA performs at least as well as standard sparse BSS on mixtures of a high number of sources and most of the experiments even show dramatically enhanced separation performances. As a byproduct, the proposed block-based strategy yields a significant decrease of the computational cost of the separation process.

Figure 1 :

 1 Figure 1: Evolution of the mixing matrix criterion (whose computation is detailed in sec. 3.1) of four standard BSS algorithms for an increasing n. For comparison, the results of the proposed bGM CA algorithm is presented, showing that its use allows for the good results of GMCA for low n (around 160 dB for n = 3) to persist for n < 50 and to stay much better than GMCA for n > 50. The experiment was conducted using exactly sparse sources S, with 10% non-zero coefficients, the other coefficients having a Gaussian amplitude. The mixing matrix A was taken to be orthogonal. Both A and S were generated randomly, the experiments being done 25 times and the median used to draw the figure.

  for 0 k < n max do Choose a set of indices I Estimation of S with a fixed A: S (k)I = prox G(.) (A (k-1) † I R I )Estimation of A with a fixed S: A (k)I = prox J (.) (R I S (k) † I ) Choice of a new threshold Λ (k)heuristic -see section 2.2.3 end for Refinement step while ∆ > τ and k < n max do Choose a set of indices I S (k)

- 2 -

 2 Source matrix S: the sources are sparse in the sample domain without requiring any transform (the results would however be identical for any source sparse in an orthogonal representation). The sources in S are exactly sparse and drawn randomly according to a Bernoulli-Gaussian distribution: among the t samples (t = 1, 000), a proportion p (called sparsity degree-unless specified, p = 0.1) of the samples is taken non-zero, with an amplitude drawn according to a standard normal distribution. Mixing matrix A: the mixing matrix is drawn randomly according to a standard normal distribution and modified to have unit columns and a given condition number C d (unless specified, C d = 1).

  where A * is the true mixing matrix and A is the solution given by the algorithm, corrected through P for the permutation and scale factors indeterminacies. I d is the identity matrix. This criterion quantifies the quality of the estimation of the mixing directions, that is the columns of A. If they are perfectly estimated, |PA † A * | is equal to I d and C A = 0. The data matrices being drawn randomly, each experiment was performed several times (typically 25 times) and the median of -10 log(C A ) over the experiments will be displayed. The logarithm is used to simplify the reading of the plots despite the high dynamics.

  is performed at each iteration from the residual R I = X -A I C S I C . If the estimation were perfect, the residual would be equal to the part of the data explained by the true sources in the current block indexed by I, which would read: R I = A * I S * I , A * and S * being the true matrices. It is nevertheless mandatory to take into account the noise N, as well as a variety of flaws in the estimation by adding a term E to model the estimation error. This entails:

  ): with more sources in I C , the energy of A I C , A * I C , S * I C and I C increases, yielding bigger error terms (A * I C -A I C )S * I C and -A I C I C . Therefore the errors E become higher, deteriorating the results.

  (a) Number of sources. (b) Condition number.

Figure 2 :

 2 Figure 2: Left: mixing matrix criterion as a function of r for different n. Right: mixing matrix criterion as a function of r for different C d .

Figure 3 :

 3 Figure 3: Mixing matrix criterion as a function of r for different sparsity degrees.

Figure 4 :

 4 Figure 4: Right: number of iterations in logarithmic scale as a function of r.

  It is designed in two parts: the first columns have relatively spaced Gaussian means while the others have a larger overlap to simulate compounds for which the physical separation is less discriminative. More precisely, an index m ∈ [1, m] is chosen, with m > m/2 (typically, m = 0.75m ). A set of n/2 indices (m k ) k=1... n/2 is then uniformly chosen in [0, m] and another set of n/2 indices (m k ) k= n/2 ...n is chosen in [ m + 1, m]. Each column of A is then created as a Gaussian whose mean is m k . Monte-carlo simulations have been carried out by randomly assigning the sources and the mixing matrix columns. The median over the results of the different experiments will be displayed.

Figure 5 :

 5 Figure 5: Exemple of A matrix with 8 columns: the four first columns have spaced means, while the last ones are more correlated

Figure 6 :

 6 Figure 6: Left: mixing criterion on realistic sources, with and without a non-negativity constraint. Right: example of a retrieved source, which is almost perfectly superimposed on the true source, therefore showing the quality of the results.

Figure 7 :

 7 Figure 7: Mixing criterion on realistic sources, using a non-negative constraint with r = 10

National Institute of Advanced Industrial Science and Technology (AIST), Spectral database for organic compounds: http://sdbs.db.aist.go.jp

Depending on the instrumentation, high SNR values can be reached in such an experiment

Acknowledgement

This work is supported by the European Community through the grant LENA (ERC StG -contract no. 678282).

Appendix

Definition of proximal operators

The proximal operator of an extended-valued proper and lower semicontinuous convex function f : R n → (-∞, ∞] is defined as:

Definition of the soft thresholding operator

The soft thresholding operator S λ (.) is defined as:

Definition of the projection of the columns of a matrix M on the 2 ball

Definition of the projection of a matrix M on the positive orthant