L Le Coq

INFLUENCE ON PERMEABILITY OF THE STRUCTURAL PARAMETERS OF HETEROGENEOUS POROUS MEDIA

Keywords: Mercury porosimetry, image analysis, pore size, porosity, tortuosity

Predicting the macroscopic properties of porous media used in treatment processes is a complex task regarding 3-D structures at micro-, meso-and macro-levels. Currently, information is scarce concerning the influence, at a microscopic level, of the 3-D structure of fibrous media on the physical laws governing their macroscopic behaviour. Nevertheless, the relationship between macroscopic properties (pressure drop, treatment efficiency) and microstructure can be assessed thanks to suitable structure modelling theories. In this context, the present study proposes and compares different methods (mercury porosimetry and image analysis) for the structure characterization at a microscopic level of filtering fibrous media, such as nonwoven and woven fabrics. The results obtained show a porous structure gradient in the thickness of the nonwoven media studied in terms of porosity, pore size and tortuosity factor. Moreover, the influence on structural parameters of media compression, when submitted to friction forces exerted by flow during filtration tests, is established. A model for the determination of multi-level pore size distributions from mercury porosimetry data is proposed. The "equivalent pore" model is used to estimate the tortuosity factor. The influence of measured structural parameters on fibrous media permeability is studied in a classical model for flow through fibrous media.

INTRODUCTION

Porous materials are widely used in treatment processes for removing pollution from air and water for environmental purposes or for fluid purification in industry. In most cases, arrangements of porous material present different levels of pore size from macropores, ranging from millimetres to micrometres, to micropores, which can include pore sizes from micrometres to nanometres according to the porous material considered. Macropores are formed by porous material units that are either packed in a reactor, constituting an unconsolidated porous medium, or bonded together in the form of a nonwoven or woven membrane (consolidated porous medium). Macropores are also called inter-pores as they are involved in inter-material porosity. Usually, mesoand micropores represent, respectively, the pores at the surface and inside porous material, thus they develop the intra-material porosity.

For treatment process optimization and design purposes, it is a complex task to consider the multi-levels of pore size for the estimation of process macroscopic behaviour with regard to flow, mass and thermal transfer. The state of the art concerning the influence of hierarchical porous material on the physical laws governing the process in which the porous structure is involved mostly deals with homogenization techniques [START_REF] Cushman | Theory and Applications of Transport in Porous Media -The Physics of fluids in hierarchical Porous Media: Angstroms to Miles[END_REF][START_REF] Kfoury | Upscaling fractured heterogeneous media: Permeability and mass exchange coefficient[END_REF].

In this context, an accurate characterization of structural parameters at the local level is of interest and some studies related to paper are available in the literature [START_REF] Holmstad | Quantification of the three-dimensional paper structure: Methods and potential[END_REF][START_REF] Silvy | The paper and its unsuspected dimensions[END_REF].

The present study deals with the characterization of the structural parameters of fibrous media exhibiting a pore size gradient in their thickness. The two characterization methods compared (mercury porosimetry [START_REF] Giesche | Mercury porosimetry: A general (practical) overview[END_REF] and image analysis [START_REF] Coq | Characteristics of fixed-beds packed with anisotropic particles -Use of image analysis[END_REF]) are developed to obtain different information about the porous structure of eight papers made of a glass fibre blend.

On the one hand, mercury porosimetry measures porosity and pore hydraulic diameter distribution at the different pore scales (macro-to micropore) and for different compression levels of the media. Moreover, the tortuosity factor in the thickness of the media is deduced from mercury porosimetry measurements thanks to the "equivalent pore" model [START_REF] Silvy | La transformation conforme du pore équivalent, méthode d'homogénéisation de la texture des milieux polyphasés[END_REF]. On the other hand, image analysis of media cross-sections is performed to determine porosity, mean pore size and tortuosity factors for different sub-layers of the media (upstream, centre and downstream layers) [START_REF] Coq | Characteristics of fixed-beds packed with anisotropic particles -Use of image analysis[END_REF]. From structural parameter data, permeability coefficients are calculated and compared with those deduced from permeametry measurements.

MATERIALS AND METHODS

Materials

Operating conditions of fibrous media production

Eight glass fibre nonwovens made by a papermaking process were studied. Fibrous media are mixtures of Glass Fibres (GF) of various diameters and lengths fabricated by using a standard hand sheet machine [START_REF] Paulapuro | Papermaking Part 1: Stock preparation and wet end In[END_REF]. Three grades of glass fibre were used in this study corresponding to mean diameters (d) of 0.6, 2.7 and 6 µm, for fibre lengths (L) ranging between 500 d < L < 4000 d. The fibre composition of each medium is expressed considering the weight percentage of each fibre grade dispersed in the fibrous admixture.

Stock glass fibres were dispersed in water at an acidic pH of 2.5, at a concentration of 10 g l -1 using a propeller mixer.

The dispersion operation was performed at 2500 to 3000 rd min -1 for 10 to 40 min according to fibre diameter and length, to avoid flocking of fibre. The required energy for a good dispersion increases with the length-to-diameter ratio for fibres of diameter greater than 1 µm [START_REF] Huber | Suspension crowding for a general fiber length distribution: application to flocculation of mixtures of short and long papermaking fibres[END_REF]. For fibres of diameter less than 1 µm, the required energy increases with decreasing diameter. The acidic pH of the suspension improves fibre repelling [START_REF] White | Wet formed non-woven webs from high performance fibres[END_REF].

The fibrous suspension was diluted in the mixing vessel of the hand sheet machine at a concentration of 0.5 g l -1 and mixed with a chemical binder (5 to 10 % in weight) for 30 s.

The binder ensures the mechanical cohesion of the fibrous web. The suspension was then percolated through a forming fabric thanks to a suction system providing a vacuum pressure of 0.1 10 5 Pa. The fibrous mat obtained was drained for 0.5 to 4 min according to mat permeability and was slightly pressed (2000 Pa) onto blotting paper. A drying operation at 110 °C for 35 min produces strengthened fibrous media.

Considering the process used to make paper prototypes, the fibrous media obtained can be considered isotropic in the sheet plane but anisotropic in the thickness because of a fibrous composition gradient due to filtration on the forming fabric.

Fibrous media obtained

In

Characterization of Structural Parameters

The structure of porous media can be characterized considering two scalar quantities, porosity and pore hydraulic diameter, and a vector quantity, the tortuosity factor [START_REF] Silvy | The paper and its unsuspected dimensions[END_REF][START_REF] Coq | Fibrous media plugging modelling for liquid filtration[END_REF]. In the present study, tortuosity is characterized according to the mean flow path direction, that is to say according to medium thickness.

The mean pore hydraulic diameter (<φ >) is defined as the ratio of the pore volume where fluid flows to the pore surface wetted by this fluid. It is commonly used to represent the characteristic dimension of a non-circular channel.

< φ > = 4 ε S v (i)
where S v (m -1 ) is the volume specific surface area of pores, which can be calculated as:

S v = (1-ε) fibre surface fibre volume (ii)
Two methods were developed to determine these structural parameters. The first one is based on postprocessing of mercury porosimetry results. This technique gives access to the structural parameters of the overall medium at different levels of compression. The second method is based on stereological analysis of images obtained from media cross-sections and can characterize local structural parameters [START_REF] Montillet | Image analysis: A useful tool for porous media characterization[END_REF]. Indeed, in the case of heterogeneous or anisotropic porous media for environmental engineering processes, considering the overall structural parameters could lead to design errors or to deficient treatment performances. So, for particle filtration, an accurate prediction of fibrous media clogging can be achieved by considering the pore size gradient in the medium thickness.

Mercury porosimetry

Mercury porosimetry measurements were performed using a Micromeritics porometer (pore size range: 150 to 0.4 µm) on unstrained and compressed samples of the studied media. For each sample, three reproducibility experiments were performed. The mercury volume penetrating the porous medium was measured at different pressures so that a porosity versus pore size curve can be obtained according to the Washburn theory [START_REF] Washburn | The dynamics of capillary flow[END_REF]. For data interpretation, a contact angle of 153 ° (glass -mercury) and a surface tension of 0.48 Nm _1 were considered.

The experimental data of porosity versus pore size were modelled considering the different levels of pore size distribution constituting a porous medium. In the case of the glass fibre prototypes, two levels of pore size can be distinguished: surface roughness and media internal pores. In the particular case of woven fabrics made of multi-fibre yarns, there are three levels of pore size: cloth surface roughness, inter-yarn pores due to yarn weaving and intra-yarn pores due to fibre twisting [START_REF] Benesse | Study of particle capture mechanisms in woven activated carbon filters: image analysis[END_REF].

Considering the maximum porosity (ε i max

) and the modal pore size (φ i mod ) of each level (i) of pore size distribution, the following equation can be proposed to model the experimental probability distribution of porosity versus pore size:

ε(φ j ) = P(φ < φ j ) = ε i max (φ j / φ i mod ) α i +1 ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ i=l k ∑ (iii)
α i is the slope attenuation coefficient at the inflection point of the cumulative porosity curve for the pore size distribution i.

It is representative of how narrow a mono-modal distribution is: the smaller the α i value, the narrower the pore size distribution.

Maximum porosity (ε i max ), modal pore hydraulic diameter (φ i mod

) and attenuation coefficient (α i ) are determined for each level (i) of pore size by an optimization algorithm.

The iterative algorithm is based on a Newton method used to solve a least squares problem.

From the intrusion volume versus pore size data, a mercury saturation curve can be plotted considering, on the one hand, the dimensionless saturation (s) defined as the ratio between the cumulative mercury volume at a given pore size and the total incremental mercury volume and, on the other hand, the dimensionless hydraulic diameter (φ 0 ) equal to the ratio between the considered pore size and the modal pore size.

Silvy [START_REF] Silvy | La transformation conforme du pore équivalent, méthode d'homogénéisation de la texture des milieux polyphasés[END_REF] models the porous structure using a homogenization approach that considers the pore/fibre interface at microscale as an ellipsoidal equivalent pore at macroscale characterized by its ellipticity coefficients b/a and c/a as well as by its surface area, equal to the measured pore/fibre interface at microscale. In order to estimate the flow through the homogenized porous model, he defines the flow direction that minimizes the pressure drops, as shown in Figure 1.

From the "equivalent pore" theory, the mercury saturation curve represents the progressive filling of the equivalent ellipsoid so as to maximize the hydraulic diameter at each step, as a result of the minimum loss of energy principle. In the case of a fibrous medium of isotropic structure in the paper plane and for c/a <1, the ellipticity coefficient (c/a) can be estimated by fitting the modelled saturation curve from equation (iv) with the experimental mercury saturation data.

φ° = 3 2 (1-s) 2/ 3 1 c/a 1+(c/a) 2 1 1-(c/a) 2 Ln 1 c/a 1-(c/a) 2 +1 ( ) ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ 1+ 1 c/a ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ 2 -1 ⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ (1-s) 2/3 (iv)
Figure 1.

Representation of the equivalent pore model [START_REF] Silvy | La transformation conforme du pore équivalent, méthode d'homogénéisation de la texture des milieux polyphasés[END_REF].

The tortuosity factor is then calculated from equation (v).

t z 2 = 1+ 1 c / a ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ 2 (v)

Image analysis

The discriminate pores and fibres. The grey level images were binarized using a mean square criteria algorithm for thresholding [START_REF] Otsu | An automatic threshold selection method based on discriminant and least square criteria[END_REF] so as to separate fibres from pores. The binarized images were then analysed using Visilog software in order to measure stereological parameters [START_REF] Underwood | Quantitative Stereology[END_REF][START_REF] Montillet | Image analysis: A useful tool for porous media characterization[END_REF].

Stereological analysis provides the secant number in the fibres per unit of measured length (N L , mm -1 ) and the total surface area of the pores per unit of measured surface area The stereological analysis was performed considering different viewing angles of the cross-section observed, from 0 to 180 ° and every 15 °, in order to access the structure anisotropy and to estimate the tortuosity factor (eq. v) of the porous path according to the "equivalent pore" model [START_REF] Silvy | La transformation conforme du pore équivalent, méthode d'homogénéisation de la texture des milieux polyphasés[END_REF].

The structural parameters of porous media are calculated from the stereological analysis and the equivalent pore approach as follows: ε = volume des pores volume total = A A (vi)

S V = 4 <N L > (vii) < φ > = 4 A A S v = A A < N L > (viii)
where c / a = N Lx N Lz

Estimation of permeability

The permeability coefficient (B, m 2 ) can be deduced from the above structural parameters considering equation (ix) [START_REF] Carman | The Flow of Gases in Porous Media[END_REF]. This permeability coefficient represents an intrinsic property of the porous media.

B = ε (< φ > /4) 2 2t z 2 = ε 3 2t z 2 S v 2 (ix)

RESULTS AND DISCUSSION

Structural Parameters

The overall structural parameters of the unstrained fibrous media are measured using mercury porosimetry and image analysis techniques. Local structural parameters are also obtained from stereological analysis of the three layers constituting the media (see Figure 2). Moreover, mercury porosimetry is used to characterize the influence of medium deformation, when submitted to pressure drop, on the modification of structural parameters.

Unstrained fibrous media

The characterization of structural parameters from mercury porosimetry and stereological analysis is presented in Table 2 for media not submitted to the action of any friction forces.

From mercury intrusion volume curves and applying equation (iii), two characteristic pore size distributions can be distinguished for all the fibrous media, as shown in Figure 3a for the 33% 0.6 + 67% 2.7 medium. The wider pore size distribution corresponds to the roughness of the fibrous media and represents pores located on the upstream and downstream surfaces of the media. The distribution of smaller pores characterizes the internal porous structure inside the medium thickness. Moreover, the two 100% 2.7 media exhibit similar mean pore sizes, which is in agreement with the fact that medium pore size mainly depends on fibre diameter.

Compression influence on media structure

In order to evaluate the influence of medium compression, when submitted to pressure drops generated by fluid flow through a porous structure, some mercury porosimetry measurements were performed on compressed media. Table 3 presents the variations in structural parameters and permeability coefficients for different applied pressures.

The relative variation in thickness (∆e/e) resulting from applied pressure is also reported in 3. The corresponding permeability coefficients are deduced from equation (ix) and from permeametry measurements (equation x).

A good agreement between the two permeability coefficients (when permeametry measurements are available)

is observed. This tends to validate the estimation of structural parameters from the "equivalent pore" model applied to the mercury porosimetry technique developed in the present study for compressed media. Moreover, the structural parameters obtained indicate that the relative variation in medium thickness is a result of a modification of internal structural parameters. Indeed, Table 3 shows a significant decrease in porosity and pore size values and a large increase in the tortuosity factor following an increase in medium deformation.

Considering that the range of pressures studied is of the same order of magnitude as pressure drops experienced by fibrous media during liquid filtration, the results obtained indicate that porous medium deformation must be considered for filtration optimization or prediction [START_REF] Coq | Fibrous media plugging modelling for liquid filtration[END_REF].

CONCLUSIONS

Stereological analysis is of interest to measure local structural parameters when porous media present a porosity and/or pore size gradient. The results obtained using this image analysis technique on anisotropic fibrous media show the same variation as structural parameters deduced from mercury porosimetry data post-processed by means of the "pore equivalent" model. Nevertheless, stereological analysis leads to overestimated structural parameters compared with mercury porosimetry. This is mainly due to the fact that mercury porosimetry data are post-processed to exclude media roughness information (which is not of interest for media permeability prediction) while image analysis data are not.

Permeability coefficients calculated from structural parameters deduced from mercury porosimetry experiments are in good agreement with those obtained from permeametry measurements for unstrained media as well as for compressed media. This can be explained by the major influence of the internal porous structure of media on flow behaviour through porous media.

Finally, a pore size distribution model has been developed to predict mercury porosimetry data. This can be used to quantify the major influence on structural parameters of fibrous media deformation when submitted to flow.

  use of image analysis provides accurate local variations in the radial and longitudinal structural properties, such as porosity and mean pore size. This technique is based on the consolidation of the porous structure by the use of a resin (araldite CY212), which fills the porous space of the fibrous medium. In order not to modify noticeably the porous structure, the fibrous sample and the resin were degassed separately until a low-vacuum was reached. Then, the resin was poured into the crucible containing the sample and atmospheric pressure was slowly reached to allow the slow percolation of the resin into the pores. Next, using a rotary microtome and polishing disks (from 60 to 0.25 µm paper grades), longitudinal or transversal slices of consolidated samples were obtained for Scanning Electron Microscope (SEM) observations using backscattered electron imaging at a magnification of 450 times in order to

(A A , mm 2 /

 2 mm 2 ). These two stereological parameters are determined from images of cross-sections x0z and x0y and for each of the layers constituting the fibrous medium, called upstream, centre and downstream respectively as shown in Figure2[START_REF] Coq | Characteristics of fixed-beds packed with anisotropic particles -Use of image analysis[END_REF]. The number of images analysed was chosen so as to reach 2 % precision for the intercept number measurements and 0.2 % for the pore surface area respectively, which corresponds to approximately 40 images and a surface area of 1 mm 2 according to the sample.

Figure 2 .

 2 Figure 2. Schematic view of consolidated sample slices.

Figure 3 .

 3 Figure 3. Cumulative porosity and saturation curves deduced from mercury porosimetry experiments for the medium 33% 0.6 + 67% 2.7.

Figure 4 .

 4 Figure 4. Permeability coefficients obtained from equation (ix) from mercury porosimetry and image analysis experiments: comparison with air permeametry measurements.

Figure 5

 5 Figure5presents the influence of the sheet formation process on structural parameters at different levels in the thickness of three studied media. Pore size and porosity gradients can be observed between the upstream and the downstream layer and are sharper for the 33% 0.6 + 67% 2.7 medium than for the 100% 2.7 samples. This can be explained by the fact that, during fibre filtration on the forming fabric, fibres of smaller diameter and length percolate through the fibre cake and concentrate in the downstream side of the medium (sieved by the forming fabric). In the case of 100% 2.7 samples, diameter and length dispersions are weak, resulting in a slight pore gradient. In contrast, for a blend of two different fibre compositions (mean diameters of 0.6 and 2.7 µm), characteristic dimension dispersions are significant and lead to a fibre composition gradient in the medium thickness.

Figure 5 .

 5 Figure 5.Porosity and pore hydraulic diameter measured from stereological analysis at the three levels in the thickness of three media.

Figure 6 .

 6 Figure 6. Differential porosity curves deduced from mercury porosimetry experiments: influence of strain (∆e/e) on the porous structure of the 33% 0.6 + 67% 2.7 medium.

2 L

 2 fibre length, mm N Li secant number per unit of measured length in the i direction, mm-1

1 P

 1 <N L > mean secant number per unit of measured length (mean value for all directions of plane), mm -probability distribution of porosity versus pore size W specific weight, g m -2S v specific surface area per unit of volume, m -1

Table 2

 2

shows good adequacy between the mean hydraulic diameters and porosity measured by the two methods. The results indicate the influence of fibrous composition on the structural parameters of the paper. Indeed, for the same mean fibre diameter, fibre length exhibits a significant influence on porous structure comparing data obtained for the 100% 2.7 -L = 0.3 mm and 100% 2.7 -L = 0.9 mm media. Porosity and hydraulic pore diameter increase with fibre length whereas media made of a blend of the two fibre lengths, 100% 2.7 furnished, leads to an intermediate porous structure. Considering the influence of fibre diameter on structural parameters, the two media made of one set of fibres (100% 0.6 and 100% 2.7 µm furnished) present a decrease in hydraulic pore diameter with decreasing mean fibre diameter. This result is consistent with the

Table 2 .

 2 Structural parameters and permeability coefficients deduced from mercury porosimetry and image analysis for unstrained fibrous media.

	Fibrous media		Mercury porosimetry			Image analysis
		ε	<φ>	t 2 z	B	ε	<φ>	t 2 z	B
			(µm)		.10 -12 (m 2 )		(µm)		.10 -12 (m 2 )
			eq. (i)	eq. (iv)	eq. (ix)	eq. (vi)	eq. (viii)	eq. (vi)	eq. (ix)
	100% 0.6	0.926	7	2.2	0.6
	33% 0.6 + 67% 2.7	0.910	10	2.7	1.1	0.916	13	1.9	2.5
	100% 2.7 furnished	0.930	25	2.4	7.6
	100% 2.7 -L=0.3 mm	0.925	22	2.1	6.7	0.925	24	1.8	9.2
	100% 2.7 -L=0.9 mm	0.936	28	2.8	8.2	0.939	28	2.1	11.0
	67% 0.6 + 33% 6	0.917	8	2.4	0.8
	50% 0.6 + 50% 6	0.920	12	2.4	1.7	0.912	18	2.0	4.6
	33% 0.6 + 67% 6	0.921	14	3	1.9

Table 3 .

 3 Structural parameters and permeability coefficients deduced from mercury porosimetry and permeametry for fibrous media under pressure.

					Mercury porosimetry		Permeametry
	Media	∆P	∆e/e	e	<φ>	T 2 z	B	B
		(10 5 Pa)			(µm)	.10 -12 (m 2 )	.10 -12 (m 2 )
		0	0	0.93	10	2.7	1.1	1.4
	33% 0.6	0.18	0.15	0.89	8.2	4.0	0.47	0.49
	+ 67% 2.7	0.53	0.30	0.86	7.6	4.7	0.33	0.37
		1.42	0.45	0.83	4.9	6.1	0.10
		0	0	0.94	28	2.8	8.2	10.0
	100% 2.7	0.14	0.15	0.92	17.0	4.1	2.0	2.1
	L = 0.9 mm	0.27	0.25	0.89	16.5	5.5	1.4
		0.80	0.45	0.86	16.0	41.6	0.17


LPC Pierre Courrieu
email: pierre.courrieu@univ-provence.fr

Fast solving of Weighted Pairing Least-Squares systems

Keywords: weighted pairing least-squares, generalized inverses, generalized Cholesky factor

This paper presents a generalization of the "weighted least-squares" (WLS), named "weighted pairing least-squares" (WPLS), which uses a rectangular weight matrix and is suitable for data alignment problems. Two fast solving methods, suitable for solving full rank systems as well as rank deficient systems, are studied. Computational experiments clearly show that the best method, in terms of speed, accuracy, and numerical stability, is based on a special {1, 2, 3}-inverse, whose computation reduces to a very simple generalization of the usual "Cholesky factorization-backward substitution" method for solving linear systems.

Introduction

In this paper, we consider weighted least-squares problems in the following generalized form:

E(V;X,Y,W )  i1 m1  w ij X i,: V  Y j,: 2 j1 m 2  , C  arg min V R n1n 2 E(V;X,Y,W ), (1)
where three matrices are given: X  R m1n1 , Y  R m 2 n 2 , and W  (w ij ) is a rectangular m1  m2 "weighted pairing" matrix whose coefficients are non-negative real numbers. X i,:

denotes the ith row of X , and Y j,: denotes the jth row of Y . Note that this generalization is not the same as the so-called "Generalized Least Squares" [START_REF] Yuan | Numerical methods for generalized least squares problems[END_REF]. In the special case where m1  m2  m and W is a diagonal matrix, the above problem clearly reduces to an ordinary weighted least-squares problem, that is:

C  arg min V R n1n 2
w ii X i,: V  Y i,:

2 i1 m   arg min V R n1n 2 W 1/ 2 XV  W 1/ 2 Y 2 . ( 2
)
In Section 2, we show that, in fact, every problem having the form (1) can be reduced to a problem having the form [START_REF] Bube | Hybrid l 1 /l 2 minimization with applications to tomography[END_REF]. In such problems, each equation of the least-squares system receives a specific weight that typically depends on some estimate of the reliability of the data used in that equation. The usual non-weighted case corresponds to W  I (identity matrix).

Ordinary weighted least-squares [START_REF] Bube | Hybrid l 1 /l 2 minimization with applications to tomography[END_REF] are commonly used to solve regression problems with noisy data [START_REF] Zhou | A weighted least squares method for scattered data fitting[END_REF], and in "iteratively re-weighted least-squares" procedures for computing robust regression statistics such as M-estimators [START_REF] Bube | Hybrid l 1 /l 2 minimization with applications to tomography[END_REF] [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF]. The generalization (1) is potentially relevant in "data alignment" problems, where there is no given one-to-one correspondence between X data points and Y data points (rows), but one has some non-negative "adequacy"

or "plausibility" measure for each possible data pair, which is represented by W . Data alignment is a hard to solve problem commonly encountered in image processing and pattern recognition [START_REF] Lai | Robust image matching under partial occlusion and spatially varying illumination change[END_REF]. In this paper, ||.|| denotes the Euclidean norm for vectors, and the Frobenius norm for matrices. As in Matlab, the notation "a:b" denotes an index interval of bounds "a"

and "b", and if the bounds are not specified (:), this corresponds to the whole index range.

Whenever W 1/ 2 X is of full column rank in [START_REF] Bube | Hybrid l 1 /l 2 minimization with applications to tomography[END_REF], the solution to this problem is unique and the normal equations lead to the well-known result:

C  (X'WX) 1 X'WY , (3)
where X' denotes the transpose of X (or the conjugate transpose in the complex case).

One of the fastest ways of numerically obtaining the factor (X'WX) 1 that appears in

(3) consists of computing a Cholesky factorization LL' of the positive definite Gram matrix X'WX , then one inverts the upper triangular factor L' by a simple backward substitution method, and one obtains (X'WX) 1  L' 1 L 1 . However, if W 1/ 2 X is not of full column rank, then the above method does not work because the matrix X'WX is singular, and in this case, the weighted least-squares system is said to be rank deficient. The solution of rank deficient systems requires more robust methods, which are also slower than the above mentioned, in general. Among the fastest methods, we can consider those based on the use of suitable generalized inverses, such as the popular Moore-Penrose inverse [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF][4] [START_REF] Rakha | On the Moore-Penrose generalized inverse matrix[END_REF]. A solution to (2) is then:

C  (W 1/ 2 X) † W 1/ 2 Y  (X'WX) † X'WY , (4)
where the Moore-Penrose inverse is known to provide the least-squares solution C whose each column has the minimum Euclidean norm ( [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF], p. 109).

However, one must note that the solution of least-squares problems does not specifically require the use of the Moore-Penrose inverse, and that other types of generalized inverses, such as {1, 3}-inverses whose numerical computation is possibly faster, can as well be used. According to ( [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF], pp. 104-105), one has always a solution to (2) with:

C  (W 1/ 2 X) (1, 3) W 1/ 2 Y , (5)
where A [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF][START_REF] Courrieu | Straight monotonic embedding of data sets in Euclidean spaces[END_REF] denotes any {1, 3}-inverse of the matrix A (see Section 3.2).

In fact, the problem of the computational cost is crucial in many practical applications, where one must repeatedly solve large least-squares systems. On the other hand, most practical problems lead to full rank systems that could be solved fast using (3), however, rank deficient systems can occasionally appear, and it is commonly not acceptable to obtain a "fatal error" diagnostic at run time. Thus, in order to optimize the performance of applications, we present in Section 3.2 a quite fast solution of type (4), and in Section 3.3 a solution of type (5) whose computational cost is similar to that of (3), which has the advantage of being fast while providing a suitable least-squares solution in all cases, even if the system is rank deficient. These solutions apply to problem (2) and to problem (1) as well.

The weighted pairing least-squares problem

In this section, we consider the generalization of the weighted least-squares (WLS) problem stated in [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF], which we refer to as the "weighted pairing least-squares" (WPLS) problem.

Theorem 1. Every WPLS problem of type (1) reduces to a WLS problem of type (2) since:

arg min

V R n1n 2 i1 m1  w ij X i,: V  Y j,: 2 j1 m 2   arg min V R n1n 2 h ii X i,: V  Z i,: 2 i1 m1  ,
where:

H  (h ij ) is a diagonal matrix with diagonal coefficients h ii  w ij , 1 i  m1 j1 m 2  , Z  H † WY , H †  (h ij † ) is the Moore-Penrose inverse of H , with h ii †  1/h ii if h ii  0 , and h ij †  0 if h ij  0.
Proof.

Set d ik  w ij y jk 2 j1 m 2   h ii † ( w ij y jk ) 2 , 1 i  m1, 1 k  n2 j1 m 2  . ( 6
)
Then one has:

5/21 (h ii X i,: V  Z i,: 2 i1 m1   d ik k1 n 2  )  h ii (X i,: k1 n 2  i1 m1  V :,k  h ii † w ij y jk ) 2  d ik j1 m 2   (h ii 1/ 2 X i,: k1 n 2  i1 m1  V :,k  h ii 1/ 2 h ii † w ij y jk ) 2  d ik j1 m 2   h ii (X i,: k1 n 2  i1 m1  V :,k ) 2  2h ii h ii † (X i,: V :,k ) w ij y jk  h ii h ii †2 ( w ij y jk j1 m 2  ) 2  d ik j1 m 2   h ii (X i,: k1 n 2  i1 m1  V :,k ) 2  2(X i,: V :,k ) w ij y jk  h ii † ( w ij y jk j1 m 2  ) 2  d ik j1 m 2   ( w ij j1 m 2  (X i,: k1 n 2  i1 m1  V :,k ) 2 )  2( w ij y jk (X i,: V :,k ))  ( w ij y jk 2 j1 m 2  ) j1 m 2   w ij j1 m 2  ((X i,: k1 n 2  i1 m1  V :,k ) 2  2(X i,: V :,k )y jk  y jk 2 )  i1 m1  w ij X i,: V  Y j,: 2 j1 m 2  .
Noting that the additional terms ( d ik ) given by ( 6) are independent of V , one obtains Theorem 1.☐ Corollary 1.

(i) If H 1/ 2 X is of full column rank, then (1) has the unique solution: C  (X' HX) 1 X' HZ  (X' HX) 1 X'WY .
(ii) No matter H 1/ 2 X is not of full column rank, (1) has the minimum norm solution:

C  (H 1/ 2 X) † H 1/ 2 Z  (X' HX) † X'WY .
(iii) No matter H 1/ 2 X is not of full column rank, (1) has all solutions of the form:

C  (H 1/ 2 X) (1, 3) H 1/ 2 Z , 6/21
where H and Z are defined as in Theorem 1.

Proof. This directly follows from Theorem 1 and equation ( 3) for (i), equation ( 4) for (ii), and equation ( 5) for (iii).☐

Fast solutions based on generalized inverses

Generalized Cholesky factors

Several generalizations of the Cholesky factorization can be found in the literature. A well-known generalized Cholesky factorization for solving the so-called "augmented linear systems" is available in [START_REF] Wang | Perturbation analysis for the generalized Cholesky factorization[END_REF] and [START_REF] Zhao | The generalized Cholesky factorization method for saddle point problems[END_REF]. Another type of generalization of the Cholesky factorization has been proposed in [START_REF] Courrieu | Straight monotonic embedding of data sets in Euclidean spaces[END_REF], and this approach has been successfully used to define a fast numerical method for computing the Moore-Penrose inverse [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF]. The fundamental result for the generalized Cholesky factorization is:

Theorem 2 (from [START_REF] Courrieu | Straight monotonic embedding of data sets in Euclidean spaces[END_REF]). Let G be a symmetric positive semi-definite matrix of order n  n .

Then there is an upper triangular matrix R such that R' R  G , r ii  0, 1  i  n , and if for an index i one has r ii  0, then r ij  0, 1  j  n. Moreover, the matrix R with these properties is unique.

Proof. A proof of this is available in ([3], Theorem 4).☐

The corresponding algorithm for computing the generalized Cholesky factor R defined in Theorem 2 is a very simple variant of the usual Cholesky factorization algorithm, and its computational complexity is the same. However, the generalization has the advantage of providing a suitable factor in all cases, even if the matrix G is singular.

Algorithm 1. {Generalized Cholesky factor R of the given matrix G}

r ij  0, 1  i, j  n {initialization of R} r 11  g 11 for j  2 to n for i  1 to j if i  j then r ii  g ii  r ki 2 k1 i1  else if r ii  0 then r ij  (g ij  r ki r kj k1 i1  )/r ii {else r ij  0 as a result of the initialization}.
By construction, the output of Algorithm 1 is an upper triangular factor R with r nonzero rows, and n  r zero rows, where r is the rank of G. The algorithm complexity is in O(n 3 ) , but the exact operations count depends on r and the indices of zero rows. When r  n, this count is maximum, and it is equal to that of the classical Cholesky factorization (plus n(n 1) /2 low cost tests).

Fast Moore-Penrose inverse based solution

Using Algorithm 1, one can define a fast method for computing the Moore-Penrose inverse of every finite matrix. Before examining this method, we rapidly recall some definitions and notations concerning generalized inverses.

Every finite matrix A has (possibly an infinite number of) generalized inverses (hereafter denoted B) that satisfy one or several of the following four Penrose equations:

ABA  A (P1) BAB  B (P2) (AB)' AB (P3) (BA)' BA (P4)
Every matrix B that satisfies the equation set {Pi, Pj,...} is said to be a {i, j,...}-inverse of A, and it is usually denoted A (i, j,...) . The Moore-Penrose inverse of A is the unique matrix [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF] . For a complete explanation, the reader can see [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF].

A †  A (1,2,3 ,
There are several methods for computing the Moore-Penrose inverse, the most usual being based on the singular values decomposition (SVD). This method is numerically very stable, however it is computationally heavy and hardly usable in many practical applications.

Another usual method is based on Gram-Schmidt orthonormalization, which is clearly faster than SVD. However, the classical Gram-Schmidt orthonormalization (CGS or GSO) is known to be numerically instable. A simple and effective remediation to this drawback has been proposed in the form of a re-orthogonalization additional step, leading to the CGS2 method [START_REF] Giraud | Rounding error analysis of the classical Gram-Schmidt orthogonalization process[END_REF]. However, the additional step in CGS2 slows down the process, while it has been observed that CGS is not the fastest method for computing the Moore-Penrose inverse [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF]. In fact, it turned out that among the most usual methods, including Greville's method, SVD, CGS/GSO, and iterative methods, the fastest known numerical method for computing the Moore-Penrose inverse is based on Algorithm 1 and on the following result [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF]:

Theorem 3 (from [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF]). Let A be a m  n matrix, with m  n, set G  A' A , compute the generalized Cholesky factorization G  R' R using Algorithm 1, remove all zero rows from R,

9/21

which results in a full row rank matrix S of size r  n, with r  n, and such that S' S  G .

Then:

A †  S'(SS') 1 (SS') 1 SA'.
Proof. The proof is available in [START_REF] Courrieu | Fast computation of Moore-Penrose inverse matrices[END_REF]. Since it is short, we provide it hereafter.

We start with equation (3.2) from [START_REF] Rakha | On the Moore-Penrose generalized inverse matrix[END_REF], that is:

(EF) †  F'(E' EFF') † E'. ( 7
)
Setting E  A, and F  I in [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF], one obtains:

A †  (A' A) † A' G † A'.
Setting E  S', and F  S in [START_REF] Maronna | Robust Statistics: Theory and Methods[END_REF], one obtains:

G †  (S' S) †  S'(SS' SS') † S  S'(SS') 1 (SS') 1 S , (8)
since SS' is invertible because S is of full row rank. ☐

If A is a m  n matrix, with m  n, it suffices to use the relation A †  ((A') † )'. Note also that (8) provides a simple formula for the Moore-Penrose inverse of any symmetric positive semi-definite matrix, and that if S is of full rank r  n, then G †  G 1 .

Corollary 2. Set A  H 1/ 2 X in Theorem 3, then the minimum norm solution of ( 1) is:

C  S'(SS') 1 (SS') 1 SX'WY ,
where S is defined as in Theorem 3.

Proof. This immediately follows from Theorem 3 and Corollary 1 (ii). ☐

Fast {1, 2, 3}-inverse based solution

Although Corollary 2 provides a fast solution to (1), this is not necessarily the fastest way of solving this problem. Moreover, observing equation [START_REF] Rakha | On the Moore-Penrose generalized inverse matrix[END_REF], one can suspect a potential numerical instability whenever the matrix SS' is ill-conditioned, worsened by the fact that the factor (SS') 1 is repeated. In this section, we describe a {1, 2, 3}-inverse based solution to [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF] whose computational complexity is similar to that of (3), using Algorithm 1 and a simple variant of the backward substitution method for inverting upper triangular matrices. The simplicity of this solution allows us to expect not only faster computation, but also better numerical stability than the Moore-Penrose inverse based approach. We first define the generalization of the backward substitution method for computing generalized inverses of generalized Cholesky factors as they are defined in Theorem 2. The computational complexity of the generalized algorithm is the same as that of the original backward substitution method ( O(n 3 ) ). The algorithm is designed to solve in U the following equation:

RU  I R , (9)
where I R is a diagonal n  n matrix whose ith diagonal coefficient is equal to 0 if the ith row of R is zero, else this diagonal coefficient is equal to 1. The algorithm to solve ( 9) is then:

Algorithm 2. {{1, 2, 3}-inverse U of the given generalized Cholesky factor R}

u ij  0, 1  i, j  n {initialization of U } for j  n downto 1 if r jj  0 then {note: this test is optional} for i  j downto 1 if r ii  0 then if i  j then u ii  1/r ii else u ij  ( r ik u kj k i1 j  )/r ii .
First, one chooses the size parameters m1, n1, m2, n2, the rank parameter r  n1, and the maximum ratio, denoted  r , of non-zero eigenvalues of the Gram matrix X' HX to be built. For practical reasons, one must choose the size parameters such that m2  m1  n1. Then one builds two orthogonal Householder matrices:

M  I  2 uu' u' u
, with a random column vector u  R m1 ,

N  I  2 vv' v'v , with a random column vector v  R n1 ,
where the identity matrices ( I) have the appropriate size in each case. One also builds a r  r

diagonal matrix D, whose ith diagonal coefficient is equal to  r (r  i ) / 2(r  1) , 1  i  r.
Then one selects the first r columns of M and the first r rows of N , and one builds the matrix:

A  M :,1:r DN 1:r,: .

The m1  n1 matrix A is of rank r , the greatest eigenvalue of A' A is equal to  r , and the lowest non-zero eigenvalue of A' A is equal to 1. Thus, we can set X' HX  A' A, that is

H 1/ 2 X  A.
For the next step, one builds a random (m1  r)  n2 real matrix F , and one sets:

P  M :, (r 1):m1 F ,
where we note that the columns of P are orthogonal to those of A.

One can now build a suitable diagonal matrix H  (h ii ),1  i  m1, by taking:

h ii  max | a ij j1 n1  |,| p ij j1 n 2  |       2 ,
which guarantees that both A and P can be factorized with H 1/ 2 as the first factor, in order to build a coherent problem, and one obtains the first matrix of problem (1) by:

X  (H 1/ 2 ) † A .
For the next step, one builds a random n1  n2 real matrix V , and one sets:

HZ  H 1/ 2 (AV  P).
It remains to build a m1  m2 matrix W , with non-negative coefficients, such that

w ij  h ii ,1  i  m1 j1 m 2



, and such that there is a m2  n2 matrix Y that is solution of the equation WY  HZ . Unfortunately, there is no available deterministic method for factorizing HZ in a suitable way, thus we must use a random "trials and errors" approach, as follows.

Repeat the following four steps until WY  HZ (at the working precision):

-build a m1  m2 matrix T  (t ij ) with non-negative random coefficients, -compute the diagonal matrix

K with k ii  t ij j1 m 2  ,1  i  m1, -set W  HK 1 T , -set Y  W † HZ .
The Moore-Penrose inverse W † can be obtained using an accurate (slow) SVD method. In general, one obtains a suitable solution quite rapidly when m2  m1, and  r is not too large.

However, one can observe that the above process frequently fails for large systems if  r  2 12 , which seems to be a practical limit for generating problems in common computational environments such as Matlab.

We have now suitable matrices X , Y , and W for problem [START_REF] Ben-Israel | Generalized Inverses: Theory and Applications[END_REF], and it remains to compute the exact weighted sum of quadratic residues (i.e. the minimized E function of (1))

corresponding to this problem as a reference value for testing the accuracy of solving methods. In order to do this, we use the fact that the columns of the matrix P are orthogonal to those of A, and the proof of Theorem 1. Then one obtains:

E exact  P 2  d ik k1 n 2  i1 m1  ,
where the additional terms ( d ik ) are defined as in (6).

Results

Using the procedure described in Section 4. As one can see in Table 1, the {1, 2, 3}-inverse based method is faster than the fast Moore-Penrose inverse based method, in all cases. One can also note that rank deficient systems are solved faster than full rank systems of the same size by both methods, which is a consequence of the zeroing of certain rows in Algorithm 1. Moreover, an inspection of Table 2 clearly shows that the {1, 2, 3}-inverse based method is accurate in all cases, while the fast Moore-Penrose inverse based method is less accurate and highly sensitive to the parameter  r , thus numerically instable. In summary, it seems that the {1, 2, 3}-inverse based method is globally preferable to other known methods suitable for solving problem (1), except if, for some particular reason, one requires a minimum norm solution. However, in this last case, it is certainly preferable to use an accurate and numerically stable method for computing the Moore-Penrose inverse, but the price to be paid for this is, in general, a quite long computation time.

Conclusion

We have first generalized "weighted least-squares" (WLS) to "weighted pairing leastsquares" (WPLS) problems. This generalization, which allows using a rectangular weight matrix, includes, as particular cases, the classical weighted and non-weighted least-squares problems, and it is more particularly suitable in the framework of data alignment problems.

We have shown that WPLS problems can always be reduced to problems having the same form as WLS problems, and we have studied two fast methods for solving such problems in the case of rank deficient systems as well as of full rank systems. Numerical experiments clearly showed that the best solving method, in terms of speed, accuracy, and numerical

Table 1 .

 1 3}-inverse based solution). In each case, the solving time was recorded in milliseconds (in Matlab 7.5, on a MacBook computer), and the accuracy of each method was measured by (E method  E exact )/ E exact . The mean solving times are reported in Table1, and the mean accuracy values are reported in Table2. All differences between the two methods are statistically significant (p<.01) using the sign test. Mean solving time (in milliseconds) of WPLS problems by the two methods.

	2, we generated various problems of type

Table 2 .

 2 Mean accuracy of the two methods in solving WPLS problems.

	n1 :  r :	16	128 256	4096	16	256 256	4096	16	512 256	4096
	Full rank
	WPLSdagger	8.9E-6
	WPLS123	< 10 -12
	Rank deficient
	WPLSdagger
	WPLS123

The test at the third line of Algorithm 2 is optional because it has no influence on the result. However, including this test allows to save a number of useless floating-point operations (whose result is zero) whenever R is singular.

By construction, the output of Algorithm 2 is an upper triangular matrix U that has the property that if the ith row of R is zero, then both the ith row and the ith column of U are zero. Note that the rank of U is equal to the rank of R, which is itself equal to the rank of the factorized matrix G and to the trace of I R . Note also that if R is not invertible (in the usual sense), then UR  I R , however, UR is always idempotent since URUR  UI R R  UR . Theorem 4. Let R be a generalized Cholesky factor as defined in Theorem 2, let U be the

R and has the same rank as R, then the conclusion follows from ([1], p. 46). Alternatively, one can easily verify that

(1, 2, 3) using Algorithm 2.

Then:

(ii) The matrix B  UU' A' is a {1, 2, 3}-inverse of A.

Proof.

Since (i) implies that B  UQ' , one has:

which proves (ii), and then completes the proof of Theorem 5. ☐ Corollary 3. Set A  H 1/ 2 X in Theorem 5, then a solution to ( 1) is given by:

where U is defined as in Theorem 5. Proof. This immediately follows from Theorem 5 (ii) and Corollary 1 (iii). ☐

One can note that if H 1/ 2 X is of full column rank, then U  R 1 , and the solution provided by Corollary 3 is equal to that of Corollary 1 In order to make the performance of the two tested methods comparable, we avoided the use of high level Matlab operators such as "inv()", "chol()", or "\", whose implementation is hidden and compiled.

Test problems

In order to test the performance of methods for solving (1) in terms of speed, accuracy, and numerical stability, we must build test problems in a way that allows strict control of relevant characteristics such as the size and the rank of the equation system, the exact weighted least-squares residue norm, and the ratio of extreme non-zero singular values of the system matrix (which can be seen as a kind of generalized condition number). Note that the solution ( C ) itself is not relevant for comparisons, since it is not unique in the case of rank deficient systems. Building coherent test problems having all required properties is not so easy, and we propose the following method.

stability, is based on a special {1, 2, 3}-inverse whose computation is very simple. In contrast, approaches based on the Moore-Penrose inverse lead to slow computation, or alternatively to numerical instability.

Appendix 1

The following codes are provided for example, and for academic use only. The code is not optimized and exception cases are not managed.

Matlab code (version 7.5) of the WPLSdagger function:

function [C,Emethod,time] = WPLSdagger(X,Y,W) % Moore-Penrose inverse based solution of a WPLS problem tic % start the clock [m1,n1]=size(X); [m2,n2]=size(Y); H=sum(W,2); g=X'*((H*ones(1,n1)).*X); % s = full row rank generalized Cholesky factor of g tol=n1*eps(norm(g,inf)); s=zeros(n1,n1); ii=0; for i=1:n1


Gaël Poux-Médard
email: gael.poux-medard@univ-lyon2.fr

Julien Velcin
email: julien.velcin@univ-lyon2.fr

Sabine Loudcher
email: sabine.loudcher@univ-lyon2.fr

Multivariate Powered Dirichlet-Hawkes Process

Keywords: Dirichlet Process, Multivariate Hawkes Process, Clustering, Information spread, Sequential data 01 Jan 31 Jan 15 Jan trump, china, hong, kong, climate police, woman, kill, arrest, charge china, arrest, kill, woman, abuse police, woman, yellow, arrest, protest government, shutdown, worker, court, family trump, shutdown, house, border, government brazil, president, bolsonaro, right, india climate, change, people, brazil, world trump, woman, text, syria, saudi saudi, woman, china, flee, family trump, russia, mueller, report,

The publication time of a document carries a relevant information about its semantic content. The Dirichlet-Hawkes process has been proposed to jointly model textual information and publication dynamics. This approach has been used with success in several recent works, and extended to tackle specific challenging problems -typically for short texts or entangled publication dynamics. However, the prior in its current form does not allow for complex publication dynamics. In particular, inferred topics are independent from each other -a publication about finance is assumed to have no influence on publications about politics, for instance. In this work, we develop the Multivariate Powered Dirichlet-Hawkes Process (MPDHP), that alleviates this assumption. Publications about various topics can now influence each other. We detail and overcome the technical challenges that arise from considering interacting topics. We conduct a systematic evaluation of MPDHP on a range of synthetic datasets to define its application domain and limitations. Finally, we develop a use case of the MPDHP on Reddit data. At the end of this article, the interested reader will know how and when to use MPDHP, and when not to.

Introduction

Understanding the data publication mechanisms on online platforms is of utmost importance in computer science. The amount of user-generated content that flows on social networks (e.g. Reddit) daily appeals for efficient and scalable approaches; they should provide us detailed insights within these mechanisms. A favoured approach to this problem is to cluster published documents according to their semantic content [START_REF] Blei | Dynamic topic models[END_REF][START_REF] Blei | Latent dirichlet allocation[END_REF][START_REF] Yin | Model-based clustering of short text streams[END_REF].

In the specific case of data flowing on social networks, time also carries a valuable information about the underlying data flow generation process [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF][START_REF] Tan | The indian buffet hawkes process to model evolving latent influences[END_REF]. Typically, we expect given publications to trigger ulterior publications within a Probability clusters Intensity red cluster on red cluster Intensity red cluster on blue cluster Intensity blue cluster on red cluster Intensity blue cluster on blue cluster Fig. 1. Illustration of the Multivariate Powered Dirichlet-Hawkes Process prior -A new event appears at time t = 3 from a cluster which is yet to be determined. The a priori probability that this event belongs to a given cluster c red depends on the sum of the red intensity functions at time t = 3. Similarly, the a priori probability that this event belongs to a cluster c blue depends on the sum of the blue lines at time t = 3. In previous models, this prior probability depends on each cluster self-stimulation only. short time range. This effect has been studied by considering spreading agents, who are individually influenced by contacts' publications [START_REF] Gomez-Rodriguez | Uncovering the temporal dynamics of diffusion networks[END_REF][START_REF] Gomez-Rodriguez | Modeling information propagation with survival theory[END_REF][START_REF] Poux-Médard | Dirichlet-survival process: Scalable inference of topic-dependent diffusion networks[END_REF].

In previous works, the understanding of large data flows boils down to sorting data pieces (documents) into independent topics (clusters). However, it has been underlined on several occasions that online publication mechanisms are more complex than that. Typically, a correct description should involve clusters that interact with each other [START_REF] Myers | Clash of the contagions: Cooperation and competition in information diffusion[END_REF][START_REF] Poux-Médard | Interactions in information spread: quantification and interpretation using stochastic block models[END_REF][START_REF] Zarezade | Correlated cascades: Compete or cooperate[END_REF]. We illustrate the implications of this claim in Fig. 1. In most existing works that explicitly model both text and time, a given topic is assumed to only trigger observations from the same topic [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF] -the red cluster can only trigger observations from the red cluster. Instead, we must allow clusters to trigger publications in any other cluster.

Therefore, we extend a previous class of models (DHP [4] and PDHP [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]) to account for cluster interaction mechanisms. We show that technical challenges arise when considering topical interaction, and solve them. This results in the Multivariate Powered Dirichlet-Hawkes Process (MPDHP). We conduct systematic experiments to test the limits of MPDHP and define its application domain.

In particular, we show that it performs well in cases when textual data is scarce and when the number of coexisting clusters is large. Finally, we investigate a real-world use case on a Reddit dataset.

Background

The original Dirichlet-Hawkes process (DHP) [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF] merges Dirichlet processes and Hawkes processes. It is used as a prior in Bayesian clustering along with a main model -typically a language model. The prior expresses the assumption that a new event from a given topic appears conditionally on the presence of older events from this same topic. The conditional probability is encoded into the intensity function of a Hawkes process. One such Hawkes process is associated to each topic. The temporal (Hawkes) intensity of a topic c is written λ c (t|H c ); it depends on the history of all previous events associated to topic c, written H c . If no Hawkes intensity manages to explain well enough the presence of a new observation happening further in time, the DHP a priori guess is that the new observation belongs to a new cluster. DHP have first been used for automated summary generation [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF]. A list of textual documents appear in chronological order and are treated as such; the DHP infers clusters of documents that are based both on their textual content and their publication date, and studies their auto-stimulated publication dynamics. This process knew several developments, that essentially consider alternative Dirichlet-based priors combined with Hawkes processes -Hierarchical DHP [START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF], Indian Buffet Hawkes process [START_REF] Tan | The indian buffet hawkes process to model evolving latent influences[END_REF] and powered DHP [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF].

However, in [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF][START_REF] Yin | Model-based clustering of short text streams[END_REF], the authors underline several limits of the standard Dirichlet-Hawkes processes and of the extensions mentioned earlier. For instance, DHP fails in cases where publications content carry few information: when textual content is short (e.g. tweets [START_REF] Yin | Model-based clustering of short text streams[END_REF]) or when vocabularies overlap significantly (e.g., topic-specific datasets). In [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF], the problem is alleviated by considering a Powered Dirichlet process [START_REF] Poux-Médard | Powered dirichlet process for controlling the importance of "rich-get-richer" prior assumptions in bayesian clustering[END_REF] instead of a standard Dirichlet process. This process is merged with a univariate Hawkes process to make the Powered Dirichlet-Hawkes process. The authors retrieve better results in challenging clustering situations (large temporal and textual overlaps).

However, none of these works allow clusters to interact with each other, despite clues pointing in that direction [START_REF] Myers | Clash of the contagions: Cooperation and competition in information diffusion[END_REF][START_REF] Poux-Médard | Interactions in information spread: quantification and interpretation using stochastic block models[END_REF][START_REF] Zarezade | Correlated cascades: Compete or cooperate[END_REF]. Indeed in [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF][START_REF] Tan | The indian buffet hawkes process to model evolving latent influences[END_REF], the considered Hawkes processes are univariate: a cluster can only be used to trigger events within itself. Exploring how clusters interact with each other would significantly extend our comprehension of the publication mechanisms at stake in various datasets -such as social media or scientific articles. Identifying which topics trigger the publication of other seemingly unrelated topics might be interesting in the study of fake news spreading. Understanding the dynamics at stake may help to surgically inhibit the spread of such topics using the right refutation. Another possible use case would be nudging users towards responsible behaviours regarding environment, health, tobacco, etc.

In this paper, we extend the (univariate) Powered Dirichlet-Hawkes Process to its multivariate version. There are several reasons why it has not been done in prior works: first of all, the adaptation to the multivariate case is not trivial and poses several technical challenges. As a first contribution, we detail the challenges that arise when developing the Multivariate Powered Dirichlet-Hawkes Process (MPDHP). We propose methods to overcome them while retaining a linear time complexity O(N ). Doing so, we also relax the near-critical Hawkes process hypothesis made in [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]. The second reason why the multivariate extension has not been developed in prior works, is that it greatly raises the number of parameters to estimate. The inference task might become harder, and the results irrelevant. Our second contribution consists in a systematic evaluation of the MPDHP on a variety of synthetic situations. Our goal is to identify the limits of MPDHP regarding textual overlap, computation time, the amount of available data, the number of co-existing clusters, etc. We show that MPDHP is perfectly fit for solving a variety of challenging situations. Finally, we illustrate the new insights on topical interaction obtained by running MPDHP on a real-world Reddit dataset.

3 The Multivariate Powered Dirichlet-Hawkes process

Powered Dirichlet Process

The Dirichlet process can be expressed using the Chinese Restaurant Process metaphor. Consider a restaurant with an infinity of empty tables. A first client enters the restaurant and sits to any of the empty tables with a probability proportional to α -the concentration parameter. Another client then enters the restaurant, and sits either at one of the occupied tables with a probability linearly proportional to the number of clients already sat at the table, or to any of the empty tables with a probability proportional to α. The process is then iterated for an infinite number of clients. The resulting clients distribution over the tables is equivalent to a draw from a Dirichlet distribution. The Powered Dirichlet process is intended as a generalisation of the Dirichlet process [START_REF] Poux-Médard | Powered dirichlet process for controlling the importance of "rich-get-richer" prior assumptions in bayesian clustering[END_REF].

The probability for a new client to sit at one of the occupied tables is now proportional to the number of clients at the power r. Let C i be the table chosen by the i th client and H the history of table allocation. Formally, the probability for the i th to choose a table reads:

P DP (C i = c|α, r, H) = N r c α+ c N r c if c = 1, ..., K α α+ c N r c if c = K+1 (1)
where N c is the number of people that already sat at table c, K is the total number of tables, and r a hyper-parameter. Note that when r = 1 we recover the regular Dirichlet process, and when r = 0 we recover the Uniform Process [START_REF] Wallach | An alternative prior process for nonparametric bayesian clustering[END_REF].

Multivariate Hawkes process

A Hawkes process is a temporal point process where the appearance of new events is conditional on the realisation of previous events. It is fully characterised by an intensity function, noted λ(t|H) that depends on the history of previous events. It is interpreted as the instantaneous probability of a new observation appearing:

λ(t)dt = P (e ∈ [t, t + dt]
) with e an event and dt an infinitesimal time interval. For simplicity of notation, we omit the term H which is implicit anytime the intensity function λ is mentioned. As in the DHP [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF], we define one Hawkes process for each cluster. However in DHP each of them is associated to a univariate Hawkes process, that depends only on the history of events comprised in this cluster. In our case, instead, we associate each cluster to a multivariate Hawkes process that depends on all the observations previous to the time being. Let t c i be the time of realisation of the i th event belonging to cluster c. We write the intensity function for cluster c at all times as:

λ c (t) = t c ′ i <t α c,c ′ • κ(t -t c ′ i ) ; κ l (∆t) = e - (∆t-µ l ) 2 2σ 2 l 2πσ 2 l ( 2
)
In Eq. 2, α c,c ′ is a vector of L parameters to infer, and κ(t -t c ′ i ) is a vector of L temporal kernel functions depending only on the time difference between two events. In our case, we consider a Gaussian RBF kernel, that allows to model a range of different intensity functions.

The log-likelihood of a multivariate Hawkes process for all observations up to a time T reads:

log L(α|H) = c T 0 λ c (t)dt + t c i λ c (t c i ) (3)

Multivariate Powered Dirichlet-Hawkes Process

The Multivariate Powered Dirichlet-Hawkes Process (MPDHP) arises from the merging of the Powered Dirichlet Process (PDP) and the Multivariate Hawkes Process (MHP), described in the previous sections. As in [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF], the counts in PDP are substituted with the intensity functions of a temporal point-process, here MHP. The a priori probability that a new event is associated to a given cluster no longer depends on the population of this cluster, but on its temporal intensity at the time the new observation appears. This is illustrated in Fig. 1, where two events from two different clusters c red and c blue have already happened at times t 0 = 0 and t 1 = 1. A new event appears at time t = 3. The a priori probability that this event belongs to the cluster c red depends on the sum of the intensity functions of observations at t 0 and t 1 on cluster c red at time t = 3 -sum of the red dotted lines. Similarly, the a priori probability that this event belongs to the cluster c blue depends on the sum of the blue dotted lines at time t = 3. Let t i be the time at which the i th event appears. The resulting expression reads:

P (C i = c|t i , r, λ 0 , H) =    λ r c (ti) λ0+ c ′ λ r c ′ (ti) if c≤K λ0 λ0+ c ′ λ r c ′ (ti) if c=K+1 (4)
In Eq. 4, λ c in defined as in Eq. 2, and the parameter λ 0 is the equivalent of the concentration parameter described in Eq. 1. Taking back the illustration in Fig. 1, this parameter corresponds to a time-independent intensity function. It has a chance to get chosen typically when the other intensity functions are below it (meaning they do not manage to explain the dynamic aspect of a new event). In this case, a new topic is opened, and gets associated to its own intensity function.

Language model

Similarly to what has been done in [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF], the MPDHP must be associated to a Bayesian model given it is a prior on sequential data. Since we study applications on textual data, we choose to side the MPDHP prior with the same Dirichlet-Multinomial language model as in previous publications [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]. According to this model, the likelihood of the i th document belonging to cluster c reads:

L(C i = c|N <i,c , n i , θ 0 ) = Γ (N c + θ 0 ) Γ (N c + n i + θ 0 ) v Γ (N c,v + n i,v + θ 0,v ) Γ (N c,v + θ 0 ) (5)
where N c is the total number of words in cluster c from observations previous to i, n i is the total number of words in document i, N c,v the count of word v in cluster c, n i,v the count of word v in document i and θ 0 = v θ 0,v . Note that for any empty cluster, the likelihood is computed using N cempty = 0 for every empty cluster c empty .

Implementation

Base algorithm SMC algorithm

We use a Sequential Monte Carlo (SMC) algorithm for the optimisation. The base algorithm is the same as in [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]] -a graphical representation of the SMC algorithm is provided in [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]-Fig. 1 and as Supplementary Material. The goal of the SMC algorithm is to jointly infer textual documents' clusters and the dynamics associated with them. It runs as follows. First, the algorithm computes each cluster's posterior probability for a new observation by multiplying the temporal prior on cluster allocation (see Eq. 4, illustrated Fig. 1) with the textual likelihood (see Eq. 5). It results in an array of K + 1 probabilities, where K is the number of non-empty clusters. A cluster label is then sampled from this probability vector. If the empty (K + 1) th cluster is chosen, the new observation is added to this cluster, and its dynamics are randomly initialised (i.e. a (K + 1) th row and a (K + 1) th column are added to the parameters matrix α). If a non-empty cluster is chosen, its dynamics are updated by maximising the new likelihood Eq. 3. The process then goes on to the next observation. This routine is repeated N part times in parallel. Each parallel run is referred to as a particle. Each particle keeps track of a series of cluster allocation hypotheses. After an observation has been treated, we compute the particles likelihood given their respective cluster allocations hypotheses. Particles that have a likelihood relative to the other particles' one below a given threshold ω thres are discarded and replaced by a more plausible existing particle.

Sampling the temporal parameters The parameters α are inferred using a sampling procedure. A number N sample of precomputed vectors is drawn from a Dirichlet distribution with probability P (α|α 0 ), with α 0 a concentration parameter. As the SMC algorithm runs, within each existing cluster, each of these candidate vectors is associated to a likelihood computed from Eq. 3, noted P (H|α), where H represents the data. The sampling procedure returns the average of each of the N sample precomputed α, weighted by the posterior distribution associated to them P (α|H) ∝ P (H|α)P (α|α 0 ). The so-returned matrix is guaranteed to be a good statistical approximation of the optimal matrix, provided the number of sample matrices N sample is large enough.

Limits This algorithm described here works well for the univariate case, but fails for the multivariate case. In particular, updating the multivariate intensity function of each cluster requires knowing the number of already existing clusters, which vary over time. Therefore, we cannot precompute the sample matrices in advance -they must be updated as the algorithm runs to account for the right number of non-empty clusters. Moreover, the number of parameters to estimate evolves linearly with the number of active clusters K, instead of remaining constant as in DHP and variants [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]. Because the number of parameters is not constant anymore, their candidate values cannot be sampled from a Dirichlet distribution anymore. In the following, we review these challenges and present our solutions to overcome them. We manage to preserve a constant time complexity for each observation.

Optimisation challenges

Temporal horizon A first problem that has been answered in [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF] is that, for each new observation, the algorithm has to run through the whole history of events to compute the DHP prior. However, carefully choosing the kernel vector κ(•) described in Eq. 2 allows to perform this step in constant time. If the chosen kernel vanishes as time goes, it happens a point where old events have a nearzero chance to have any influence on new observations, according to our model. These events can be discarded from the computation for new events.

Updating the triggering kernels In the univariate case [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Mavroforakis | Modeling the dynamics of learning activity on the web[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF], the coefficients α c ∈ R L are sampled from a collection of existing sample vectors computed at the beginning of the algorithm (where L is the size of the kernel vector). However, we must now infer a matrix instead. We recall that matrix α c represents the weights given to the temporal kernel vector of every cluster influence on c -see Eq. 2. The likelihood Eq. 3 can be updated incrementally for each sample matrix. A given cluster c has a likelihood value associated to each of those N sample sample matrices, which represents how fit one sample matrix is to explain one cluster's dynamics. The final value of the parameters matrix is sampled simply by averaging the samples matrices weighted by their likelihood for a given cluster times the prior probability of these vectors being drawn in the first place.

Such sampling was possible in the univariate case, where each sample matrix was in fact a vector of fixed length. In our case, because Hawkes processes are multivariate, each entry α c ∈ R K×L is now a matrix. Moreover, the number of existing clusters K increases over with time, and can grow very large. Each time a new cluster is added to the computation, a row is appended to the α matrix -it accounts for the influence of this new cluster regarding c.

However, some older events can be discarded from the computation. When an event is older than 3σ with respect to the longest range entry of the RBF kernel, it can be safely discarded. Clusters whose last observation has been discarded thus have a near-zero chance to get sampled once again. These clusters' contribution to the likelihood Eq. 3 will not change anymore. Therefore, they do not have a role in the computation of the parameters matrix α c . The row corresponding to each of these clusters in the parameters matrix can then be discarded in every sample matrix. Put differently, the last sampled value for their influence on c will remain unchanged for the remaining of the algorithm. The dimension of α c only depends on the number of active clusters, whose intensity function has not faded to zero. For a given dataset, the number of active clusters typically fluctuates around a constant value, making one iteration running in constant time O(1).

A beta prior on parameters Another problem inherent to the multivariate modelling is the prior assumption on sample vectors. In [START_REF] Du | Dirichlet-hawkes processes with applications to clustering continuous-time document streams[END_REF][START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF], each sample vector is sampled from a Dirichlet distribution. This choice is to infer Hawkes processes that are nearly-unstable: the spectral radius of the temporal kernel function λ c (t) is close to 1. However in our case, such assumption is not possible because the size of each sample matrix can vary as the number of active clusters evolve. Drawing one Dirichlet vector for each entry α c,c ′ would force the spectral radius of α c to equal K = |α c |, which transcribes a highly-unstable Hawkes process. Our solution is to consider the parameters as completely independent from each other. Each entry of the matrix α is drawn from an independent β distribution of parameter β 0 . In this way, we make no assumption on the spectral radius of the Hawkes process, and samples rows/columns corresponding to new clusters can be generated one after the other.

Experiments

Setup

We design a series of experiments to determine the use cases of the Multivariate Powered Dirichlet-Hawkes Process 1 . We list the parameters we consider in our experiments. When a parameter does not explicitly vary, it takes a default value given between parenthesis. These parameters are: the textual overlap (0) and the temporal overlap (0) discussed further in the text, the temporal concentration parameter λ 0 (0.01), the strength of temporal dependence r (1), the number of synthetically generated clusters K (2), the number of words associated to each document n words (20), the number of particles N part [START_REF] Posadas Duran | Detection of fake news in a new corpus for the spanish language[END_REF] and the number of sample matrices used for sampling α, noted N sample (2 000). For the detail of these parameters, please refer to Eq. 4. The interplay between the parameters that are not part of MPDHP (N sample and N part ) is studied in Supplementary Material. Note that the overlap o(f 1 , ..., f N ) between N functions is defined as the sum over each function f i of its intersecting area with the largest of the N -1 other functions, divided by the sum of each function's total area [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]. This value is bounded between 0 (perfectly separated functions) and 1 (identical functions):

o(f 1 , ..., f N ) = i R min(f i (x), max({f j (x)} j̸ =i ))dx
For each combination of parameters considered, we generate 10 different datasets. In all datasets, we consider a fixed size vocabulary V = 1000 for each cluster. All datasets are made of 5 000 observations. Observations for each cluster c are generated using a RBF temporal kernel κ(t) weighted by a parameter matrix α c . We set κ(t) = [G(3; 0.5); G(7; 0.5); G(11; 0.5)] where G(µ; σ) is a Gaussian of mean µ and standard deviation σ. We note L = 3 the number of entries of κ.

The inferred entries of α determine the amplitude (weight) of each kernel entry. The generation process is as follows. First, we draw a random matrix α ∈ R K×L and normalise it so that its spectral radius equals 1 -near unstable Hawkes process. We repeat this process until we obtain the wanted temporal overlap. Then, we simulate the multivariate Hawkes process using the triggering kernels α • κ(t), where κ(t) is the RBF kernel as defined earlier. Given the Hawkes process is multivariate, each event is associated to its class it has been generated from among K possible classes. For each event, we draw n words words from a vocabulary of size V . The vocabularies are drawn from a multinomial distribution and shifted over this distribution so that they overlap to a given extent.

Baselines

We evaluate our clustering results in terms of Normalised Mutual Information score (NMI). This metric is standard when evaluating non-parametric clustering models. It compares two cluster partitions (i.e. the inferred and the ground truth ones); it is bounded between 0 (each true cluster is represented to the same extent in each of the inferred ones) and 1 (each inferred partition comprises 100% of one true cluster). The standard error is computed on 100 runs. We compare our approach to 3 closely related baselines. Powered Dirichlet-Hawkes process (PDHP) [START_REF] Poux-Médard | Powered hawkes-dirichlet process: Challenging textual clustering using a flexible temporal prior[END_REF]: in this model, clusters can only self replicate. It means that the intensity function of a cluster c Eq. 2 only considers past events that happened in the same cluster c. r is set to 1. Dirichlet process (DP): this prior is standard in clustering problems. It corresponds to a special case of Eq. 1 where r = 1. The prior probability for an observation to belong to a cluster depends on its population. Uniform process (UP) [START_REF] Wallach | An alternative prior process for nonparametric bayesian clustering[END_REF]: this prior corresponds to a special case of Eq. 1 where r = 0. It assumes that the prior probability for an observation to belong to a cluster does not depend on any information about this cluster (neither population nor dynamics).

Results on synthetic data

MPDHP outperforms state-of-the-art In Fig. 2, we plot our results for datasets that have been generated using a Multivariate Hawkes process (clusters have an influence on each other) and a Univariate Hawkes process (clusters can only influence themselves). We compare MPDHP to our baselines for various values of textual overlap -we provide a similar study that considers various temporal overlaps as Supplementary Material. MPDHP systematically outperforms the baselines on multivariate data -when clusters interact with each other. Considering that clusters interacts with each other improves our description of the datasets. MPDHP performs at least as good as PDHP on univariate datawhen clusters can only self-stimulate. The complexity of MPDHP does not make it unfit to simpler tasks. PDHP performs better than MPDHP when the textual overlap is large (textual overlap of 0.8) due to its reduced complexity. Increasing the number of observations would fix this.

Highly interacting processes

We test when a large number of clusters coexist simultaneously. The rate at which new clusters get opened is mainly controlled by the λ 0 hyperparameter (see Eq. 4), which we vary to see whether MPDHP is robust against it. In Fig. 3 (left), we plot the performances of MPDHP according to these two parameters. We draw two conclusions: MPDHP can handle a large number of coexisting clusters and still correctly identify to which one each document belongs, and MPDHP is robust against large variations of λ 0 . In this case, results are similar for λ 0 varying over 5 orders of magnitude. It means MPDHP does not have to be fine-tuned according to the number of expected clusters in cases where this number is not known in advance. An extended discussion on the choice of the parameter λ 0 is provided as Supplementary Material.

Handling scarce textual information In this paragraph, we determine how much data should be provided to MPDHP to get satisfying results. In Fig. 3 (right), we plot the performances of MPDHP with respect to the number of words generated by each observation and to the clusters' vocabulary overlap. MPDHP needs a fairly low number of words to yield good results over 5 000 observations. For reference, the overlap between topics can be estimated around 0.25 ( [START_REF] Posadas Duran | Detection of fake news in a new corpus for the spanish language[END_REF], in Spanish). Similarly, we can estimate an average of ∼10-20 named entities per Twitter post (240 characters). These results support the application of MPDHP to model real-world situations.

Real-world application on Reddit

Data We conclude this work with an illustration of MPDHP in a real-world situation. We investigate the interplay between topics on news subreddits, that is how much influence a topic can exert on other ones. The dataset is collected from the Pushshift Reddit repository [START_REF] Baumgartner | The pushshift reddit dataset[END_REF]. We limit our study to headlines from popular English news subreddits in January 2019: inthenews, neutralnews, news, nottheonion, offbeat, open news, qualitynews, truenews, worldnews. From these, we remove posts that have a popularity (difference between upvotes and downvotes) lesser than 20, as they are of lesser influence in the dataset and only add to the modeling. We remove headlines that contain less than 3 words as they only add noise to the modelling. After curating the dataset in the way described above, we are left with roughly 8,000 news headlines, which makes a total of 65,743 tokens drawn from a vocabulary of size 7,672. Additional characteristics of the dataset are provided as Supplementary Material.

Results

We present the results of MPDHP on real-world data in Fig. 4. Fig. 4a. illustrates a typical output from the model. The transparency in the representation of λ(t) accounts for the number of times such interaction has effectively been observed; transparency of the intensity function λ c,c ′ (t) of c ′ on c is proportional to t c t c ′ <t c λ c,c ′ (t). We can make two interesting observations from this figure. Firstly, the interaction strength between clusters seems to fade as time passes (Fig. 4b.). Cluster interactions are more likely to happens within short time ranges. Secondly, the first two clusters seem to be consistently used across the whole month (Fig. 4c.). When we look at their composition, we notice that the first cluster is made of 75% of articles from the subreddit r/worldnews, which is +20% from what one would expect from chance. Similarly, the second cluster comprises 46% of r/news articles, which is also roughly +20% from expected at random. These two clusters therefore significantly account for publications from either of these subreddits, independently from the textual content. Both are general news forums with a large audience; an article that gets posted on other subreddits is likely to also appear on these. Other clusters follow a bursty dynamic, which concurs with [START_REF] Haralabopoulos | Lifespan and propagation of information in on-line social networks: A case study based on reddit[END_REF]. More details on this experiment can be found in [START_REF] Poux-Médard | Properties of reddit news topical interactions[END_REF].

Conclusion

In this paper, we extended the Powered Dirichlet-Hawkes process so that it can consider multivariate processes, resulting in the Multivariate Powered Dirichlet-Hawkes process (MPDHP). This new process can infer temporal clusters interaction networks from textual data flow. We overcome several optimisation challenges to preserve a time complexity that scales linearly with the dataset. We showed that MPDHP outperforms existing baselines when clusters interact with each other, and performs at least as well as the PDHP baseline when clusters do not (which PDHP is designed to model). MPDHP can handle cases where textual content is lesser informative better than other baselines. It is robust against tuning of the temporal concentration parameter λ 0 , which allows to handle highly intricate processes. We finally showed that MPDHP performs well with scarce textual data. Our results suggest that MPDHP can be applied in a robust way to a broad range of problems, which we illustrate on a real-world application, that provides insights in topical interactions mechanisms between news published on Reddit.

Fig. 2 .

 2 Fig.2. Numerical results on synthetic data -MPDHP consistently outperforms other baselines designed for the univariate case on both univariate and multivariate data. The standard error has been computed using 100 independent runs.

Fig. 3 .

 3 Fig.3. MPDHP can handle a large number of coexisting clusters and scarce textual information -MPDHP yields good results when a large number of clusters coexist simultaneously (left) and when texts are short or little informative (right). It is robust against variations of λ0 over 5 order of magnitude.

ParametersFig. 4 .

 4 Fig. 4. Real-world application on Reddita. Examples of clusters along with their inferred reproduction dynamics. b. Visualisation of interaction patterns at different times as a network; each dot is a cluster, each edge accounts for the value of λ(t) at a given time t ∈ [0; 2; 4; 6; 8] hours. c. Most used clusters represented over real time.

Data and implementations are available at https://github.com/GaelPouxMedard/ MPDHP