Fredrik Johansson
email: fredrik.johansson@gmail.com

ON A FAST AND NEARLY DIVISION-FREE ALGORITHM FOR THE CHARACTERISTIC POLYNOMIAL

We review the Preparata-Sarwate algorithm, a simple O(n 3.5 ) method for computing the characteristic polynomial, determinant and adjugate of an n×n matrix using only ring operations together with exact divisions by small integers. The algorithm is a baby-step giant-step version of the more well-known Faddeev-Leverrier algorithm. We make a few comments about the algorithm and evaluate its performance empirically.

Introduction

Let R be a commutative ring. Denote by ω > 2 an exponent of matrix multiplication, meaning that we can multiply two n × n matrices using O(n ω ) ring operations (additions, subtractions and multiplications). Given a matrix A ∈ R n×n , how fast can we compute its characteristic polynomial, determinant and adjugate (or where applicable, inverse), without dividing by elements in R?

The obvious division-free algorithm is cofactor expansion, which uses O(n!) operations. It is mainly interesting for n ≤ 4 and for sparse symbolic matrices. The first published efficient method is the Faddeev-Leverrier algorithm [START_REF] Verrier | Sur les variations séculaire des élements des orbites pour les sept planétes principales[END_REF][START_REF] Helmberg | On Faddeev-Leverrier's method for the computation of the characteristic polynomial of a matrix and of eigenvectors[END_REF][START_REF] Csanky | Fast parallel matrix inversion algorithms[END_REF][START_REF] Hou | Classroom note: A simple proof of the Leverrier-Faddeev characteristic polynomial algorithm[END_REF]; Alg. 2.2.7 in [START_REF] Cohen | A Course in Computational Algebraic Number Theory[END_REF]), which solves the problem using roughly n matrix multiplications, for a complexity of O(n ω+1 ). The Faddeev-Leverrier algorithm requires some divisions, but it is nearly division-free in the sense that all divisors are fixed small integers 1, 2, . . . , n rather than general elements of R, and the divisions are exact in the sense that the quotients remain in R.

The Berkowitz algorithm [START_REF] Stuart | On computing the determinant in small parallel time using a small number of processors[END_REF] achieves the same complexity using O(n 4 ) operations, or O(n ω+1 log n) with fast matrix multiplication, without performing any divisions. In practice, the Berkowitz algorithm is faster than the Faddeev-Leverrier algorithm by a constant factor and is now widely used in computer algebra systems for linear algebra over rings where divisions are problematic.

The complexity of the Faddeev-Leverrier algorithm can be improved using a baby-step giant-step technique, leading to a method with O(n ω+0.5 +n 3 ) complexity which is asymptotically better than the Berkowitz algorithm. This method was apparently first discovered by S. Winograd who did not publish the result, and then independently discovered and published by Preparata and Sarwate [START_REF] Preparata | An improved parallel processor bound in fast matrix inversion[END_REF].

Unfortunately, many references to Faddeev-Leverrier algorithm in the literature present it as an O(n 4 ) or O(n ω+1 ) algorithm without mentioning the improvement of Preparata and Sarwate. Berkowitz [START_REF] Stuart | On computing the determinant in small parallel time using a small number of processors[END_REF] claims that a baby-step giant-step Faddeev-Leverrier algorithm is possible but only cites private communication with S. Winograd without giving an explicit algorithm. The present author is not aware of any software using the Preparata-Sarwate algorithm.

There are perhaps three reasons for the relative obscurity of the method. 1 First of all, the Faddeev-Leverrier algorithm is numerically unstable [START_REF] Rehman | La Budde's method for computing characteristic polynomials[END_REF][START_REF] Hardy | The algebraic eigenvalue problem[END_REF] making it virtually useless for ordinary numerical computation, and it is in any case much slower than standard O(n 3 ) methods such as Gaussian elimination and Hessenberg reduction which are numerically stable. Second, faster methods are known for the most commonly used exact fields and rings such as F q , Q and Z. Third, Kaltofen and Villard [START_REF] Kaltofen | On computing determinants of matrices without divisions[END_REF][START_REF] Kaltofen | On the complexity of computing determinants[END_REF][START_REF] Villard | Kaltofen's division-free determinant algorithm differentiated for matrix adjoint computation[END_REF] have achieved an even lower division-free complexity of O(n 3.2 ) for the determinant or adjugate and O(n 3.29 ) for the characteristic polynomial with classical multiplication (ω = 3), and about O(n 2.72 ) and O(n 2.84 ) for the respective tasks using the fastest currently known matrix multiplication algorithm due to Le Gall [START_REF] François | Powers of tensors and fast matrix multiplication[END_REF]. The Kaltofen-Villard algorithm is far more complicated than the Preparata-Sarwate algorithm, however.

The contribution of this note is twofold. First, we give explicit pseudocode for a version of the Preparata-Sarwate algorithm (Algorithm 2) that may serve as a convenient reference for future implementers. The code makes a superficial change to the algorithm as it is usually presented [START_REF] Preparata | An improved parallel processor bound in fast matrix inversion[END_REF][START_REF] Codenotti | Parallelism and fast solution of linear systems[END_REF][START_REF] Abdeljaoued | Méthodes matricielles-Introduction à la complexité algébrique[END_REF], halving memory usage (this is not an entirely negligible change since the Preparata-Sarwate algorithm uses O(n 2.5 ) memory). Second, we discuss the applicability of the algorithm and perform computational experiments comparing several algorithms over various rings.

The Faddeev-Leverrier algorithm and its refinement

We first recall the original Faddeev-Leverrier algorithm (Algorithm 1) for the characteristic polynomial

p A (x) = c n x n + • • • + c 1 x + c 0 of a matrix A ∈ R n×n .

Algorithm 1 Faddeev-Leverrier algorithm

Input: A ∈ R n×n where n ≥ 1 and R is a commutative ring, R having a unit element and characteristic 0 or characteristic coprime to 1, 2, . . . , n Output: (p A (x), det(A), adj(A))

1: c n = 1, B ← I, k ← 1 2: while k ≤ n -1 do 3: B ← AB 4: c n-k ← -1 k Tr(B) 5: B ← B + c n-k I 6: k ← k + 1 7: end while 8: c 0 ← -1 n Tr(AB) 9: return (c 0 + c 1 x + . . . + c n x n , (-1) n c 0 , (-1) n+1 B)
Algorithm 1 is based on the recursion c n-k = -1 k k j=1 c n-k+j Tr(A j ): we compute a sequence of matrices (stored in the variable B) through repeated multiplication by A, and in each step extract a trace. The determinant and the adjugate 1 Indeed, the present author also reinvented the same method. The first version of this note presented it simply as a "baby-step giant-step Faddeev-Leverrier algorithm", mentioning Berkowitz's attribution of the idea to S. Winograd, but without citing the concrete algorithm already given in [START_REF] Preparata | An improved parallel processor bound in fast matrix inversion[END_REF]. The author is grateful to Eric Schost for pointing out the relevant prior art. matrix appear as byproducts of this process: (-1) n det(A) is the last coefficient c 0 , and (-1) n+1 adj(A) is the penultimate entry in the matrix sequence. 2It is easy to see that Algorithm 1 performs O(n) matrix multiplications and O(n 2 ) additional arithmetic operations. The condition on the characteristic of R ensures that we can divide exactly by each k, i.e. (xk)/k = x holds for x ∈ R, k ≤ n.

Algorithm 2 Preparata-Sarwate algorithm (slightly modified) Input: A ∈ R n×n with the same conditions as in Algorithm 1 Output: (p A (x), det(A), adj(A))

1: m ← ⌊ √ n⌋ 2: Precompute the matrices A 1 , A 2 , A 3 , . . . , A m 3: Precompute the traces t k = Tr(A k ) for k = 1, . . . , m 4: c n = 1, B ← I, k ← 1 5: while k ≤ n -1 do 6: m ← min(m, n -k) 7: c n-k ← -1 k Tr(A, B) 8: for j ← 1, 2, . . . , m -1 do 9:
c n-k-j ← Tr(A j+1 , B) ⊲ Using precomputed power of A 10:

for i ← 0, 1, . . . , j -1 do 11: c n-k-j ← c n-k-j + t j-i c n-k-i 12:
end for 13:

c n-k-j ← c n-k-j /(-k -j) 14:
end for

15:

B ← A m B ⊲ Using precomputed power of A 16:

for j ← 0, 1, . . . , m -1 do 17:

B ← B + c n-k-j A m-j-1 ⊲ Using precomputed power, or A 0 = I 18:
end for

19: k ← k + m 20: end while 21: c 0 ← -1 n Tr(A, B) 22: return (c 0 + c 1 x + . . . + c n x n , (-1) n c 0 , (-1) n+1 B)
Algorithm 1 computes a sequence of n matrices but only extracts a small amount of unique information (the trace) from each matrix. In such a scenario, we can often save time using a baby step giant-step approach in which we only compute O(

√ n) products explicitly (see [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF][START_REF] Brent | Fast algorithms for manipulating formal power series[END_REF][START_REF] Stuart | On computing the determinant in small parallel time using a small number of processors[END_REF][START_REF] Johansson | A fast algorithm for reversion of power series[END_REF] for other examples of this technique).

Preparata and Sarwate [START_REF] Preparata | An improved parallel processor bound in fast matrix inversion[END_REF] improve Algorithm 1 by choosing m ≈ √ n and precomputing the powers A 1 , A 2 , . . . , A m and A m , A 2m , A 3m , . . .. The key observation is that we can compute Tr(AB) using O(n 2 ) operations without forming the complete matrix product AB, by simply evaluating the dot products for the main diagonal of AB. We denote this product trace operation by Tr(A, B). Algorithm 2 presents a version of this method. It is clear from inspection that this version performs roughly m + n/m ≈ 2 √ n matrix multiplications of size n × n, and O(n 3 ) arithmetic operations in the remaining steps.

Algorithm 2 is a small modification of the Preparata-Sarwate algorithm as it is usually presented. Instead of computing the giant-step powers A m , A 2m , A 3m , . . . explicitly, we expand the loop in Algorithm 1 to group m iterations, using multiplications by A m to update the running sum over both the powers and their traces. This version performs essentially the same number of operations while using half the amount of memory.

As an observation for implementations, the matrix-matrix multiplications and product traces are done with invariant operands that get recycled O( √ n) times.

This can be exploited for preconditioning purposes, for instance by packing the data more efficiently for arithmetic operations. We also note that the optimal m may depend on the application, and a smaller value will reduce memory consumption at the expense of requiring more matrix multiplications.

Applicability and performance evaluation

When, if ever, does it make sense to use Algorithm 2? If R is a field, then the determinant, adjugate and characteristic polynomial can be computed in O(n ω ) operations or O(n 3 ) classically allowing divisions [START_REF] Dumas | Efficient computation of the characteristic polynomial[END_REF][START_REF] Pernet | Faster algorithms for the characteristic polynomial[END_REF][START_REF] Abdeljaoued | Efficient algorithms for computing the characteristic polynomial in a domain[END_REF][START_REF] Neiger | Deterministic computation of the characteristic polynomial in the time of matrix multiplication[END_REF]. Most commonly encountered rings are integral domains, in which case we can perform computations in the fraction field and in some cases simply clear denominators. This does not automatically render an O(n 3.5 ) algorithm obsolete, since naively counting operations may not give the whole picture. Nevertheless, we can immediately discard some applications:

• For computing over R and C in ordinary floating-point arithmetic, the Faddeev-Leverrier algorithm is slower and far less numerically stable than textbook techniques such as reduction to Hessenberg form and Gaussian elimination with O(n 3 ) or better complexity, as we already noted in the introduction. • For finite fields, classical O(n 3 ) methods using divisions have no drawbacks, and linear algebra with O(n 2.81 ) Strassen complexity is well established [START_REF] Dumas | Computational linear algebra over finite fields[END_REF]. Over rings with small characteristic, the applicability of Algorithm 2 is in any case limited due to the integer divisions. Generally speaking, division-free or nearly division-free algorithms are interesting for rings and fields R where dividing recklessly can lead to coefficient explosion (for example, Q) or in which testing for zero is problematic (for example, exact models of R). The optimal approach in such situations is usually to avoid computing directly in R or its fraction field, for example using modular arithmetic and interpolation techniques, but such indirect methods are more difficult to implement and must typically be designed on a case by case basis. By contrast, Algorithm 2 is easy to use anywhere. We will now look at some implementation experiments.

3.1. Integers. For exact linear algebra over Z and Q, the best methods are generally fraction-free versions of classical algorithms (such as the Bareiss version of Gaussian elimination) for small n, and multimodular or p-adic methods for large n (see for example [START_REF] Dumas | Efficient computation of the characteristic polynomial[END_REF][START_REF] Pernet | Faster algorithms for the characteristic polynomial[END_REF]). We do not expect Algorithm 2 to beat those algorithms, but it is instructive to examine its performance. Table 1 shows timings for computing a determinant, inverse or characteristic polynomial of an n × n matrix over Z with random entries in -10, . . . , 10, using the following algorithms:

• FFLU: fraction-free LU factorization using the Bareiss algorithm. • FFLU2: as above, but using the resulting decomposition to compute A -1 (equivalently determining adj(A)) by solving AA -1 = I. • ModDet: a multimodular algorithm for the determinant.

• ModInv: a multimodular algorithm for the inverse matrix.

• ModCP: a multimodular algorithm for the characteristic polynomial.

• Berk: the Berkowitz algorithm for the characteristic polynomial.

• Alg1: original Faddeev-Leverrier algorithm, Algorithm 1.

• Alg2: modified Preparata-Sarwate algorithm, Algorithm 2. We implemented Alg1 and Alg2 on top of Flint [START_REF] Hart | Fast Library for Number Theory: An introduction[END_REF], while all the other tested algorithms are builtin Flint methods.

3.1.1. Observations. Alg2 clearly outperforms both Alg1 and Berk for large n, making it the best algorithm for computing the characteristic using direct arithmetic in Z (the modular algorithm is, as expected, superior). Alg2 is reasonably competitive for computing the inverse or adjugate matrix, coming within a factor 2-3× of FFLU2 and ModInv in this example. For determinants, the gap to the FFLU algorithm is larger, and the modular determinant algorithm is unmatched.

Number fields.

Exact linear algebra over algebraic number fields Q(a) is an interesting use case for division-free algorithms since coefficient explosion is a significant problem for classical O(n 3 ) algorithms. As in the case of Z and Q, modular algorithms are asymptotically more efficient than working over Q(a) directly, but harder to implement. Here we compare the following algorithms:

• Sage: the charpoly method in SageMath [START_REF]The Sage Developers[END_REF], which implements a specialpurpose algorithm for cyclotomic fields based on modular computations and reconstruction using the Chinese remainder theorem. • Hess: Hessenberg reduction for the characteristic polynomial • Dani: Danilevsky's algorithm for the characteristic polynomial.

• LU: LU factorization to compute the determinant.

• FFLU: fraction-free LU factorization using the Bareiss algorithm.

• LU2 and FFLU2: as above, but using the resulting decomposition to compute A -1 (equivalently determining adj(A)) by solving AA -1 = I. • Berk (Berkowitz), Alg1 and Alg2 as in the previous section. With the exception of the Sage function, we implemented all the algorithms using Antic [START_REF] Hart | ANTIC: Algebraic number theory in C[END_REF] for number field arithmetic and Flint for other operations. We perform fast matrix multiplication by packing number field elements into integers and multiplying matrices over Z via Flint. Our implementations of LU, LU2, Alg1 and Alg2 benefit from matrix multiplication while Hess, Dani, FFLU, FFLU2 and Berk do not. The benchmark is therefore not representative of the performance that ideally should be achievable with these algorithms, although it is fair in the sense that the implementation effort for Alg1 and Alg2 was minimal while the other algorithms would require much more code to speed up using block strategies.

Table 2 compares timings for two kinds of input: random matrices over a fixed cyclotomic field, and discrete Fourier transform (DFT) matrices which have special structure. Choosing cyclotomic fields allows us to compare with the dedicated algorithm for characteristic polynomials in Sage; the corresponding method for generic number fields in Sage is far slower. All the other algorithms make no assumptions about the field.

3.2.1.

Observations. There are no clear winners since there is a complex interplay between operation count, multiplication algorithms, matrix structure and coefficient growth. Modular algorithms are the best solution in general for large n, but implementations for number fields are complex and less readily available in current software than for Z and Q. Among the non-modular algorithms, the O(n 3 ) division-heavy Hessenberg and Danilevsky algorithms are nearly useless due to coefficient explosion for generic input, but both perform well on the DFT matrix. The LU and FFLU algorithms have more even performance but alternate with each other for the advantage depending on the matrix. Alg2 has excellent average performance for the determinant, characteristic polynomial as well as the adjugate matrix considering the large variability between the algorithms for different input. It is highly competitive for computing the inverse or adjugate of the random matrix.

3.3. Ball arithmetic. Division-free algorithms are useful when computing rigorously over R and C in interval arithmetic or ball arithmetic. The reason is that we cannot test whether elements are zero, so algorithms like Gaussian elimination and Hessenberg reduction fail when they need to branch upon zero pivot elements or zero vectors. Although zeros will not occur for random input, they are likely to occur for structured matrices arising in applications. A posteriori verification of approximate numerical solutions or perturbation analysis is in principle the best workaround [START_REF] Rump | Verification methods: Rigorous results using floating-point arithmetic[END_REF], but it is sometimes useful to fall back to more direct division-free methods, especially when working in very high precision.

Table 3 shows timings for computing a determinant or characteristic polynomial with 100-digit accuracy using the following algorithms implemented in ball arithmetic. The input is taken to be an n × n real matrix with uniformly random entries in [0, 1]. For this experiment, we only focus on the determinant and characteristic polynomial (the conclusions regarding matrix inversion would be similar to those regarding the determinant).

• Hess: Hessenberg reduction using Gaussian elimination.

• Hess2: Hessenberg reduction using Householder reflections.

• LU: LU factorization using Gaussian elimination.

• LU2: approximate computation of the determinant using LU factorization followed by a posteriori verification.

• Eig: approximate computation of the eigenvalues using the QR algorithm followed by a posteriori verification and reconstruction of the characteristic polynomial from its roots. • Berk (Berkowitz), Alg1 and Alg2 as in the previous section. All algorithms were implemented in Arb [START_REF] Johansson | Arb: Efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF] which uses the accelerated dot product and matrix multiplication algorithms described in [START_REF] Johansson | Faster arbitrary-precision dot product and matrix multiplication[END_REF]. The LU, LU2, Alg1 and Alg2 implementations benefit from fast matrix multiplication while Hess, Hess2, Eig and Berk do not.

The methods LU2 and Eig are numerically stable: the output balls are precise to nearly full precision for well-conditioned input. All other algorithms are unstable in ball arithmetic and lose O(n) digits of accuracy. At least on this example, the rate of loss is almost the same for Hess2, Berk, Alg1 and Alg2, while LU is more stable and Hess and Dani are less stable. To make the comparison meaningful, we set the working precision (shown in the table) to an experimentally determined value so that all algorithms enclose the determinant with around 100 digits of accuracy.

3.3.1.

Observations. For computing the characteristic polynomial in high-precision ball arithmetic, it seems prudent to try Hessenberg reduction and fall back to a division-free algorithm when it fails due to encountering a zero vector. The Berkowitz algorithm is the best fallback for small n, while Alg2 wins for large n (n ≈ 50, although the cutoff will vary). On this particular benchmark, Alg2 runs only about 4× slower than Hessenberg reduction, making it an interesting onesize-fits-all algorithm. The verification method (Eig) gives the best results if the precision is constrained, but is far more expensive than the other methods.

For computing the determinant alone, all the division-free methods are clearly inferior to methods based on LU factorization in this setting. The only advantage of the division-free algorithms is that they are foolproof while LU factorization requires some attention to implement correctly.

Better methods for computing the characteristic polynomial in ball arithmetic or interval arithmetic are surely possible. For the analogous problem of computing the characteristic polynomial over Q p , see [START_REF] Caruso | Characteristic polynomials of p-adic matrices[END_REF].

3.4. Polynomial quotient rings. At first glance Algorithm 2 seems to hold potential for working over multivariate polynomial quotient rings. Such rings need not be integral domains and division can be very expensive (requiring Gröbner basis computations). Unfortunately, in most examples we have tried, Algorithm 2 performs worse than both the Berkowitz algorithm and Algorithm 1, presumably because repeated multiplication by the initial matrix A is much cheaper than multiplication by a power of A which generally will have much larger entries. There may be special classes of matrices for which the method performs well, however.

Discussion

We find that the Preparata-Sarwate algorithm often outperforms the Berkowitz algorithm in practice, in some circumstances even being competitive with O(n 3 ) algorithms. We can therefore recommend it for more widespread adoption.

Galil and Pan [START_REF] Galil | Parallel evaluation of the determinant and of the inverse of a matrix[END_REF][START_REF] Abdeljaoued | Méthodes matricielles-Introduction à la complexité algébrique[END_REF] have further refined the Preparata-Sarwate algorithm to eliminate the O(n 3 ) complexity term which dominates asymptotically with a sufficiently fast matrix multiplication algorithm. We have not tested this method since those O(n 3 ) operations are negligible in practice.

Table 1 .

 1 Time in seconds to compute characteristic polynomial (C), determinant (D), adjugate/inverse (A) of an n×n matrix over Z with random elements in -10, . . . , 10, using various algorithms.

	n	FFLU	ModDet	FFLU2	ModInv ModCP	Berk	Alg1	Alg2
		D	D	DA	A	CD	CD	CDA	CDA
	10	0.0000060 0.000021 0.000015 0.00012 0.000016 0.000035 0.000015 0.000030
	20	0.000036	0.000078	0.00043	0.00096	0.00011	0.0010	0.00086	0.00061
	50	0.0023	0.0012	0.011	0.016	0.0039	0.048	0.052	0.017
	100	0.039	0.0068	0.14	0.18	0.055	0.84	1.1	0.22
	200	0.64	0.044	2.3	1.8	0.89	16	27	4.1
	300	3.4	0.15	13	9.0	4.6	94	174	20
	400	12	0.38	45	22	15	321	696	66
	500	32	0.77	127	52	37	900	2057	150

Table 2 .

 2 Time in seconds to compute characteristic polynomial (C), determinant (D), adjugate/inverse (A) of a matrix over a cyclotomic number field.

	n	Sage	Hess	Dani	LU	FFLU	LU2	FFLU2 Berk		Alg1	Alg2
		CD	CD	CD	D	D	DA	DA	CD		CDA	CDA
		Input: n × n matrix over Q(ζ 20 ), entries k (p/q)ζ k 20 , random |p| ≤ 10, 1 ≤ q ≤ 10.
	10	0.038	0.31	0.16	0.024	0.0059	0.21	0.11	0.010 0.0073	0.010
	20	0.12	19	6.7	0.22	0.067	2.6	1.4	0.28		0.15	0.16
	30	0.39	200	67	0.93	0.31	12	6.8	2.0		1.1	0.8
	40	1.1		353	2.8	0.9	37	22	7.5		3.7	2.6
	50	1.9			7.0	2.3	88	56	22		8.7	5.7
	60	3.4			15	4.7	182	119	54		19	12
	70	5.1			29	8.6	337	230	114		39	22
	80	7.5			53	15	581	409	208		67	34
	90	11			89	24			397		144	54
	100	15			144	41			608		670	130
	120	24			322	83			1439		3013	420
									(i-1)(j-1) n	.
	10	0.010 0.0018 0.0016 0.00017 0.00022 0.00076 0.0014 0.00061 0.00075 0.00059
	20	0.039 0.0019 0.0024 0.0017 0.0046 0.0071	0.038	0.020	0.020	0.0070
	50	1.3	0.17	0.13	0.065	0.80	0.28	6.0	8.2		2.0	0.49
	100	22	5.4	22	0.89	43	5.3	335	803		223	29
	150	78	22	7.9	4.4	214	19	1423	7259		933	138
	200	333	1928	140	31	1655	192					1687

Input: n × n DFT matrix over Q(ζn), entries A i,j = ζ

Table 3 .

 3 Time in seconds to compute characteristic polynomial (C), determinant (D), adjugate/inverse (A) of a random n × n matrix in real ball arithmetic. The respective algorithms were run with 333 + p bits of precision, with p chosen to give roughly 100digit output accuracy.

	n	Hess	Hess2	Dani	Eig	LU	LU2	Berk	Alg1	Alg2
		CD	CD	CD	CD	D	D	CD	CDA	CDA
	10	0.00021 0.00038 0.00022 0.017 0.000068 0.00023 0.00080 0.0010 0.00078
	20	0.0021	0.0032	0.0020	0.18	0.00052	0.0015	0.0039	0.017	0.0092
	50	0.045	0.057	0.048	4.6	0.0078	0.019	0.22	0.61	0.21
	100	0.64	0.69	0.61	56	0.062	0.15	6.3	9.7	2.4
	150	3.5	3.5	3.0	245	0.23	0.44	52	52	11
	200	12	11	10		0.59	1.0	224	176	34
	250	31	29	25		1.4	1.9	687	460	73
	300	66	59	53		2.5	3.2	1804	1075	160
	350	135	115	110		4.4	5.0	4033	2107	306
	p	10n	6n	10n	0	n	0	6n	6n	6n

The code can be tightened assuming n ≥ 2, in which case the line before the start of the loop can be changed to {cn = 1, c n-1 = -Tr(A), B ← A + c n-1 I, k ← 2}. We can change the loop condition to k ≤ n and remove the line after the loop if we omit returning adj(A).

An interesting problem is whether it is possible to design a division-free or nearly division-free algorithm with better than O(n 4 ) classical complexity that minimizes the observed problems with growing entries, particularly in multivariate rings.


Meritxell Genovart
email: m.genovart@csic.es

Olivier O Gimenez

Albert Bertolero

Rémi Choquet

Daniel Oro

Roger Pradel

Decrease in social cohesion in a colonial seabird under a perturbation regime

Keywords: Audouin's gull, colonial species, decision-making, non-linear response, perturbations, resilience, social cohesion, social network analyses

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

Introduction

Ecosystems are subject to perturbations, both natural and human induced, affecting individuals, populations and communities. When they are strong or are maintained through time, these perturbations may cause a shift between stable states at the level of both individual and population and even lead to population collapses and extinctions [START_REF] Dai | Relation between stability and resilience determines the performance of early warning signals under different environmental drivers[END_REF][START_REF] Dakos | Resilience indicators: prospects and limitations for early warnings of regime shifts[END_REF] .

Understanding how individuals and populations will respond to these perturbations is critical both from a 'pure' ecological standpoint and also from an applied point of view to mitigate the effects of global change [START_REF] Colchero | The diversity of population responses to environmental change[END_REF][START_REF] Coulson | Data from: Modeling adaptive and nonadaptive responses of populations to environmental change[END_REF][START_REF] Donohue | Navigating the complexity of ecological stability[END_REF] . Population dynamics may be directly affected by these perturbations through a decrease in demographic parameters such as survival or fecundity, or by a change, immediate or delayed, on individual behaviour, such as an increase in dispersal [START_REF] Fernandez-Chacon | When to stay, when to disperse and where to go: survival and dispersal patterns in a spatially structured seabird population[END_REF] . We define population resilience as the maximal pulse perturbation a population can tolerate or absorb without going extinct [START_REF] Dai | Relation between stability and resilience determines the performance of early warning signals under different environmental drivers[END_REF][START_REF] Holling | Resilience and Stability of Ecological Systems[END_REF] . In social animals, social behavioural processes, such as information sharing and decisionmaking, add another dimension to understanding the resilience of populations facing perturbations. For instance, the amount of social information can be enhanced not only by positive density-dependence, but also by social cohesion [START_REF] Barrett | Taking sociality seriously: the structure of multi-dimensional social networks as a source of information for individuals[END_REF][START_REF] Centola | How Behavior Spreads: the Science of Complex Contagions[END_REF][START_REF] Firth | Considering Complexity: Animal Social Networks and Behavioural Contagions[END_REF][START_REF] Kerth | Bats are able to maintain long-term social relationships despite the high fission-fusion dynamics of their groups[END_REF] . Social cohesion favor the exchange of private information and consequently reduce uncertainty in resource acquisition (e.g. shelter against predators, food, mates) or in decision-making in the face of disturbances, such as dispersal to non-perturbed or less perturbed sites [START_REF] Fernandez-Chacon | When to stay, when to disperse and where to go: survival and dispersal patterns in a spatially structured seabird population[END_REF][START_REF] Rosenthal | Revealing the hidden networks of interaction in mobile animal groups allows prediction of complex behavioral contagion[END_REF][START_REF] Snijders | Animal Social Network Theory Can Help Wildlife Conservation[END_REF][START_REF] Webber | An evolutionary framework outlining the integration of individual social and spatial ecology[END_REF] . Thus, the structure of a group may affect social interactions, information transfer, and collective decisions [START_REF] Sueur | Social Interaction in Animals: Linking Experimental Approach and Social Network Analysis[END_REF] . Some recent studies also show that spatial cohesion may be risk sensitive, and that responses may differ depending on the risk and the social group [START_REF] Labarge | Reactive and pre-emptive spatial cohesion in a social primate[END_REF] . However, little is known about the effects of environmental perturbations on the cohesion of social groups in empirical studies of social animals [START_REF] Oro | Perturbation, social feedbacks, and population dynamics in social animals[END_REF] .

The analysis of social relationships in animal populations may include a range of social dynamics, from simple and ephemeral contacts, to permanent and strong bonds between individuals [START_REF] Kerth | Bats are able to maintain long-term social relationships despite the high fission-fusion dynamics of their groups[END_REF][START_REF] Firth | Experimental manipulation of avian social structure reveals segregation is carried over across contexts[END_REF][START_REF] Genton | How Ebola impacts social dynamics in gorillas: a multistate modelling approach[END_REF][START_REF] Leu | Environment modulates population social structure: experimental evidence from replicated social networks of wild lizards[END_REF] . Coloniality is a life-history strategy where individuals show a social link among conspecifics by breeding in large and dense groups [START_REF] Brown | The ecology and evolution of colony-size variation[END_REF][START_REF] Rolland | The evolution of coloniality in birds in relation to food, habitat, predation, and life-history traits: a comparative analysis[END_REF] . However, many colonial species are philopatric, thus this link may not necessarily reflect individual social ties but a shared tendency to breed in the same birthplace. This tendency may result from the need to share information about resources, especially when they are patchy and more unpredictable, or it may result from the advantages of social defence against predators [START_REF] Clode | Colonially breeding seabirds: Predators or prey?[END_REF][START_REF] Hoogland | Advantages and Disadvantages of Bank Swallow (Riparia riparia) Coloniality[END_REF] . A challenge lies in disentangling whether annual association between individuals is only due to philopatry, or also due to the existence of a social ties within groups of individuals over time [START_REF] Shizuka | Across-year social stability shapes network structure in wintering migrant sparrows[END_REF] . Social ties between neighbouring pairs in breeding colonies are rarely considered in behavioural and ecological studies [START_REF] Brandl | Wild zebra finches that nest synchronously have long-term stable social ties[END_REF] and, if true, such associations may suggest the evolution of social cohesion for exploiting the evolutionary advantages of social living (including social information sharing) for individual fitness prospects. Social network theory, which originated in sociology to study human relationships and social organization [START_REF] Centola | How Behavior Spreads: the Science of Complex Contagions[END_REF][START_REF] Moreno | Who shall survive?: A new approach to the problem of human interrelations[END_REF][START_REF] Scott | Social Network Analysis[END_REF] now provides both a conceptual framework and the analytical tools to explore social cohesion and social processes in animal populations [START_REF] Croft | Exploring Animal Social Networks[END_REF][START_REF] Farine | Structural trade-offs can predict rewiring in shrinking social networks[END_REF][START_REF] Farine | Constructing, conducting and interpreting animal social network analysis[END_REF][START_REF] Ward | Sociality: The Behaviour of Group-Living Animals[END_REF][START_REF] Whitehead | Analyzing Animal Societies Quantitative Methods for Vertebrate Social Analysis[END_REF] .

Network theory is now being simultaneously developed in a number of fields, including statistical physics, sociology, molecular biology, and computer science. As a result, the field is changing at a rapid pace. While not all developments can or should be applied toward the study of animal societies [START_REF] James | Potential banana skins in animal social network analysis[END_REF] , this rush of novel ideas from outside disciplines is enriching behavioural ecology [START_REF] Hasenjager | Chapter Three -Social Network Analysis in Behavioral Ecology[END_REF] .

To assess the existence of social cohesion and its dynamics under perturbations, we studied an ecological system involving a colonial, social vertebrate (the Audouin's gull Larus audouinii) living in a site experiencing a shift to a perturbed regime. Interestingly from a social point of view, the species breeds aggregated in spatially-discrete patches within large colonies. Each breeding season, some patches go extinct and some are colonized [START_REF] Payo-Payo | Predator arrival elicits differential dispersal, change in age structure and reproductive performance in a prey population[END_REF] , forcing individuals to breed in patches different from the ones they were born in or they bred in the previous year. These colonization-extinction processes may allow us to disentangle whether social aggregation among individuals is an annual random association, or it rather results from social cohesion among individuals over the years.

An extensive long-term individual monitoring program has been carried out since 1988 at the Ebro Delta, including the main breeding site for the species, the Punta de la Banya [START_REF] Genovart | Immature survival, fertility, and density dependence drive global population dynamics in a long-lived species[END_REF] . At the study site, population dynamics has undergone different phases: an initial growing phase after site colonization, a stable phase of dynamic equilibrium, and a final transition phase to population collapse [START_REF] Oro | Perturbation, social feedbacks, and population dynamics in social animals[END_REF][START_REF] Almaraz | Size-mediated non-trophic interactions and stochastic predation drive assembly and dynamics in a seabird community[END_REF] (Figure 1). This collapse was due to the arrival of terrestrial predators, and led this colony to hold from 70% to only 3% of the total world population in only a decade (32% mean annual decline) [START_REF] Payo-Payo | Predator arrival elicits differential dispersal, change in age structure and reproductive performance in a prey population[END_REF][START_REF] Genovart | Immature survival, fertility, and density dependence drive global population dynamics in a long-lived species[END_REF] (Figure 1). Most predators were foxes, but also badgers and other mesocarnivores. Predators invaded the site likely due to their increasing densities in recent decades [START_REF] Martínez-Abraín | refuge abandonment' and spread of fearless behavior in a reconciling world[END_REF] and the attractiveness of the site in terms of food availability and lack of competition. The perturbation regime caused changes in the spatial distribution of patches at the site, changes in age structure, decrease in fecundity and a progressively increase of dispersal to other sites [START_REF] Fernandez-Chacon | When to stay, when to disperse and where to go: survival and dispersal patterns in a spatially structured seabird population[END_REF][START_REF] Payo-Payo | Predator arrival elicits differential dispersal, change in age structure and reproductive performance in a prey population[END_REF] . The response of this population to predators has been not immediate probably due to strong philopatry, high site-suitability inertia and social behavioural processes, such as conspecific attraction [START_REF] Oro | Perturbation, social feedbacks, and population dynamics in social animals[END_REF] . This raise in dispersal was caused by social processes, as social copying [START_REF] Oro | Perturbation, social feedbacks, and population dynamics in social animals[END_REF] , however it remains to assess how social cohesion among individuals, if occurred, was affected by dispersal processes. One possibility is that dispersal would broke social cohesion by individual heterogeneity in the willingness to disperse [START_REF] Shizuka | How demographic processes shape animal social networks[END_REF] . In contrast, social cohesion can be maintained over time when dispersal occurs collectively at the scale of social groups to the same sites [START_REF] Shizuka | Across-year social stability shapes network structure in wintering migrant sparrows[END_REF][START_REF] Francesiaz | Familiarity drives social philopatry in an obligate colonial breeder with weak interannual breeding-site fidelity[END_REF] .

Taking advantage of the long-term monitoring of this long-lived species, the knowledge of its population dynamics, and the use of tools recently developed in the Social Network Analysis (SNA) framework, we specifically addressed the following questions: 1) is there any long-term social ties between individuals breeding in the same patch? and 2) have perturbations, in this case a perturbation regime, affected social cohesion? We finally discuss the role and consequences of social cohesion in population dynamics and resilience in social species.

Results

We analysed a total of 1,610,922 dyadic interactions during the first period (2002-2011) and 368,142 during the second period (2012-2017).

When assessing the social ties with the contingency table approach during the period of stability, the assumption that breeding aggregations in Audouin's gull were at random was rejected, and those individuals that bred together during the sub-period 2002-2006 had a higher probability of breeding together during the sub-period 2007-2011 (𝜒 = 64.685, P < 0.0001). When randomly reducing sample size of the data set, results were still statistically significant in more than 95% of the cases (1000 randomizations).

Accordingly, when assessing the social ties with dependent regression terms in the AME function, we showed that the probability of breeding together during the second sub-period (2007-2011) depended on whether they have bred previously together in the first sub-period (2002-2006), with a statistically significant coefficient of regression parameter (Table 1; Figure 2). When analysing this data set with permutated data, we concluded that there was indeed a non-random association of individuals within the patches, with our statistic being among the 5% extreme values.

When we analysed the social ties during the transition to collapse phase, we observed that the probability of breeding together during the period 2012-2017 did not depend on whether they have bred together the five previous years (χ 2 =1.957, p-value = 0.162) and we could not reject the hypothesis of a random association between individuals. In addition, the SNA approach showed that breeding aggregations in Audouin's gull during the transition phase did not depend on whether they have bred together the five previous years (Table 1; Figure 2).

Stable period

Transition to collapse period group have been identified [START_REF] Ginsberg | Measuring association between individuals or groups in behavioural studies[END_REF][START_REF] Cairns | A comparison of association indices[END_REF] .

Discussion

By studying a particular ecological system of a colonial long-lived species that experiences a perturbation regime, we showed that social ties among individuals persist over unperturbed years and that perturbations may decrease social cohesion in animal populations.

The characteristic breeding behaviour of the study species that aggregates in patches that change annually, allowed us to show that individuals do not annually breed aggregated at random but rather there is some group stability, with individuals establishing social ties that persist over time. Our study system resembles what it was recorded for Slender-billed gulls (Chroicocephalus genei), a colonial breeder with weak inter-annual breeding-site fidelity: some individuals bred together across breeding seasons and some social groups showed tenacity despite the colony often moving every year [START_REF] Francesiaz | Familiarity drives social philopatry in an obligate colonial breeder with weak interannual breeding-site fidelity[END_REF] . Group stability can emerge as a product of network self-organization, but may provide the necessary conditions for the evolution of other social processes [START_REF] Cantor | Simple foraging rules in competitive environments can generate socially structured populations[END_REF][START_REF] Cantor | Animal social networks: revealing the causes and implications of social structure in ecology and evolution[END_REF] . Our results would support the idea that social aggregation during breeding would provide other advantages than the mere defence against predators [START_REF] Anderson | Predator Behavior Favors Clumped Nesting in an Oceanic Seabird[END_REF][START_REF] Oro | Colonial Seabird Nesting in Dense and Small Sub-Colonies: An Advantage against Aerial Predation[END_REF] , such as social information sharing (e.g. [START_REF] Gil | Social Information Links Individual Behavior to Population and Community Dynamics[END_REF][START_REF] Lewanzik | Insectivorous bats integrate social information about species identity, conspecific activity and prey abundance to estimate cost-benefit ratio of interactions[END_REF] ). Social information sharing is crucial for decision-making in risky behaviours, such as dispersal, and previous studies showed that the perturbed regime in our study site caused dispersal to other sites, including colonization of new habitats [START_REF] Payo-Payo | Predator arrival elicits differential dispersal, change in age structure and reproductive performance in a prey population[END_REF][START_REF] Payo-Payo | Colonisation in social species: the importance of breeding experience for dispersal in overcoming information barriers[END_REF] . In our case study, sociality may have played a major role driving dispersal and thus population dynamics, both during the exponential growth after colonization and the collapse after the perturbation regime. This idea is also reinforced with a mechanistic dynamical model that shows that population dynamics of Audouin's gulls at the study site can only be explained by dispersal runaway caused by social copying [START_REF] Oro | Perturbation, social feedbacks, and population dynamics in social animals[END_REF] .

The importance of social information compared to private information increases under perturbations, even when the quality of social information does not increase compared to a non-perturbed regime [START_REF] Arganda | A common rule for decision making in animal collectives across species[END_REF][START_REF] Pérez-Escudero | Adversity magnifies the importance of social information in decision-making[END_REF] . For instance, Maldonado et al. [START_REF] Maldonado-Chaparro | Transient LTRE analysis reveals the demographic and trait-mediated processes that buffer population growth[END_REF] show that experimental disturbances applied to a social bird may impact its foraging efficiency by social instability caused by the split of social groups. In colonial birds, breeding failure, which is a proxy of environmental stress, may trigger splitting of the social groups (e.g. [START_REF] Francesiaz | Familiarity drives social philopatry in an obligate colonial breeder with weak interannual breeding-site fidelity[END_REF] ). At demographic level, the alteration of social network structure may affect the behaviour of populations. For instance, under stress conditions, sociality may operate through feedback loops such as social copying for dispersal, causing non-linear population dynamics and playing a critical role on the resilience of populations (e.g. [START_REF] Oro | Perturbation, social feedbacks, and population dynamics in social animals[END_REF] ). We showed here that after a perturbation, not only the number of individuals in the population may decrease (by increased mortality or dispersal) but also its social cohesion, likely reducing but also altering the information transfer within the social network composed by those individuals that remain in the site where perturbation occurs. Among other demographic processes, dispersal may alter social connections of both individuals remaining and those dispersing, with consequences for social network structure [START_REF] Shizuka | How demographic processes shape animal social networks[END_REF] . The perturbation regime suffered by the study population has likely triggered a social transition [START_REF] Pruitt | Social tipping points in animal societies[END_REF] in collective behaviour from philopatric to dispersal and with the fast diffusion of innovations such as the colonization of harbours by large number of individuals, a habitat safe from predators never occupied before [START_REF] Payo-Payo | Colonisation in social species: the importance of breeding experience for dispersal in overcoming information barriers[END_REF] . Previous studies have shown that responses of populations to perturbations may also depend on individual personalities in the population [START_REF] Dall | The behavioural ecology of personality: consistent individual differences from an adaptive perspective[END_REF][START_REF] Doering | Social tipping points in animal societies in response to heat stress[END_REF][START_REF] Wolf | Life-history trade-offs favour the evolution of animal personalities[END_REF] . For example, dispersers are different from non-dispersing individuals for a suite of phenotypic traits, including their behavioural profile [START_REF] Clobert | Informed dispersal, heterogeneity in animal dispersal syndromes and the dynamics of spatially structured populations[END_REF][START_REF] Cote | Personality-dependent dispersal: characterization, ontogeny and consequences for spatially structured populations[END_REF][START_REF] Fogarty | Social Personality Polymorphism and the Spread of Invasive Species: A Model[END_REF] . Heterogeneities in personalities for dispersal decision-making may have also played a role in our studied population, with most individuals dispersing to other sites after a period of disturbance, while some individuals remaining philopatric. This change may have also further consequences for social network stability, as performance in social groups may improve with heterogeneity in individual personalities [START_REF] Fogarty | Social Personality Polymorphism and the Spread of Invasive Species: A Model[END_REF][START_REF] O'shea-Wheller | Variability in individual assessment behaviour and its implications for collective decision-making[END_REF] .

Our study opens new research avenues about resilience of social populations under perturbations; if perturbations affect social cohesion and heterogeneity in personalities in the population, we may wonder whether this population would be equally resilient to future perturbations. Additionally, in our study population, sociality seemed to operate not right after the first perturbation episode but after cumulative maintained perturbation 50

; would the type of perturbation, either pulse or in regime [START_REF] Nimmo | Vive la résistance: reviving resistance for 21st century conservation[END_REF] influence the response of social groups? We have shown here that a regime perturbation may decrease social cohesion in animal populations, but further studies should be carried out to improve our understanding on the demographic consequences of the breaking down of social ties under perturbations for population dynamics and resilience in social species.

Material and Methods

Study species and study area

The Audouin's gull is a long-lived seabird with more than 80% of the global population breeding in the western Mediterranean (http://www.iucnredlist.org/details/22694313/0; 37 . The species was critically endangered until the early 80's, when it colonized a new site, the Punta de la Banya in the Ebro Delta (Figure 2). Here, the large availability of both suitable breeding habitat and food resulted in a rapid and exponential growth, ending with the site holding more than 70% of the total world population in 2006. The global population dynamics was mainly driven by this colony and after the exponential growth, the species was downgraded to a conservation category of "least concern" [START_REF] Larus | BirdLife International: The IUCN Red List of Threatened Species 2018: e[END_REF] . However, the Punta de la Banya colony is now collapsing and even if the species is colonizing new sites, the global population is decreasing at a 5% annual rate [START_REF] Genovart | Immature survival, fertility, and density dependence drive global population dynamics in a long-lived species[END_REF] . In 1997, first carnivores (mainly foxes, but also badgers, beech martens and least weasels) arrived at Punta de la Banya, and since then the site has been perturbed by the presence of carnivores.

Annual censuses of breeding pairs at every patch within colonies at the Ebro Delta area have been carried out since colonization in 1981 to 2017 (Figure 1, 3 and Table S1). In the Ebro Delta there are three colonies: Punta de la Banya, colonized in 1981 and occupied throughout the study period, Sant Carles de la Ràpita harbour, occupied since 2011 to 2015, and Sant Antoni, occupied from 2013 up to now (Figure 3, Figure S1).

Within colonies, individuals are distributed in patches [START_REF] Genovart | Offspring sex ratios in subcolonies of Audouin's gull, Larus audouinii , with differential breeding performance[END_REF] . As patch location may change from one year to another, we annually geolocalized, mapped and defined the breeding patches (Figure 3 and S2). The species is monogamous, there is assortative mating by age, and from an evolutionary point of view is a bet-hedger, laying commonly 3 eggs, although few chicks survive, except in years with high food availability, when the strength of densitydependence is low [START_REF] Oro | Audouin's gull account[END_REF] .

Individual data

During 1988-2017 a total of 30,290 individuals were captured and ringed as chicks at the Punta de la Banya [START_REF] Genovart | Exploiting uncertain ecological fieldwork data with multi-event capture-recapture modelling: an example with bird sex assignment[END_REF][START_REF] Oro | Comparing demographic parameters for philopatric and immigrant individuals in a long-lived bird adapted to unstable habitats[END_REF] . From 2002 to 2017, resightings were made using spotting scopes from a distance all over the western Mediterranean with a total of 63,106 resights in the study area and 5,593 different individuals resighted. Each year we recorded the breeding patch for each individual. To make sure that individuals were breeding and that they did so in a particular patch, we only selected those individuals seen during the breeding season in a particular patch showing unequivocal breeding behaviour. Specifically, individuals making alarm calls, incubating eggs or with chicks.

After this selective filter, our final database included 3,548 individuals.

SNA framework

Our social network, defined as the observed pattern of breeding association, was constructed taking individuals (N=3548) as the nodes of the network and each edge dyad (i.e. pair of individuals) representing the fact that individuals breed in the same patch. We ended with a global sociomatrix, i.e. the matrix representation of the dyadic relationships among individuals, of 3,548* 3,548. Edges showed if two individuals bred in the same patch at least once in a certain period (see below). The network was not directional. Based on previous results on population dynamics, and on the population size of this species and colony [START_REF] Genovart | Immature survival, fertility, and density dependence drive global population dynamics in a long-lived species[END_REF][START_REF] Payo-Payo | Colonisation in social species: the importance of breeding experience for dispersal in overcoming information barriers[END_REF] , we divided our dataset in two main periods: a period defined as "stable phase" from 2002 to 2011, and a period of "transition phase to collapse", from 2012 to 2017 (Figure 1).

We used the recently developed AME function from the AMEN package [START_REF] Hoff | Additive and multiplicative effects network models[END_REF][START_REF] Minhas | Inferential Approaches for Network Analyses: AMEN for Latent Factor Models[END_REF] , that can be applied to binary, ordinal, and continuous network data. This new approach is a random-effects regression model; uses an iterative Markov chain Monte Carlo (MCMC)

algorithm that provides Bayesian inference of the parameters in the social relations regression model (SRM; [START_REF] Warner | A new round robin analysis of variance for social interaction data[END_REF] ) using additive and multiplicative effects and combining the linear regression model with the covariance structure of the SRM [START_REF] Minhas | Inferential Approaches for Network Analyses: AMEN for Latent Factor Models[END_REF] . The AME method, not currently used in research on animal social networks is also able to cope with missing and censored data, our data set complying with the assumption that individuals are missing at random. Coping with missing data is highly relevant when analysing sociality on wild populations, as detection rate for individuals is almost always imperfect, and properly controlling for missed observations is a very important step in social network analysis [START_REF] Gimenez | Inferring animal social networks with imperfect detection[END_REF][START_REF] Hoppitt | Association indices for quantifying social relationships: how to deal with missing observations of individuals or groups[END_REF] . To create and visualize our networks we used the packages Amen 68 , Asnipe [START_REF] Farine | Animal social network inference and permutations for ecologists in R using asnipe[END_REF] , gdata [START_REF] Warnes | gdata: various R programming tools for data manipulation[END_REF] and igraph in R [START_REF] Csardi | The igraph software package for complex network research[END_REF] .

Are there social ties that persist over time?

We investigated if individuals create social ties that persist over time longer than one breeding occasion by means of two approaches: i) contingency tables and ii) the inclusion of time dependent regression terms in the AME modelling framework [START_REF] Hoff | Additive and multiplicative effects network models[END_REF] (see previous section). We used both methods because this is the first application of the AME approach in an ecological context. We analysed data of the period of stability, from 2002 to 2011, dividing this period into two sub-periods of five years (2002-2006 and 2007-2011). In the contingency table approach, we tested if the probability of breeding together at least once during the second sub-period was independent of having bred together at least once in the first period. We built a 3x3 table of frequencies,

showing the frequencies of two individuals breeding or not together at least once during the second sub-period depending on whether they bred together or not at least once during the first sub-period, and pulling apart those dyads with missing data. We then tested for deviation of random frequencies by Chi Square test.

With the SNA approach, we analysed the social ties between individuals using the AME function provided in the Amen package in R and including data from the first period (five previous years) as predictors of association during the second sub-period. We considered that this time window was not too large to include important death events, but large enough to account for the imperfect detection of individuals. To achieve convergence, we increased the number of iterations to 100,000 from the default value of 10,000 and lengthened the burn-in period to 500.

Had perturbations affected social cohesion in this species?

To assess if perturbations affected social structure in this species, we analysed as previously, with both the contingency table approach and the SNA approach, the social ties during the period of "transition phase to collapse" (2012-2017). To do so, we tested if the probability of breeding together in this phase (2012-2017) was independent of having bred together during five previous years (2007-2011). We then compared these results from those previously observed during the "stability phase".

A potential concern was the reduced power during the collapse period because the number of individuals decreased from the stability period. In order to have a similar power in both analyses, we performed the contingency table analysis during the stability period by drawing at random a number of observed associations equal to the number of observed associations during the collapse. We did this repeatedly (1000 times) and calculated which percentage of times the resulting chi-square was significant at the 5% level.

Regarding the SNA approach, it is advised [START_REF] Farine | A guide to null models for animal social network analysis[END_REF] to do permutations on the raw data prior to the analysis and compare the result of some relevant statistic obtained with the original data to the distribution of the same statistic over the permutations. We chose the regression coefficient of the association of a dyad on the previous association of the same dyad as our statistic of social cohesion. This statistic is provided by the function ame of the package amen [START_REF] Hoff | Additive and multiplicative effects network models[END_REF] . We calculated this statistic on the original data. Then, within each year, we rearranged randomly the individuals among the patches, keeping the same number of individuals within each patch. We did this 200 times and calculated each time the regression coefficient. Then, we situated the value of the regression coefficient from the original data among the distribution of regression coefficients from the permutated data and examined how extreme it was. If it was among the 5% extreme values, we concluded that there was indeed a non-random association of individuals within the patches.

Figure 1 .

 1 Figure 1. Number of breeding pairs in the Punta de la Banya colony from colonization

Figure 2 .

 2 Graphical representation of social networks by the association between individuals of Audouin's gulls in breeding patches comparing a) the stability phase (2002-2011) and b) the transition phase to colony collapse (2012-2017) (see Figure 1). Each node represents an individual and each edge links those individuals that have bred together in the same patch. We used the halfweight association index (HWI) to estimate the strength of relationship between pairs of individuals an index more suitable when not all individuals within each

Figure 3 .

 3 Figure 3. Map of the study area comprising the 3 main colonies and the distribution of patches within colonies during the study period at a) Sant Antoni and b) Punta de la Banya. Sant Carles de la Rápita colony is considered to have only one patch.

Table 1 .

 1 Results of the AME regression function to test if there were social ties between

		pmean psd	z-stat	p-value pmean psd	z-stat	p-value
	intercept -1.467 0.087 -16.808 0.000	-0.231 0.018 -12.629 0.000
	.dyad	6.420	0.328	19.576 0.000	-0.004 0.022 -0.188	0.851

individuals while breeding during the stable period and during the transition phase to collapse period. The alternative hypothesis is that individuals aggregate annually at random for breeding. ".dyad": coefficient of the dependent regression term considering the previous dyadic relationship between individuals; "pmean": posterior mean estimate; "psd": posterior standard deviation; "z-stat ": nominal z-score.

Acknowledgements

We would like to thank all the people who have helped with the fieldwork in the Ebro delta over the years, particularly Julia Piccardo, Toni Curcó and the technical staff and volunteers at the Ebro Delta Natural Park. We would also like to thank Peter Hoff, for his advices and solving an analytical problem we encountered while using the AMEN R package. We also thank the Regional Government of Catalonia, for permits to access the study sites. Elisabeth Rochon corrected the English. Funding came from the Spanish Ministry of Science (CGL2017-85210), grant PICS INTERACT n°07699 (2016, CSIC-CNRS). MG was partially supported by the European Union (MINOUW Project, H2020-634495) and the Spanish Ministry of Science (CGL2017-85210). We have no conflict of interest to declare.

Authors' contributions MG conceived the idea; MG, OG, RP and RC designed methodology; MG, DO and AB collected the data; MG analysed the data; MG led the writing of the manuscript. All authors critically contributed to the drafts and gave final approval for publication.

Data accessibility

Data is available via CSIC repository.