Seung Gyu

Hyun Cheriton

Romain Lebreton
email: romain.lebreton@lirmm.fr

Éric Schost
email: eschost@uwaterloo.ca

Algorithms for structured linear systems solving and their implementation

ere exists a vast literature dedicated to algorithms for structured matrices, but relatively few descriptions of actual implementations and their practical performance. In this paper, we consider the problem of solving Cauchy-like systems, and its application to mosaic Toeplitz systems, in two contexts: rst in the unit cost model (which is a good model for computations over nite elds), then over Q. We introduce new variants of previous algorithms and describe an implementation of these techniques and its practical behavior. We pay a special a ention to particular cases such as the computation of algebraic approximants.

INTRODUCTION

In this paper, we discuss linear algebra algorithms inspired by the following kind of question: given polynomials such as

       t 0 = 8x 4 -8x 2 + 1 t 1 = 16x 5 -20x 3 + 5x t 2 = 32x 6 -48x 4 + 18x 2 -1,
nd that the relation t 0 -2xt 1 + t 2 = 0 holds (this is the recurrence between Chebyshev polynomials of the rst kind).

More precisely, given input polynomials (t 0 , . . . , t s-1 ) over a eld K, together with integers (n 0 , . . . , n s-1 ) and σ , the Hermite-Padé approximation problem asks to compute polynomials (p 0 , . . . , p s-1 ), not all zero, such that deg(p i ) < n i holds for all i, and such that we have p 0 t 0 + • • • + p s-1 t s-1 = O(x σ ). ere exist numerous applications to this type of question, very important particular cases being algebraic approximants (with t i = f i , for some given f ) or di erential approximants (with t i = d i f /dx i , for some given f ); see for instance [START_REF] Bostan | Algorithmes e caces en calcul formel[END_REF]Chapitre 7].

Expressed in the canonical monomial bases, the matrix of a Hermite-Padé problem has size σ × (n 0 + • • • + n s-1 ), and consists of s lower triangular Toeplitz blocks. More generally, our goal in this paper is to compute e ciently elements in the kernel of mosaic Toeplitz matrices [START_REF] Heinig | On the inverses of Hankel and Toeplitz mosaic matrices[END_REF], that is, m ×n matrices T = (T i, j ) 1≤i ≤p,1≤j ≤q with a p ×q block structure, each of the pq blocks T i, j being Toeplitz.

An m × n Toeplitz matrix T = (t i-j ) 1≤i ≤m,1≤j ≤n can be succinctly represented by the polynomial P T = t -n+1 + t -n+2 x + • • • + t m-1 x m+n-2 ; multiplication of T by a vector b = [b 0 • • • b n-1 ] t amounts to computing P T P b mod x m+n-1 and keeping the coecients of degrees n -1, . . . , m + n -2. More generally, a mosaic Toeplitz T = (T i, j ) 1≤i ≤p,1≤j ≤q can be described by a sequence of pq polynomials P = (P T i, j ) 1≤i ≤p,1≤j ≤q , together with the sequences I = (m 1 , . . . , m p ) and = (n 1 , . . . , n q ) giving the row-sizes and column-sizes of the blocks. en, our main problem is as follows.

P

A. Given P and integers I, as above, de ning a mosaic Toeplitz matrix T, nd a non-zero vector in the kernel of T.

We consider two situations, rst over an arbitrary eld K, counting all operations in K at unit cost, then over Q, taking bit-size into account. In all this paper, we closely follow the existing formalism of structured matrix computations developed in previous work by Morf [START_REF] Morf | Doubling algorithms for toeplitz and related equations[END_REF], Bitmead-Anderson [START_REF] Bitmead | Asymptotically fast solution of Toeplitz and related systems of linear equations[END_REF], Kailath and co-authors [START_REF] Kailath | Displacement ranks of matrices and linear equations[END_REF][START_REF] Kailath | Fast Reliable Algorithms for Matrices with Structure[END_REF], Pan [START_REF] Pan | On computations with dense structured matrices[END_REF][START_REF] Pan | Parametrization of Newton's iteration for computations with structured matrices and applications[END_REF], Kaltofen [START_REF] Kaltofen | Asymptotically fast solution of Toeplitz-like singular linear systems[END_REF], Cardinal [START_REF] Cardinal | On a property of Cauchy-like matrices[END_REF], etc.

We complement the existing literature as follows. First, we de ne a class of Cauchy-like matrices (see de nitions below) for which the matrix-vector product is faster by a constant factor than previous designs. Next, we show how Pan's technique of "multiplicative transformation" of operators results in matrices that have generic rank pro le (with high probability), so that further regularization is not needed in general. We then describe an improved iterative algorithm (of quadratic complexity with respect to the matrix size), that makes use of fast matrix multiplication; nally, for matrices de ned over Q, we introduce a divide-and-conquer algorithm as an alternative to Newton iteration, and we show how it can be improved in the case of algebraic approximation.

Another contribution of this paper is a discussion of the design and practical performance of a C++ implementation of these algorithms. To our knowledge, only a few papers address these methods from the practical viewpoint. An early reference is [START_REF] Sexton | Remark on a displacement-rank inversion method for Toeplitz systems[END_REF], which concludes that the divide-and-conquer (MBA) algorithm for solving Toeplitz matrices in quasi-linear time would require matrices of size 10 6 to break even with quadratic-time algorithms; a more recent article [START_REF] Huckle | Implementation of a superfast algorithm for symmetric positive de nite linear equations of displacement rank 2[END_REF] estimates the crossover point to be around 8000.

Our experiments rst consider computations over nite elds. We assess the practical impact of fast matrix multiplication for structured matrix algorithms (such as in the new algorithm mentioned above, or those in [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF][START_REF] Bostan | Solving structured linear systems with large displacement rank[END_REF]), and we estimate for what matrix size quasilinear algorithms become e ective (with much lower crossover points than above). In the context of computations over Q, we show that our divide-and-conquer algorithm outperforms Newton iteration consistently, and we demonstrate signi cant speed-ups for the case of algebraic approximants.

We start in Section 2 with a review of basic results on structured matrices, including a discussion of transformation of operators and regularization. Section 3 describes algorithms applicable in a unit cost model, with in particular a new algorithm that uses fast matrix multiplication, and a discussion of the practical performance of these algorithms. Finally, Section 4 presents li ing algorithms for structured matrices, and the corresponding implementation.

BASIC RESULTS

is section reviews background material on structured matrices; for a more comprehensive treatment, we refer the reader to [START_REF] Pan | Structured Matrices and Polynomials[END_REF].

Overview. Developed in [START_REF] Kailath | Displacement ranks of matrices and linear equations[END_REF], the displacement operator approach associates to a matrix A its displacement ∇(A), that is, the image of A under a displacement operator ∇. en, we say that A is structured with respect to ∇ if ∇(A) has a small rank compared to its size; the rank of ∇(A) is called the displacement rank of A with respect to ∇. A prominent example is the family of so-called Toeplitz-like matrices, which are structured for the Toeplitz displacement operator

ϕ : A → (Z m A -AZ n ) = (A ↓) -(A ←)
where the n × n lower shi matrix Z n is the matrix with ones below the diagonal. e displacement rank of a Toeplitz matrix for this operator is at most two; the displacement rank of a mosaic Toeplitz with a p × q block structure is at most p + q.

e key idea of most algorithms for structured matrices is summarized by Pan's mo o [START_REF] Pan | Structured Matrices and Polynomials[END_REF]: compress, operate, decompress. Indeed, for A of size m × n over a eld K, if ∇(A) has rank α, it can be represented using few elements through ∇-generators, that is, two matrices (G, H) in K m×α × K n×α , with ∇(A) = GH t ; α is the length of the generators. e main idea behind algorithms for structured matrices is to use generators as a compact data structure, involving α(m + n) eld elements instead of mn.

Cauchy-like matrices. Beyond the Toeplitz structure (and the directly related Hankel one), two other important cases are the socalled Vandermonde and Cauchy structures. While the case of Toeplitz-like matrices was the rst one to be studied in detail, we will actually focus on Cauchy-like matrices, as we will see that this particular structure is quite convenient to work with.

For a sequence u = (u 1 , . . . , u m ) in K m , let D u ∈ K m×m be the diagonal matrix with entries u 1 , . . . , u m . en, given u as above and v in K n , we will consider the operator ∇ u,v : A ∈ K m×n → D u A -AD v ; Cauchy-like matrices (with respect to the choice of u and v) are those matrices A for which ∇ u,v (A) has small rank.

Let u, v be given and suppose that u i j holds for all i, j. en, the operator ∇ u,v is invertible: given ∇ u,v -generators (G, H) of length α for A, we can reconstruct A as

A = α i=1 D g i C u,v D h i , C u,v = 1 u i -j 1≤i ≤m 1≤j ≤n , (1)
where g i and h i are the ith columns of respectively G and H, and matrix C u,v is known as a Cauchy matrix. Remark that we can equivalently rewrite A as

A = (GH t ) C u,v , (2)
where denotes the entrywise product. We will have to handle submatrices of A through their generators. e fact that D u and D v are diagonal matrices makes this easy (this is one of the aspects in which the Cauchy structure behaves more simply than the Toeplitz one). Suppose that (G, H) are generators for A, with respect to the operator ∇ u,v , and let u I = (u i ) i ∈I and v = ( j ) j ∈ be subsequences of respectively u and v, corresponding to entries of indices I and . Let A I, be the submatrix of A obtained by keeping rows and columns of indices respectively in I and , and let (G I , H ) be the matrices obtained from (G, H) by respectively keeping rows of G of indices in I , and rows of H of indices in . en, (G I , H ) is a ∇ u I ,v -generator for A I, .

Another useful property relates to inverses of Cauchy-like matrices. If a matrix A ∈ K n×n is invertible, and is structured with respect to an operator ∇ u,v , its inverse is structured with respect to ∇ v,u : if

D u A -AD v = GH t , one easily deduces that D v A -1 -A -1 D u = -(A -1 G)(A -t H) t ,
where A -t is a shorthand for (A -1 ) t .

Algorithms. In most of the paper, we give costs in an algebraic model, counting base eld operations at unit cost (in Section 4, we work over Q and use a boolean model).

We let M be such that over any ring, polynomials of degree at most d can be multiplied in M (d) base ring operations; we also assume that the super-linearity assumptions of [START_REF] Zur Gathen | Modern Computer Algebra[END_REF]Chapter 8] hold. Using the Cantor-Kaltofen algorithm [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF], we can take M (d) ∈ O(d log(d) log log(d)). We let ω be a feasible exponent for linear algebra, in the sense that matrices of size n can be multiplied in O(n ω ) base ring operations over any ring; the best bound to date is ω < 2.38 [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF][START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF]. We will have to compute rank and rank pro le of dense matrices; the cost reduces to that of matrix multiplication [START_REF] Ibarra | A generalization of the fast LUP matrix decomposition algorithm and applications[END_REF].

e notation O˜( ) indicates that we omit polylogarithmic terms. Matrix-vector multiplication with C u,v reduces to degree n polynomial interpolation at the points v and evaluation at the points u. Using fast polynomial evaluation and interpolation, this can be done in time O(M (ν ) log(ν )), with ν = max(m, n); thus, we can multiply matrix A of (1) by a vector in time O(αM (ν ) log(ν )) ⊂ O˜(αν). In [39, eorem 4.7.3], Pan shows that if the entries of both u and v are in geometric progression, one can reduce the cost of the matrix-vector multiplication by C u,v to O(M (ν)), since polynomial evaluation or interpolation at n points in geometric progression can be done in time O(M (n)) [START_REF] Bluestein | A linear ltering approach to the computation of the discrete Fourier transform[END_REF][START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF]; then, multiplication by A as above takes time O(αM (ν )).

Remark 2.1. We propose here a re nement of this idea, that allows us to save a constant factor in runtime: we require that u and v be geometric progressions with the same ratio τ . en, the Cauchy matrix C u,v has entries 1/(u ij ) = 1/(u 1 τ i-1 -1 τ j-1 ), so it can be factored as

C u,v = D τ -1 1 u 1 -1 τ j-i 1≤i ≤m 1≤j ≤n
, where D τ -1 is diagonal with entries (1, τ -1 , τ -2 , . . . , τ 1-m ), and where the right-hand matrix is Toeplitz. In the reconstruction formula (1), the diagonal matrix D τ -1 commutes with all matrices D g i , so we can take it out of the sum. Hence, we replaced α evaluations / interpolations at geometric progressions by α product by Toeplitz matrices, each of which can be done in a single polynomial multiplication. e cost for a matrix-vector product by A remains O(αM (ν )), but the constant in the big-O is lower: for m = n, using middle product techniques [START_REF] Bostan | Tellegen's principle into practice[END_REF][START_REF] Hanrot | e Middle Product Algorithm[END_REF], the cost goes down from 3αM (ν ) + O(αν ) to αM (ν ) + O(αν).

If one needs to multiply A by several vectors, further improvements are possible: we mention without giving details an algorithm from [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF], that itself follows [START_REF] Bostan | Solving structured linear systems with large displacement rank[END_REF], and which makes it possible to multiply A by α vectors in time O(α ω-1 M (ν )) instead of O(α 2 M (ν )), by reduction to a sequence of polynomial matrix multiplications.

Reduction from mosaic Toeplitz to Cauchy structure. Our primary interest lies in mosaic Toeplitz matrices. An important insight of Pan [START_REF] Pan | On computations with dense structured matrices[END_REF] shows that one can reduce questions about Toeplitz-like matrices to ones about Cauchy-like matrices (and conversely, if one wishes to), for a moderate cost overhead. In this paragraph, we give details on this transformation, following [START_REF] Pan | Structured Matrices and Polynomials[END_REF]Chapter 4.8].

To a vector u in K m , let us associate the corresponding Vandermonde matrices V u = u j-1 i 1≤i, j ≤m , W u = u m-i j 1≤i, j ≤m . For u as above, v in K n and T in K m×n , we let A = V u T W v ; we now prove that if T is Toeplitz-like, A is Cauchy-like.

e ∇ u,v -generators of A depend on the ϕ-generators of T, so let us start with those. If T = (T i, j ) 1≤i ≤p,1≤j ≤q is mosaic Toeplitz of size m × n, with T i, j of size m i × n j , we de ne the following.

For 1 ≤ i ≤ p, let T i, * be the block matrix

[T i,1 • • • T i,q ]
. Let η i be the rst row of T i, * , shi ed le once, let η i be the last row of T i, * , and nally let h i = (η i-1 -η i ) (for i = 1, η 0 is the zero vector); this is a row-vector of size n. If e i,n denotes the rowvector of dimension n with only a one at index i, we also let g i be the transpose of e m 1 +•••+m i -1 +1,m . For 1 ≤ j ≤ q, de ne T * , j similarly to T i, * ; λ j and λ j are its rst and last columns, with now a downward shi for λ j ; in addition, their entry of index m 1 +• • •+m j are set to zero. We de ne g p+j = (λ j -λ j+1 ) (λ q+1 is set to zero) and let h p+j = e n 1 +•••+n j ,n . With these de nitions, the matrices

G = [g 1 • • • g p+q ] and H = [h t 1 • • • h t p+q ] are ϕ-generators of T. To obtain A, we conjugate T by Vandermonde matrices, since D u V u -V u Z m and W v D v -Z n W v have small rank; explicitly, they are equal to respectively (u m ) t e m,m and (e 1,n ) t v n , where u m = [u m 1 • • • u m n ],
and similarly for v n . As a consequence, since ϕ(T) = GH t , we have

∇ u,v (A) = V u ϕ(T) W v + (u m ) t e m,m TW v -V u T(e 1,n ) t v n . (3)
is motivates the following de nitions. For 1 ≤ i ≤ p, let h i = h i W v and also let g i be the column of V u of index m 1 + • • • +m i-1 + 1.For 1 ≤ j ≤ q, de ne g p+j = V u g p+j and also let h p+j be the row of index n 1 + • • • + n j in W v . We de ne g p+q+1 = (u m ) t , and h p+q+1 as the last row of TW v ; we also de ne g p+q+2 as the rst column of -V u T, and

h p+q+2 = v n . By Eq. (3), G = [g 1 • • • g p+q+2 ] and H = [h 1 t • • • h p+q+2 t ] are ∇ u,v
-generators of A of length p +q + 2 (whereas T has ϕ-generators of length p +q). e bo leneck in the computation of G and H are p le -products by W v and q products by V u . If all entries of u and v are in geometric progression, it takes time O(pM (n) + qM (m)).

Regularization. In our algorithms for Cauchy-like matrices, we will assume that the input matrix has generic rank pro le, that is, that its leading principal minors of size up to its rank are invertible. Regularization for structured matrices was introduced for this purpose by Kaltofen [START_REF] Kaltofen | Asymptotically fast solution of Toeplitz-like singular linear systems[END_REF], for the Toeplitz structure. In our Cauchy context, one could apply to A as de ned above the regularization procedure from [39, Section 5.6], which consists in replacing A by A = D x C a,u AC v,b D y , for some new vectors x, a ∈ K m and y, b ∈ K n . eorem 5.6.2 in [START_REF] Pan | Structured Matrices and Polynomials[END_REF] shows that if a, u consist of 2m distinct scalars, and b, v consist of 2n distinct scalars, then there exists a non-zero polynomial ∆ in the entries of x and y, of degree at most µ = min(m, n) in each block of variables, such that the non-vanishing of ∆ implies that A has generic rank pro le (that theorem is stated for square matrices, but the result holds in the rectangular case as well). e downside of this construction is that it requires to compute a pair of generators of length p + q + 4 for A , involving the multiplication of G and H by C a,u and C v,b .

We now prove a useful property, involving only A = V u T W v as above. For the rest of this paragraph, we assume that the entries of u and v are indeterminates over K, and we show that A has generic rank pro le; this will imply the same property for a generic choice of u, v with entries in K. is construction is clearly favorable over the one above, since it involves no extra computation on A.

Following the proof of [39, eorem 5.6.2], we can use the Cauchy-Binet formula to express the minors of A in terms of those of V u , T and W v . Let i ≤ rank(T) and let I = {1, . . . , i}. e determinant δ i of the ith leading principal minor of A is the sum over all = {j 1 , . . . , j i } ⊂ {1, . . . , m} and K = {k 1 , . . . , k i } ⊂ {1, . . . , n} of α I, β , K γ K, I , where α I, is the determinant of (V u ) I, , β , K is the determinant of T , K , and γ K, I is the determinant of (W v ) K, I .

Let < be the monomial ordering on the variables u 1 , . . . , u i , 1 , . . . , i that rst sorts monomials using the lexicographic order on u 1 , . . . , u i with u 1 < • • • < u i , and then breaks the ties using the lexicographic order on 1 , . . .

, i with 1 > • • • > i . If = {j 1 , . . . , j i } with j 1 < • • • < j i then the leading monomial (denoted lm(•)) of α I, is u := u j 1 -1 1 • • • u j i -1 i . Similarly if K = {k 1 , . . . , k i } with k 1 < • • • < k i , then v K := n-k 1 1 • • • n-k i i = lm(γ K, I
). Since T is of rank greater of equal to i, at least one of its ith minor β ,K does not vanish. Let max , K max be the pair of subsets that maximizes u v K among those for which β , K 0. en we must have lm(det(A I, I )) = u max v K max , which shows that det(A I, I ) is non-zero.

e partial degree of det(A I, I ) in any variable is at most max(n, m).

OVER AN ABSTRACT FIELD

In this section, we work over a eld K, and we explain how to solve Problem A by using Pan's reduction to the Cauchy-like Problem B below. Even though our goal is to solve a linear system, the algorithms do slightly more: they compute the inverse of a given matrix (or of a maximal minor thereof); this is similar to what happens for dense matrices, where it is not known how to solve linear systems with an exponent be er than ω. en, the main question of this section is the following (the assumptions on the vectors u, v are slightly stronger than the one required for ∇ u,v to be invertible).

P B. Consider u = (u 1 , . . . , u m ) and v = ( 1 , . . . , n ), such that (u, v) has m + n distinct entries. Given ∇ u,v -generators of length α for A in K m×n , with α ≤ min(m, n), do the following.
If A does not have generic rank pro le, raise an error; else, return ∇ v ,u -generators for the inverse of the leading principal minor A r of A of order r , with v = ( 1 , . . . , r ), u = (u 1 , . . . , u r ) and r = rank(A).

To solve an instance of Problem A, with input matrix T, we apply the transformation of the previous section, to obtain a Cauchy-like matrix A = V u TW v . We compute generators of A -1 r , where A r is the maximal leading minor of A, and we return a vector b with

b = W v -A -1 r Bc c
, where A r, * = A r B and c ∈ K n-r random. ere exist two major classes of algorithms for handling Problem B, iterative ones, of cost Θ(mn) (for xed α), and divide-andconquer algorithms, of quasi-linear cost in m + n. We stress that having a fast quadratic-time algorithm is actually crucial in practice: as is the case for the Half-GCD, fast linear algebra algorithms, etc, the divide-and-conquer algorithm will fall back on the iterative one for input sizes under a certain threshold, and the performance of the la er will be an important factor in the overall runtime.

Iterative algorithms that solve a size n Toeplitz system in time O(n 2 ) have been known for decades [START_REF] Durbin | e ing of time series models[END_REF][START_REF] Levinson | e Wiener RMS error criterion in lter design and prediction[END_REF][START_REF] Trench | An algorithm for the inversion of nite Toeplitz matrices[END_REF]; extensions to structured matrices were later given, as for instance in [START_REF] Kailath | Fast Gaussian elimination with partial pivoting for matrices with displacement structure[END_REF]. A er a section of preliminary results, we give in Section 3.2 an algorithm inspired by [START_REF] Mouilleron | Algorithmes rapides pour la résolution de problèmes algébriques structurés[END_REF]Algorithme 4], for the speci c form of our Problem B. In this reference, Mouilleron sketches an algorithm that solves Problem B in time O(αn 2 ), in the case where m = n and A is invertible (but without the rank pro le assumption); he credits the origin of this algorithm to Kailath [29, §1.10], who dealt with symmetric matrices. Our algorithm follows the same pa ern, but reduces the cost to O(α ω-2 mn).

In Section 3.3, we review divide-and-conquer techniques. Kaltofen [START_REF] Kaltofen | Asymptotically fast solution of Toeplitz-like singular linear systems[END_REF] gave a divide-and-conquer algorithm that solves the analogue of Problem B for Toeplitz-like matrices, li ing assumptions of (strong) non-singularity needed in the original Morf and Bitmead-Anderson algorithm [START_REF] Bitmead | Asymptotically fast solution of Toeplitz and related systems of linear equations[END_REF][START_REF] Morf | Doubling algorithms for toeplitz and related equations[END_REF]; a generalization to the Cauchy case is in [START_REF] Cardinal | On a property of Cauchy-like matrices[END_REF][START_REF] Pan | Superfast algorithms for Cauchy-like matrix computations and extensions[END_REF]. A further improvement due to Jeannerod and Mouilleron [START_REF] Jeannerod | Computing speci ed generators of structured matrix inverses[END_REF], following [START_REF] Cardinal | On a property of Cauchy-like matrices[END_REF], allows one to bypass costly compression stages that are needed in Kaltofen's algorithm and its extensions, by predicting the shape of the generators we have to compute. For the case of square Cauchy-like matrices of size n, this results in an algorithm of cost O(α 2 M (n) log(n) 2 ), but we will point out that be er estimates are available by choosing u, v suitably.

Cauchy generators of the inverse

Let (G, H) ∈ K m×α × K n×α be ∇ u,v -generators of a matrix A, with u = (u 1 , . . . , u m ) and v = ( 1 , . . . , n ).
Let further r be the rank of A. Our goal is to decide if A has generic rank pro le, and if so, to return generators (Y, Z) ∈ K r ×α × K r ×α of the inverse of the leading principal minor of A. Below, we write µ = min(m, n).

For 0 ≤ i ≤ µ, write A = A (i ) 0, 0 A (i ) 0, 1 A (i ) 1, 0 A (i ) 1, 1 with A (i) 0,0 ∈ K i×i . If A (i) 0,0 is invertible, set S (i) = S (i) 0,0 S (i) 0,1 S (i) 1,0 S (i) 1,1 =       A (i) 1,1 -A (i) 1,0 A (i) 0,0 -1 A (i) 0,1 A (i) 1,0 A (i) 0,0 -1 -A (i) 0,0 -1 A (i) 0,1 A (i) 0,0 -1       (
this is similar to [START_REF] Cardinal | On a property of Cauchy-like matrices[END_REF], where the order of block rows and columns was however di erent). e sequence of matrices (S (i) ) i=0...µ starts from A = S (0) and ends at A -1 = S (µ) , at least when A is square and invertible. To understand where these matrices come from, we introduce the matrix R (0) = A I µ,n in K (m+µ)×n , where I µ,n ∈ K µ×n is the rectangular matrix with ones on the main diagonal. If one applies to R (0) the column operations that reduce its rst ith rows

A (i ) 0, 0 A (i ) 0, 1 to [ 0 I i,i ] using A (i)
0,0 as pivot, we get

R (i) = S (i) I µ,n = P i R (0) -A (i) 0,0 -1 A (i) 0,1 A (i) 0,0 -1 I n-i,n-i 0 , (4)
P being a cyclic permutation matrix than moves the rows up once. Given integers a, b, u a:b denotes the sequence (u a , . . . , u b ), and similarly for v a:b ; we then de ne u (i) = (u i+1:m , v 1:i ) and v (i) = (v i+1:n , u 1:i ). A key result for the sequel is Lemma 3.1, which is [14, Proposition 1]; it gives ∇ u (i ) ,v (i ) -generators of S (i) (remark that this operator is invertible, in view of our assumption on u and v).

For i as above, decompose G, H as G = G (i ) 0 G (i ) 1 , H = H (i ) 0 H (i ) 1 with G (i) 0 , H (i) 0 ∈ K i×α and de ne (Y (0) , Z (0) ) = (G, H) and
Y (i) = -A (i) 1,0 A (i) 0,0 -1 G (i) 0 + G (i) 1 -A (i) 0,0 -1 G (i) 0 , Z (i) = -A (i) 0,1 t A (i) 0,0 -t H (i) 0 + H (i) 1 A (i) 0,0 -t H (i) 0 . L 3.1. (Y (i) , Z (i) ) are ∇ u (i ) ,v (i ) -generators for S (i) .
If A has generic rank pro le, this shows that for r = rank(A), (Y (r )

0 , Z (r ) 0 ) are ∇ v ,u -generators for A (r ) 0,0
-1 , for u , v as in Problem B, so they solve our problem; it remains to explain how to compute these matrices. Let i, j be non-negative integers with 0 ≤ i + j ≤ µ, and A (i) 0,0 invertible. Decompose S (i) , Y (i) and Z (i) into blocks

S (i) = S (i, j) 0,0 S (i, j) 0,1 S (i, j) 1,0 S (i, j) 1,1
,

Y (i) = Y (i, j) 0 Y (i, j) 1 , Z (i) = Z (i, j) 0 Z (i, j) 1 with S (i, j) 0,0 ∈ K j×j and Y (i, j) 0 , Z (i, j) 0 ∈ K j×α , write u (i) = [ u (i, j) 0 u (i, j) 1 ]
with u (i, j) 0 ∈ K j and the same for v (i) . Operations that derive S (i+j) from S (j) and S (i) from A = S (0) are similar; a direct calculation shows

S (i+j) =       S (i, j) 1,1 -S (i, j) 1,0 S (i, j) 0,0 -1 S (i, j) 0,1 S (i, j) 1,0 S (i, j) 0,0 -1 -S (i, j) 0,0 -1 S (i, j) 0,1 S (i, j) 0,0 -1       . ( 5
)
is implies the following formulae on the generators:

Y (i+j,m-j) 1 = -S (i, j) 0,0 -1 Y (i, j) 0 , Z (i+j,n-j) 1 = S (i, j) 0,0 -t Z (i, j) 0 , Y (i+j,m-j) 0 = -S (i, j) 1,0 Y (i+j,m-j) 1 + Y (i, j) 1 , (6)
Z (i+j,n-j) 0 = -S (i, j) 0,1 t Z (i+j,n-j) 1 + Z (i, j)
1 ; this generalizes previous formulae for (Y (i) , Z (i) ). e idea behind Eq. ( 6) is that by decomposing the column operations that reduce [ A (i +j) 0, 0 A (i +j) 0, 1 ] to [ 0 I i +j, i +j ] (see Eq.( 4)) into two operations that respectively reduce the rst ith rows using A (i) 0,0 as pivot, then the next j rows using S (i, j) 0,0 as pivot, we split the computation of R (i+j) into two similar operations that map R (0) to R (i) and R (i) to R (i+j) .

Finally, to solve Problem B, we need to detect when A has generic rank pro le. By construction, S (i, j) 0,0 is the Schur complement of A (i) 0,0 , seen as a submatrix of A (i+j) 0,0 . is implies the following result.

A faster iterative algorithm

Using the previous discussion, we describe in Algorithm 1 a quadratictime iterative algorithm for Problem B; we will use a step size β ∈ {1, . . . , α }, given as a parameter. e formulae in (6) that compute (Y (i+j) , Z (i+j) ) from (Y (i) , Z (i) ) only involve S (i, j) 0,0

-1 , S (i, j) 0,1 , S (i, j) 1,0 . Hence, line 8 recovers S (i, j) 0,1 from its ∇ u (i, j) 0 ,v (i, j) 1 -generators (Y (i, j) 0 , Z (i, j)
1 ) and does the same for S (i, j) 1,0 (Lemma 3.1); this allows us to apply [START_REF] Bostan | Algorithmes e caces en calcul formel[END_REF].

For the correction of the algorithm, note that A has generic rank pro le i the rank of A is r max = max{r | ∀ ≤ r , det(A ( ) 0,0 ) 0}.

1 i = 0, Y (0) = G, Z (0) = H, r = µ 2 while i r do 3 j = min(β, r -i) 4
Recover S (i, j) 0,0 from its generators (Y (i, j)

0 , Z (i, j) 0 ) 5
Compute the rank ρ, rank pro le, inverse of S (i, j) 0,0 6 if S (i, j) 0,0 does not have generic rank pro le then raise error

7 if ρ < j then j = ρ, r = i + ρ 8 Recover S (i, j) 0,1 , S (i, j)
1,0 from their generators 9

Compute the generators (Y (i+j) , Z (i+j) ) of S (i+j)

10 i = i + j
11 Recover S (r ) from its generators (Y (r ) , Z (r ) )

12 Read in S (r ) the Schur complement of A (r ) 0,0 in A 13 if the Schur complement is non zero then raise error 14 else return (r ,

Y (r,m-r ) 1 , Z (r,n-r ) 1
) e while loop of the algorithm will nd the maximal leading principal minor of A that is non-zero and has generic rank pro le. If at some step i, we know that A (i) 0,0 is invertible and has generic rank pro le, we use Lemma 3.2 to test if A (i+j) 0,0 has the same properties. e loop exits on one of two conditions that both ensures that r = r max . Either rank(S (i, j) 0,0 ) is always j, or rank(S (i, j) 0,0 ) < j at some point and a er se ing r = i + rank(S (i, j) 0,0 ), we have det(A ( ) 0,0 ) 0 for all ≤ r and det(A (r +1) 0,0 ) = 0; we exit the while loop right a er. A er the while loop, it remains to test if rank(A) = r , which is done by checking that the Schur complement of A (r ) 0,0 in A is zero. Computing the rank pro le at line 5 costs O(β ω ). Using block matrix multiplication in Eq. ( 2), we recover the matrices of line 8 in time O(β ω-2 α(m + n)), since j ≤ β ≤ α. To compute Y (i+j) at line 9, we rst compute S (i, j) 0,0

-1 in time O(β ω ), S (i, j) 0,0 -1 Y (i, j) 1 in time O(β ω-1 α)
; then all other operations take time O(β ω-2 αm). Similarly, computing Z (i+j) takes time O(β ω-2 αn). We recover S (r ) at line 11 from its ∇ u ( r ) ,v ( r ) -generators (Y (r ) , Z (r ) ) by means of (2), with a cost O(α ω-2 mn) using block matrix multiplication.

One iteration of the while loop costs O(β ω-2 α(m + n)); we iterate O(µ/β) times, for a total cost of O(β ω-3 α µ(m + n)), which is O(β ω-3 αmn). is dominates the cost of line 11, so that the whole running time is O(β ω-3 αmn). e algorithm of [START_REF] Mouilleron | Algorithmes rapides pour la résolution de problèmes algébriques structurés[END_REF] uses β = 1, for which the cost is O(αmn); choosing β = α, we bene t from fast matrix multiplication, as the cost drops to O(α ω-2 mn).

e divide-and-conquer algorithm

We nally review the divide-and-conquer approach to solving Problem B, as a basis for the discussion of the next subsection. Algorithm 2 follows [START_REF] Jeannerod | Computing speci ed generators of structured matrix inverses[END_REF], recast in our framework, with the minor di erence that do not assume A invertible, and that we explicitly check if A satis es the generic rank pro le assumption. Recursive calls on lines 3, 8 will raise an error if the respective input matrix A (i) 0,0 and its Schur complement

A (i) 1,1 -A (i) 1,0 A (i) 0,0 -1 A (i) 0,1
are not in generic position, so we assume that we are not in that Algorithm 2: Divide-And-Conquer algorithm DAC Input : Generators (G, H) of A, threshold ν 0 Output :r = rank(A), generators of A linear algebra. us, we focused on comparisons between our own implementations of the various techniques seen so far, and did not a empt comparisons with other systems. e main subroutines we need are polynomial and matrix computations. NTL already o ers e cient FFT-based polynomial arithmetic; matrix multiplication over small elds Z/pZ, for p < 2 23 , is now extremely e cient, comparable to reference implementations such as FFLAS-FFPACK [START_REF]FFLAS-FFPACK: Finite Field Linear Algebra Subroutines / Package[END_REF]. On top of this, we implemented a fast polynomial matrix multiplication, using the cyclotomic TFT of [START_REF] Arnold | A Truncated Fourier Transform middle product[END_REF].

We rst discuss Problem B. Our rst tests compare our new O(α ω-2 n 2 ) algorithm to that of [START_REF] Mouilleron | Algorithmes rapides pour la résolution de problèmes algébriques structurés[END_REF], with runtime O(αn 2 ). Except for very small values of n, say n < 50, for which the behavior uctuates rapidly, we found that the new algorithm becomes more e cient for rather small values of α: our crossover points are α = 8 or 9 for primes less than 2 23 , and α = 13 or 14 for larger word-size primes.

e following graph shows the time ratio between the algorithm of [START_REF] Mouilleron | Algorithmes rapides pour la résolution de problèmes algébriques structurés[END_REF] and our algorithm, for α = 30; the larger value of p used here is the FFT prime p = 82705526964617217 = 7 2 2 54 + 1. We see that for small primes, for which matrix multiplication is very e cient, our algorithm brings a substantial improvement. We consider next the divide-and-conquer algorithm. e key factor for the e ciency of this algorithm is the cost of multiplying an n × n Cauchy-like matrix of displacement rank α by α vectors. We compare the approach of cost O(α 2 M (n)) of Remark 2.1 to the algorithm of [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF], with cost O(α ω-1 M (n)), with a view of determining for what values of α (if any) the la er becomes useful.

For the former algorithm, because we are able to cache several FFTs, we found it slightly more advantageous to use NTL's FFT rather than TFTs. e runtime of the rst algorithm then displays the typical FFT staircase behavior, so that as n grows, the crossover value for α uctuates, roughly between 30 and 55. e following graph shows the time ratio between the algorithm of [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF] and the direct one, for α = 60; at best, the new algorithm wins by a factor of 2. e results do not depend much on the nature of the prime (since polynomial arithmetic is a signi cant part of the runtime, and behaves essentially in the same manner, independently of p). Using these results, we determined empirical crossover values ν 0 (α) to end the recursion in the divide-and-conquer algorithm and switch to the iterative algorithm. We expect ν 0 (α) to grow with α, since it behaves like the solution of αn 2 = α 2 M (n) log(n) (assuming here we do not use fast linear algebra, for simplicity). e following table reports some values for ν 0 (obtained by searching for ν 0 in increments of 100). e threshold is higher for small primes, since such primes bene t more from our new iterative algorithm. ese values being set, we show runtimes for solving Problem B modulo p = 65537 and p = 882705526964617217, for increasing values of n, with α = 5 and α = 50; we also show the runtime of a dense matrix inversion in the same size. For a small displacement rank such as α = 5, the runtime is essentially the same for these two primes; with α = 50, we observe a di erence, by a factor of up to 3. In any case, there is a clear gain over dense methods. We also determined, for a given value of n, the crossover value α 0 above which dense linear algebra becomes faster than structured methods. e value α 0 (n) grows quite regularly with n, good approximations being α 0 (n) 0.2n (for p < 2 23 ) and α 0 (n) 0.25n (for larger values of p). is means that there is a wide range of inputs for which structured methods can be of use.

Our solution of Problem A is a direct reduction to the Cauchylike case. Pu ing the problem into Cauchy form by means of the formulae of Eq. ( 3) accounts for a small fraction of the total runtime: between 5% and 10% for small α (say α < 10), and less than 5% for larger values of α, in all instances we considered.

In the thousands of experiments we made, for matrix sizes such as n ≥ 5000 and small primes such as p = 65537, we observed some instances where the Cauchy matrix did not have generic rank pro le. is never happened for p having more than 50 bits.

MODULAR TECHNIQUES

We now address Problem A in the particular case where K = Q. Several modular algorithms are available to solve dense linear systems over Q. A rst option relies on the Chinese Remainder eorem, solving the system modulo several primes p 1 , p 2 , . . . before reconstructing the solution, making sure to ensure consistency of the modular solutions when ker(T) has dimension more than 1. Other approaches such as Dixon's algorithm [START_REF] Dixon | Exact solution of linear equations using p-adic expansions[END_REF], Newton iteration or divide-and-conquer algorithms use one prime p, and li the solution modulo powers of p.

Newton iteration for structured matrices goes back to Pan's article [START_REF] Pan | Parametrization of Newton's iteration for computations with structured matrices and applications[END_REF] (for Toeplitz-like matrices), and is explained in detail in [START_REF] Pan | Structured Matrices and Polynomials[END_REF]Chapter 7]. To the best of our knowledge, the other approaches mentioned above have not been discussed in the literature on structured matrices, save for the second author's PhD thesis [START_REF] Lebreton | Contributions to relaxed algorithms and polynomial system solving[END_REF]. Our goal in this section is to rst brie y present some of these techniques and analyze their complexity for the problem at hand; we will then discuss their practical performance. We will highlight in particular the case of algebraic approximants, for which we are able to obtain signi cant improvements.

We denote by I : N → R a function such that integers of bit size at most d can be multiplied in I (d) bit operations; we can take [START_REF] Fürer | Faster Integer Multiplication[END_REF][START_REF] Schönhage | Schnelle Multiplikation großer Zahlen[END_REF]. As in [START_REF] Zur Gathen | Modern Computer Algebra[END_REF], we assume that d → I (d)/d is non-decreasing.

I (d) = O(d log(d) log log(d)) or I (d) = d log(d)2 O (log * (d ))

Solving square systems by li ing

In this subsection, we are given a prime power p t , where t is a power of 2, an n × n matrix A and a vector b, both with entries modulo p t , and we assume that A is invertible modulo p t , with inverse B. We discuss algorithms that solve the equation Ax = b, and we estimate their complexity when A is a structured matrix. To simplify the cost analysis, we assume here that p = O(1).

We take C A such that, for any t, we can compute Ax mod p t in O(C A I (t)) bit operations, given a vector x with entries de ned modulo p t . Below, we will assume that A is Cauchy-like as in Remark 2.1, and given by generators of length α, so that we can take C A ∈ O(αM (n)). We will see however that in the case of algebraic approximants, be er estimates are available.

We consider two approaches: a divide-and-conquer algorithm and Newton iteration, which both feature a running time linear in the target precision t. We start with the divide-and-conquer approach. A version of it is in [START_REF] Berthomieu | Relaxed p-adic Hensel li ing for algebraic systems[END_REF] for dense matrices; the PhD thesis [START_REF] Lebreton | Contributions to relaxed algorithms and polynomial system solving[END_REF] describes this algorithm for Toeplitz-like matrices.

Assume that we have obtained generators of length α for B, for instance by applying the algorithm of the previous section. us, at the leaves of the recursion, each product Bb mod p can be computed using O(αM (n)) bit operations; using our assumption on I , the total runtime is then O(C A I (t) log(t) + αM (n)t) bit operations. With our upper bounds on C A , this simpli es further as O(αM (n)I (t) log(t)).

We turn next to Newton iteration. e matrix form of Newton iteration computes the inverses B k = A -1 mod p 2 k by B k +1 = Algorithm 3: Divide-And-Conquer algorithm DAC Q Input : A, b, B mod p, p, t as above Output : a solution of Ax = b mod p t 1 if t = 1 then return B b mod p 2 Compute x 0 =DAC Q (A, b, B, p, t/2) 3 Compute r 0 = (Ax 0b) mod p t and r 1 = r 0 /p t /2 4 Compute x 1 =DAC Q (A, r 1 , B, p, t/2) 5 return x 0 -p t /2 x 1 mod p t 2B k -B k AB k mod p 2 k +1 ; once they are known, we deduce the solution to our system by means of a matrix-vector product.

Let (G, H) be generators for A. Pan's insight was to use the relation above to compute the generators (X, Y) = (-BG, B t H) for

B. Given (X k , Y k ) = (X, Y) mod p 2 k , we can reconstruct B k , but to deduce (X k +1 , Y k +1 ), we have to multiply G and H by B k +1 or B t k +1
. is is done using the expression for B k +1 above, which gives

X k+1 = -(2B k -B k AB k ) G mod p 2 k +1 Y k+1 = (2B k -B k AB k ) t H mod p 2 k +1 .
is time, we are multiplying A and B k (and their transposes) by matrices of size n × α, all computations being done modulo p 2 k +1 . Each multiplication by one vector takes O(αM (n)I (2 k )) operations (even if C A is less than O(αM (n)), multiplications by B k are the bo leneck). Since we multiply these matrices by α vectors, taking all steps into account, we arrive at a total of O(α 2 M (n)I (t)) bit operations. Using the algorithm of [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF], this can be further reduced to O(α ω-1 M (n)I (t)); this improvement was not implemented.

Altogether, because it computes (generators of) a whole inverse, Newton iteration is slower by a factor α (or slightly less, if we use [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF]); on the other hand, it saves a factor of log(t). is is similar to what one observed when comparing these techniques for e.g. the solution of di erential equations [START_REF] Bostan | Power series solutions of singular (q)-di erential equations[END_REF][START_REF] Bostan | Fast computation of power series solutions of systems of di erential equations[END_REF].

Solving Problem A

e techniques above allow us to solve instances of Problem A, that is, nd nullspace elements for a mosaic Toeplitz matrix T, but they would require the knowledge of a maximal invertible minor of it. However, there is no guarantee that T possesses such a minor that would be mosaic Toeplitz as well, and we do not know how to nd it e ciently. Hence, we rely on the transformation to the Cauchy structure / regularization technique of Section 2.

Over Q, this approach has a certain shortcoming: the output is

a vector b = W v -A (r ) 0, 0 -1 A (r ) 0, 1 c c , with A = V u TW v , A (r )
0,0 a maximal minor of it and c a random vector in K n-r , and due to the preconditioning, the entries of b are expected to be of large height (larger than what we may expect for a solution of T). When ker(T) has dimension 1, we are not a ected by this issue. Indeed, in this case, all solutions are of the form λb 0 , for some vector b 0 ∈ Z n whose bit-size can be bounded only in terms of T. Hence, it su ces to compute the solution b of the regularized system modulo a large enough integer N and normalize it by se ing one of its entries to 1; this gives us the normalization of b 0 mod N . is idea could be extended to the case of an arbitrary nullspace dimension, but it would be tantamount to li ing a whole nullspace basis. Instead, we will reduce the nullspace dimension by adding a new block of equations; when ker(T) has moderate dimension, this barely a ects the overall runtime (for particular applications to algebraic approximants, another solution is described below).

For simplicity, we discuss here a version of the algorithm that may run forever in unlucky cases. We choose a prime p and compute ∇ u,v -generators for A = V u TW v (mod p), for u, v as in Remark 2.1. We call the algorithm of the previous section, in order to determine the rank of A; if A does not have generic rank pro le, choose another u, v. A er an expected O(1) a empts, we obtain the rank of A (as a matrix over F p ), and thus the dimension s of its nullspace. If s is greater than one, we add a block of Toeplitz matrices having s -1 rows to T, with small random entries, and update A accordingly. Heuristically the new matrix T has nullspace dimension 1 (otherwise, add another block of equations).

Let d = rank(T), A d be the d × d top-le submatrix of A, and b d be the vector of the rst d entries of the last column of A. We assume that A has generic rank pro le, so that A d is invertible. We compute

x = A -1 d b d (mod p) and y = W v [x 1 , • • • , x d , - 1
] t (mod p), we normalize y by dividing it by its rst non-zero entry, and we set t = 1. While we either cannot apply rational reconstruction to the entries of y, or a er applying rational reconstruction to y we have Ty 0, we do the following: set t = 2t, update modulo p t the quantities A d , b d , x (use one of the algorithms of Section 4.1) and y, and divide y by its rst non-zero entry.

A complete analysis of this algorithm would quantify the primes of bad reduction, give bounds that allow us to stop li ing the solutions if we reduced modulo such a bad prime and study the previous reduction to nullity one of T; we leave this to future work. In any case, if the li ing stops, we have obtained a solution to our system. Due to the doubling nature of this procedure, the runtime is proportional to that of the algorithm for solving square systems used at line 3, plus the cost of rational reconstruction. e former is O(C A I (t) log(t) + αM (n)t) bit operations using divideand-conquer, and O(α 2 M (n)I (t)) or O(α ω-1 M (n)I (t)) using Newton iteration; the la er is O(nI (t) log(t)).

Finally, we mention the cost of Chinese Remaindering techniques: instead of computing the solution modulo the t-th power of a single prime p, we might want to solve the system modulo t primes of the same magnitude. If we assume that A remains invertible with generic rank pro le modulo all these primes, and all these primes are O(1), the runtime for solving the systems is now

O(α 2 M (n) log(n)t) or O(α ω-1 M (n) log(n)t)
, by the results of Section 3 (to this, we have to add the cost O(nI (t) log(t)) of Chinese Remaindering and rational reconstruction).

We conclude this subsection with an important particular case, the computation of algebraic approximants. We are given a power series f in Q[[x]], together with degree bounds d, e; our goal is to compute a polynomial P ∈ Q[x, ], with deg(P, x) ≤ d, deg(P, ) ≤ e, such that P(x, f ) = 0. For any σ ≥ 0, nding P with the above degree bounds and such that P(x, f ) = 0 mod x σ is a Hermite-Padé approximation problem of a very special kind (all input series are powers of f ). To our knowledge, no algorithm in the framework of Section 3 can exploit this extra structure; in the case of computations over Q, we will now see that improvements are possible.

First, we show how to simplify the reduction to nullity one. éorème 7.15 in [START_REF] Bostan | Algorithmes e caces en calcul formel[END_REF] shows that if such a P exists and is irreducible, then any Q ∈ Q[x, ] with the same degree bounds as above and such that Q(x, f ) = 0 mod x 2de+1 is a multiple of P. From this, one deduces easily that if we compute two such polynomials Q 1 , Q 2 , their GCD will generically be P itself. For all primes p except a nite number, the rank of our matrix T does not change modulo p, and the GCD of two random basis elements commutes with reduction modulo p (note that P mod p may not be irreducible anymore). Hence, we can (probabilistically) nd P mod p by computing two solutions Q 1 , Q 2 to the above Hermite-Padé problem modulo p and taking their GCD. is reveals the support of P; we can then re ne the degree bounds in our Hermite-Padé problem, which in turn reduces the nullity of matrix T to 1.

Next, we show how to speed-up algorithm DAC Q in this case. e block-Toeplitz matrix T of our Hermite-Padé problem has displacement rank α = O(e), with e + 1 Sylvester blocks having O(d) columns and O(de) rows; as a result, the naive estimate on the cost of the matrix-vector product by matrix A = V u TW v is C A = O(eM (de)). However, we can reduce this cost using babysteps / giant steps techniques: the bivariate modular composition algorithm of [START_REF] Nüsken | Fast multipoint evaluation of bivariate polynomials[END_REF] shows that we can do matrix-vector product by T using O(e 1/2 M (de) + e (ω+1)/2 M (d)) operations; since multiplications by V u and W v take O(M (de)), we obtain the improved estimate C A = O(e 1/2 M (de) + e (ω+1)/2 M (d)) in this case. Algorithm DAC Q is the only algorithm we know of that takes into account the extra structure of algebraic approximants; we expect that similar improvements are possible for di erential approximants, using the evaluation algorithm of [START_REF] Bostan | Fast algorithms for di erential equations in positive characteristic[END_REF].

Experimental Results

Our experiments are dedicated to solving instances of Hermite-Padé approximation, which is a useful particular case of Problem A. We discuss two families of problems, rst the approximation of general power series, then algebraic approximants. Timings are measured on the same machine as Subsection 3.4.

In both cases, we show two graphs: on the le , we have ve blocks, each with n columns; on the right, we have n blocks, each with n columns. With the notation of the introduction, this means we are looking for approximants (p 0 , . . . , p 4 ), resp. (p 0 , . . . , p n-1 ), with degree bounds (n, n, n, n, n), resp. (n, . . . , n). e displacement rank α of these matrices if 6, resp. n + 1.

We rst examine the results for general power series. Our experiments showed that the runtime grows predictably with respect to the input bit-size; as a result, we x the input coe cients to be 10 bit integers, so the number of li ing steps depends just on n. Matrices are generated so as to have 1 less row than they have columns; since the inputs are random, they have nullspace of dimension 1 (we also generated instances with nullity up to 10; using our heuristic to reduce nullity, runtimes were almost indistinguishable). e primes we use have 59 to 60 bits. For DAC Q and Newton iteration, the sharp increases indicate an additional li ing step.

Newton iteration is slower than DAC Q , especially when the number of block grows (since its runtime is quadratic in the displacement rank α). Newton iteration should theoretically be competitive with DAC Q when the number of blocks is xed, and the size (and thus the output bit size) grow, since its theoretical runtime is be er by a log(t) factor, where t is essentially the output bit size. However, this is not noticeable on our experiments: in practice, the integer multiplication function I (d) grows like d 1+ε , for some ε > 0; in that case, the analysis of DAC Q can be re ned to O(αM (n)I (t)).

CRT, on the other hand, seems to be competitive with DAC Q when α is small but is signi cantly worse as α grows: the runtime CRT is not linear in α, while DAC Q is.

Next, we examine the computation of algebraic approximants; on the basis of the previous experiments, we consider algorithm DAC Q only. We start by generating a bivariate polynomial P(x, ) and compute one of its power series solutions f . Since we expect the coe cients of P to be smaller than the coe cients of f , we can be er control the behavior of the algorithms by choosing the bit size of P (here, it was xed to be 1000 bit integers). We compare algorithm DAC Q as in the general case, with C A = O(αM (n)), to the improved version using the bivariate modular composition (BMC) algorithm of [START_REF] Nüsken | Fast multipoint evaluation of bivariate polynomials[END_REF] described above, featuring a lower value for C A . e la er algorithm uses polynomial matrix multiplication, which we implemented using reduction modulo FFT primes and TFT polynomial multiplication.

As a result, we can see a signi cant di erence between the two algorithms as the size of the matrix grows: the theoretical speed-up predicted in Subsection 4.2 is observed in practice. Overall, the divide-and-conquer li ing algorithm turned out to be the most e cient method in all our experiments, especially as it can take extra structure into account, as in the case of algebraic approximants.

  is invertible and has generic rank pro le, A

Algorithm 1 :

 1 Iterative algorithm Iter for Problem B Input : Generators (G, H) of A and step size β Output :r = rank(A), generators of A

(r ) 0,0 -1

) and (Y (0) , Z (0) ) 5 if r 0 < i then 6 if Schur complement is non zero then raise error

) and (Y (r 0 ) , Z (r 0 ) ) case; then, the correction proof is similar to that in Section 3.2. Computations of lines 4, 9 are of the same essence; using Eq. ( 5) and ( 6), we recover the full generators (Y (i+j) ,

) of S (i, j) 0,0 -1 and previous full generators (Y (i) , Z (i) ) of S (i) . Rather than reconstructing the submatrices of S (i, j) as in the previous subsection, we use them through their generators: the formulae of ( 6) boil down to O(1) multiplications of Cauchy-like matrices (for which we have generators of length α) by α vectors.

In general, these multiplications cost O(α 2 M (ν ) log(ν )) operations, with ν = max(m, n). As pointed out in [39, eorem 5.3.1] (see also [START_REF] Chen | An e cient solution for Cauchy-like systems of equations[END_REF]), if u and v are geometric progressions, the cost for the matrix-vectors multiplication drops to O(α 2 M (ν )). If u and v have the same ratio (if this is the case at the top-level, this will remain the case for all recursive calls), by Remark 2.1, we can save a further constant factor. For large values of α, we can also apply the algorithm of [START_REF] Bostan | On matrices with displacement structure: generalized operators and faster algorithms[END_REF], which uses fast polynomial matrix multiplication to reduce the cost to O(α ω-1 M (ν )).

On line 6, the Schur complement is zero i Y Z = 0. is is tested by nding a minimal set of independent rows in Z (this takes time O(α ω-1 ν )) and multiplying their transposes by Y . ere are at most α such rows, so this takes time O(α ω-1 ν ) as well. Taking all recursive steps into account, the total cost of DAC is a factor log(ν ) times that of the Cauchy-like matrix products.

Experimental results

We implemented the algorithms described so far, together with all subroutines they rely on, in a C++ library available at h ps: //sea le.lirmm.fr/f/01b53dc420/. Our implementation is based on Shoup's NTL [START_REF] Shoup | NTL: A library for doing number theory[END_REF][START_REF] Shoup | A new polynomial factorization algorithm and its implementation[END_REF] version 10.3.0, and is dedicated to word-size primes (NTL's lzz p class); the divide-and-conquer algorithm of Subsection 3.3 actually requires FFT primes. In all that follows, timings are measured on an Intel i7-4790 CPU with 32 GB RAM; only one thread is used throughout. Some goals of our experiments are to assess whether fast matrix multiplication can bring practical improvements that re ect the theoretical ones, what are reasonable crossover points between iterative and divide-and-conquer algorithms, and what is the range of applicability of these structured methods compared to dense


Kevin Primicerio

Damien Challet

Stanislao Gualdi

Wisdom of the institutional crowd

come

Introduction

The collective ability of a crowd to accurately estimate an unknown quantity is known as the "Wisdom of the Crowd" [START_REF] Surowiecki | The wisdom of crowds[END_REF] (WoC thereafter). In many situations, the median estimate of a group of unrelated individuals is surprisingly close to the true value, sometimes significantly better than those of experts [START_REF] Galton | Vox populi (the wisdom of crowds)[END_REF][START_REF] Hill | Expert stock picker: the wisdom of (experts in) crowds[END_REF][START_REF] Landemore | Why the many are smarter than the few and why it matters[END_REF][START_REF] Nofer | Are crowds on the internet wiser than experts? the case of a stock prediction community[END_REF]. WoC may only hold under some conditions [START_REF] Surowiecki | The wisdom of crowds[END_REF][START_REF] Davis-Stober | When is a crowd wise?[END_REF]: for example social imitation is detrimental as herding may significantly bias the collective estimate [START_REF] Lorenz | How social influence can undermine the wisdom of crowd effect[END_REF][START_REF] Muchnik | Social influence bias: A randomized experiment[END_REF]. WoC is a reminiscent of collective rationality without explicit individual rationality: when it applies, it is a consistent aggregation of possibly inconsistent individual estimates [START_REF] Hogarth | A note on aggregating opinions[END_REF]. This is to be contrasted with the economic paradigm that collective rationality reflects individual rationality, where only a "typical" decision maker -the representative agent -is considered [START_REF] Hartley | The representative agent in macroeconomics[END_REF] or team reasoning where the individual agents explicitly optimize the collective welfare [START_REF] Colman | Collective rationality in interactive decisions: Evidence for team reasoning[END_REF].

Although almost all known examples of WoC are about a single number, there is no reason why it could not hold for whole functional relationships between several quantities. For example, Haerdle and Kirman analyze the prices and weights of many transactions in Marseille fish market: while the relationship between these two quantities is rather noisy, a local average does display the expected quantity discount [START_REF] Härdle | Nonclassical demand: A model-free examination of price-quantity relations in the Marseille fish market[END_REF]. More generically, many simple relationships found in Economic textbooks may only hold on average, but not for each agent or each transaction. Aggregation of quite diverse individual actions is still an open problem [START_REF] Kirman | Whom or what does the representative individual represent?[END_REF].

Asset price efficiency is an obvious instance of WoC in Finance: it states that current prices, determined by the actions of many traders, are the best possible estimates and fully reflect all available information [START_REF] Malkiel | Efficient capital markets: A review of theory and empirical work[END_REF][START_REF] Malkiel | The efficient market hypothesis and its critics[END_REF][START_REF] Fama | Market efficiency, long-term returns, and behavioral finance[END_REF]. Another WoC candidate in Finance is portfolios. The rationale is that many market participants, especially investment funds, strive to build optimal portfolios following their own criteria and constraints. If anything, the average portfolio structure of this population may be close to that of a rational benchmark.

In this paper we follow an approach similar to [START_REF] De Lachapelle | Turnover, account value and diversification of real traders: evidence of collective portfolio optimizing behavior[END_REF], where the benchmark is a simple mean-variance portfolio optimization with a given transaction cost structure. The relationship between the optimal number of assets in a portfolio and the portfolio value depends strongly on the transaction cost structure. Such dependence is found in the average portfolio structure of individual and professional investors over many orders of magnitude of portfolio value. In this article we focus on much larger investment funds facing additional constraints and show how the latter shape WoC.

Results

The added value of an investment fund resides in how well the latter allocates capital to securities with respect to some predefined criterion (return on investment, risk, performance tracking, etc.). Many kinds of funds exist, each with its own objective and methods. Thus, understanding the minute structure of the portfolio of a given fund is both impractical and fortunately not necessary in our approach as we do not need to assume individual rationality. Therefore, we shall focus on simple properties of portfolios whose average can be related to a rational benchmark. Generally speaking, at time t, Figure 1: Total mark-to-market value W i as a function of the number of investments n i , with a robust locally weighted regression fit (yellow line) and two linear fits (blue dashed lines) for two different ranges of n. Robust locally weighted regression fit for the simulated data (in green). fund i has capital W i (t) and may invest it into N (t) existing securities. As a result, each security α is found in m α (t) portfolios. The explicit time dependence is dropped hereafter for the sake of clarity. All the plots which are not time dependent refer to t = 2013Q1 but the same qualitative behavior is observed for other dates.

In the following we focus on the number of assets in portfolio of fund i, denoted by n i , and on its relationship with the total assets under management W i . We also investigate how the asset selection of funds shapes the relationship between the capitalization of asset α, denoted by C α , and the number of funds which have invested in this security, m α . The key point is that only n i depends on asset allocation strategies of fund i while W i , m α and C α are beyond its control. Given its scalar simplicity, n i must be regarded as a convenient and concise proxy of possibly complex asset allocation algorithms, thus as a useful simplification. WoC, if any, is to be found in the relation between n i and the three other quantities. Figure 1 plots W i versus n i in logarithmic scale: a cloud of point emerges, with a roughly increasing trend. The large amount of noise confirms the great diversity of fund allocation strategies. WoC may only appear in some average behavior. This is why we computed a locally weighted polynomial regression [START_REF] Cleveland | Local regression models[END_REF]. Remarkably, two distinct regions appear, each characterized by a roughly linear behavior. The axes of this figure being logarithmic, log W i = µ x log n i + i , where x ∈ {<, >} labels the region, and i is the deviation of W i from its local average W . In other words, defining n as the local average of n i , the average portfolio follows W ∝ n µx .

(

) 1
The boundary n * between the two regions is algorithmically determined for each quarterly snapshot (see S.I.), which allows us to measure the two exponents µ < and µ > . The latter are quite stable as a function of time (see Fig. 11 in S.I.); their time-averages µ < 2.1 ± 0.2 and µ > 0.3 ± 0.1 are markedly different, which points to distinct collective ways of building portfolios in these two regions. The cross-over point n * is also stable as time goes on (see Fig. 11 in S.I.). Let us start with the small diversification region (n i < n * ). First, the time average of exponent µ < is close to 2. The same exponent was found in the portfolios of individual and professional investors, whose portfolio values were much smaller than investment funds [START_REF] De Lachapelle | Turnover, account value and diversification of real traders: evidence of collective portfolio optimizing behavior[END_REF]. This is not entirely surprising: the value µ = 2 corresponds to optimal mean-variance portfolios with constant cost per transaction [START_REF] De Lachapelle | Turnover, account value and diversification of real traders: evidence of collective portfolio optimizing behavior[END_REF] (see also S.I.), provided that the same amount is invested in all the chosen assets (equally-weighted portfolios). One checks that large institutional portfolios are indeed very close to being equally weighted in this region. Take fund i; denoting by p iα = W iα /W i its investment fraction in security α, its diversification S i may be measured with Shannon Entropy S i =α p iα log p iα , which equals 1 and is maximal when all the non-null p iα are equal. As shown in Fig. 2, it is the case for low-diversification funds. The small deviations from S i = 1 in this region are due at least in part to the fact that the values of assets evolve after the initial investment and that portfolios are not rebalanced just before the dates of our dataset.

Since one finds the same exponent µ over many decades of portfolio values for a wide spectrum of market participants, and since µ = 2 corresponds to a realistic transaction cost per transaction, we argue that WoC is a plausible explanation of the average portfolio structure. Note that µ = 2 does not imply that funds really face constant transaction cost per transaction, only that their population acts as if it does. Finally, the fact that in this view WoC holds for a whole functional relationship, not only a single number, considerably extends its reach.

So far, bringing to light WoC only required to focus on the number of securities in a portfolio, not on how funds select securities. This implicitly assumed that funds could invest in all securities they wished, which is clearly not the case in the large diversification region: the fact that the exponent µ is much smaller in this region implies that funds need on average to split their investments into many more securities. This is most likely due to liquidity constraints: large funds cannot invest as much as they wish in some assets because there are simply not enough shares to build a position larger than a certain size without impacting too much their prices. Each fund has its own way to determine the maximal amount to invest in a given security α; a common criterion is to limit the fraction W iα /C α . Fig. 8 in S.I. strongly suggests that each fund fixes its upper bound

f (max) i ≥ max α f iα where f iα = W iα C α . (2)
It turns out that f max i is highly heterogeneous among funds log 10 (f max i ) -3.0 ± 1.0 (see Fig. 9). The existence of such limits implies that portfolios are less likely to be equally weighted in the large diversification region, which is confirmed by the decrease of the typical portfolio weights entropy (see Fig. 2).

Funds, however, do not invest in a randomly chosen security, even in the low diversification region. Figure 3 displays a scatter plot of the capitalization C α of each security α versus n α , the number of funds which have invested in this security, together with a local non-linear fit. As for W vs n, one finds a power-law relationship

log C α = γ log m α + α (3)
for large enough m (see S.I.). Hence C ∝ m γ , where m is the local average of n α . Exponent γ is stable during the period 2007-2014 (see Fig. 11 in S.I.) and its average γ 2.2 ± 0.1.

In short, one needs to introduce a model of how funds choose to invest in securities both to bring WoC to light in the large diversification region and to reproduce the average behavior of both Eqs (3) and 1. Since one sees a cross-over between two types of behavior rather than an abrupt change, we create logarithmic bins of the axis n i and denote the bin number of fund i by [n i ]. Two mechanisms must be specified: how a fund selects security α and how much it invests in it. The latter point is dictated by Fig. 8 in the large n i region where fund i invests W iα = f (max) i C α ; for the sake of simplicity, we approximate f [ni] . In the small diversification region, we assume that n i = n opt i , thus W iα = W i /n opt i to be consistent with our previous results. We choose a security selection mechanism that rests on the market capitalization C α of a security α (see S.I.) which is a good proxy of the liquidity (Fig. 10). We perform Monte-Carlo simulations from the empirical selection probabilities and f

(max) [ni]
and display the resulting W vs n and C vs m in Figs 1 and3 (continuous green lines), in good agreement with the local averages (continuous orange lines). One notices a discrepancy in the relationship C vs m for large n, which mainly comes from funds in the large diversification region. (See Fig 12 S.I).

The large diversification region illustrates how constraints may considerably modify the rational benchmark. While the above mechanism of security selection is able to reproduce adequately the behavior of well diversified funds, we could not find a rational benchmark for the dependence of f max and n i . Thus, the case for WoC in the large diversification region is not entirely closed. The main limitation of this dataset is that it provides accurate figures for long positions only. The other positions (short, bonds, ...) are most of the time only partially known. The frequency of the dataset is also inhomogeneous: data for most of the funds are quarterly updated (depending on regulations), hence we decided to restrict ourselves to 4 points in a year only. Such frequency is probably too low for investigating the dynamics of individual behavior but is not a problem for we focus on an aggregate and static representation of the investment structure.

Discussion and conclusion

While WoC is commonly applied to a population collectively guessing a single number, we investigate here a fundamentally different situation and provide evidence for a collective functional optimization of the asset ownership structure. What the reference function should be is dictated by optimality arguments. In the case of financial markets, the rational benchmark was not related to the efficient market hypothesis, but to the way a large population of professional fund managers build their portfolios. Whereas each fund has its own benchmark with respect to which the fund performance may be assessed, this, fortunately, has no discernible influence on the average structure of their portfolio. In addition, WoC is often meant as a collective guessing of non-experts; one thus may conclude that the population investigated here has decidedly more expertise than the subjects of other WoC studies. What kind of expertise the typical fund manager has is not obvious, at least when one looks at their pure performance (see e.g. [START_REF] Barras | False discoveries in mutual fund performance: Measuring luck in estimated alphas[END_REF]). In addition, the optimal relationship between the number of assets in a portfolio and the value of the latter is clearly not broadly known in these circles, as shown by the very large deviations from the ideal case in Fig. 1, and the collective expertise only appears when their decisions are suitably averaged. The presence of WoC when the subjects face strong constraints, as those of highly diversified funds, is more conjectural, and more work will be needed in that respect.

At a higher level, our results suggest that, while individuals may deviate much from the rational expectation theory, standard economic theory may hold at a collective level, without need for micro-founded individual decisions: the average decision may in some cases be approximated by a rational, representative agent. Our results however only hold on a snapshot of the system, for which individual fluctuations may be averaged out. In a dynamic setting, the very large deviations from the rational benchmark may not be neglected in the presence of feedback loops [START_REF] Gualdi | Tipping points in macroeconomic agent-based models[END_REF]. In other words, the dynamics of these fluctuations are worth investigating in their own right.

Supporting Information (SI) 2 Filtering

In order to remove inconsistencies in the dataset, we applied the following filters

Country of origin

Our dataset is sparse and heterogeneous. Indeed, the quality of the sources of data is directly related to each country's disclosure regulations. For these reasons we decided to keep only the entities which use an US based mail address. About 60% of the total market capitalization of the dataset is concentrated in US based securities. Figure 4 shows two large clouds of dots, each of them correspond to a different region of origin: green (resp. orange) cloud corresponds to non-US (resp. US) based securities. The origin of this large difference between these two regions are not clear: it could for example come from differences in regulations in non-US countries. It turns out that the ratio of the investment values in US and non-US assets varies little as a function of time (see Fig. 4), which does not affect the exponent µ in Eq. 1. As a consequence we focused on US securities.

Frequency

Large funds are requested to report their positions at a frequency which depends on the applicable regulation. As a result, reporting frequency ranges from monthly to yearly, most funds filing quarterly reports. We therefore focused of the latter.

Penny Stocks

The "penny stocks", i.e., usually securities which trade below $5 per share in the USA, are not listed on a national exchange. Since they are considered highly speculative investments and are subject to different regulations, we filtered them out.

Size

We also filtered out small founds and securities : W i > 10 5 USD, C α > 10 5 USD, n i ≥ 5, m α ≥ 10.

Output

We restricted our study to 36 quarterly snapshots starting from the first quarter of 2005 and ending with the last quarter of 2013. Figure 5 reports the evolution of the number of securities and funds in the database before and after filtering.

Asset selection

The framework we introduce in this paper follows a series of a few elementary steps described below. The aim is for the model to be sensitive to the different constraints which dominates the portfolio selection of a fund.

Finding n *

For date t, we define the cross-over point n * between the two regions which appear in the local polynomial regression. We determine this point value with a likelihood maximization of the model

W = µ < n + (µ > -µ < )(n -n * )θ(n -n * ), (4)
where θ(x) is the Heaviside function. We use a recursive method to find parameters µ < , µ > and n * [START_REF] Muggeo | Estimating regression models with unknown break-points[END_REF]. Figure 11 shows that n * is stable as a function of time.

Selection algorithm

Small diversification region n i < n *

In this region, we consider the equally weighted portfolio hypothesis to be true. Each position has a size Wi n opt i , where

n opt i
is the optimal number of position computed with eq 1. The funds select their asset randomly with a probability proportional to C α . Also, in order to construct an equally-weighted portfolio, a position is valid only if it is of size Wi n opt i .

Large diversification region n i ≥ n *

In this region, the liquidity constraints make it harder for funds to keep an equally weighted portfolio and portfolio values are thus spread on a larger number of assets. We model this process by introducing both a maximum fraction f which materializes the lack of liquidity and a probabilistic asset selection mechanism based on the rescaled rank of the security in terms of capitalization. In particular, we extract the rescaled rank from a beta distribution and take the corresponding security. The parameters of the distribution are fitted from the empirical density function (see later) Fig. 6. 1 1 We do not use the same rank-based selection mechanism in the low-diversification region because in this case it is harder to have a good fit with the beta distribution. This is however only a minor point since the capitalization is approximately power-law distributed and the two selection mechanisms are basically equivalent (the rank is proportional to a power of the capitalization). Results are indeed very similar in both cases. We chose a security selection mechanism which rests on the scaled rank of capitalization of security α, defined as ρα = rα N where rα is

Importance of capitalization ranking

Figure 6 shows that the distribution of investment of a fund is sensitive to its diversification n i . The Beta distribution, which is limited to a [0, 1] interval, is flexible enough to describe the asset selection mechanism of a fund.

f (x; a, b) = 1 B(a, b) x a-1 (1 -x) b-1 , (5)
where a and b are the shape parameters of the distribution, and B is a normalization constant.

Maximum investment ratio

The funds limit their investment in a given asset. They seem to follow a simple rule: defining the investment ratio f i,α = Wiα Cα , one easily sees in Fig. 8 that each fund has a maximum investment ratio

f max i = max α W iα C α (6)
Since the average exchanged dollar-volume of an asset is proportional to its capitalization (Fig. 10), the existence of f max i is a way to account for the available liquidity. Although that limit is clear for an individual fund, there is a large range of empirical values f max i Fig. 9.

Simulation

The simulation is done in a few simple steps:

1. Compute n * using the segmented model Eq. 4.

the rank of capitalization Cα at a given time. The selection probability P (W iα > 0|ρα) is then obtained by parametric fit to a beta distribution in each logarithmic bin. (b) Select assets randomly following a Beta probability distribution Fig. 6 with the parameters found in Fig. 7.

By iterating those steps we obtain Fig. 1 Since the simulation outputs a portfolio for every fund, we can directly infer the number of investors m α of every security. Figure 12: We separate the contribution from the low and highly diversified region. The origin of the discrepancy observed in Fig. 3 appears to be mainly due to the highly diversified region (Green dots for the empirical data, and orange dots for the model).

Figure 2 :

 2 Figure 2: Shannon Entropy S i as a function of the number of investments n i , each dot represents a single fund.

  (max) i by the median value of f (max) i in the bin [n i ], denoted by f (max)

Figure 3 :

 3 Figure 3: Market capitalization of securities as a function of the number of investors in logarithmic scale. From the local non-linear robust fit (yellow line) we observe a linear relationship for assets with more than about 100 investors. The blue dashed line corresponds to a linear fit on that group of asset. Hence W α ∝ m γ α , with γ 2.1. Robust locally weighted regression fit for the simulated data (in green).

Figure 4 :

 4 Figure 4: Top: Market capitalization as a function of the number of investors for all securities. Bottom: Temporal evolution of the aggregated market capitalization of US over the total market capitalization.

Figure 5 :

 5 Figure 5: Temporal evolution of the number of funds N i and securities N α in the data base. Unfiltered in dashed lines and US based only in solid lines.

Figure 6 :

 6 Figure 6: Top: Empirical probability density function of investing in a security of scaled rank ρ given the diversification n i of the fund. Bottom: Probability density function of investing in a security of scaled rank ρ given the diversification n i of the fund, given by the model.

Figure 7 :Figure 8 :

 78 Figure 7: Coefficients a and b of the Beta Distribution 5 as a function of n i

Figure 9 :

 9 Figure 9: Empirical probability density function of f max i

Figure 10 :

 10 Figure10: Market capitalization as a function of the daily exchange volume dollar. We find a slope close to 1 for all the dates in our database, confirming the hypothesis that the daily exchange volume dollar of an asset is proportional to its market capitalization.

Figure 11 :

 11 Figure 11: Top: Temporal evolution of the coefficients µ < , µ > and γ. Bottom: Temporal evolution of the value of the cross-over point n * between the two regions as a function of time. It reaches a global minimum after the 2008 crisis, and it also seems to decrease again after 2013.

Data

Our dataset consists of an aggregation of the following publicly available reports (in order of reliability): the SEC Form 13F, the SEC's EDGAR system forms N-Q and N-CSR and (occasionally) the form 485BPOS. Our work focuses on the period starting from the first quarter of 2005 to the last quarter of 2013.

These forms are filled manually and are thus error prone. We partially solve this issue by cross-checking different sources (which often contains overlapping information) and by filtering data before processing (see details in S.I.).