David Dureisseix

Generalized fraction-free LU factorization for singular systems with kernel extraction

Keywords: exact factorization, linear system, symbolic computations, integral domain, CFFLU, Gaussian elimination MSC codes: 65F05, 15A23, 68W30

Linear systems are usually solved with Gaussian elimination. Especially when multiple right hand sides are involved, an efficient procedure is to provide a factorization of the left hand side. When exact computations are required in an integral domain, complete fraction-free factorization and forward-backward substitutions are useful. This article deals with the case where the left hand side may be singular. In such a case, kernels are required to test a solvability condition and to derive the general form of the solutions. The complete fraction-free algorithms are therefore extended to deal with singular systems and to provide the kernels with exact computations on the same integral domain where the initial data take their entries.

Introduction

Linear systems with entries in an integral domain Z arise from several applications such as symbolic and exact polynomial computations [START_REF] Bini | Polynomial and matrix computations[END_REF][START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF], applications to cryptography [START_REF] Odlyzko | Discrete logarithms: The past and the future[END_REF], computational geometry [START_REF] O'rourke | Computational Geometry in C[END_REF][START_REF] Schirra | Robustness and precision issues in geometric computation[END_REF][START_REF] Halperin | Iterated snap rounding[END_REF], signal processing, etc.

Though other techniques are available, see [START_REF] Wan | An algorithm to solve integer linear systems exactly using numerical methods[END_REF] and associated references, we consider herein a direct solution technique. We denote such a problem with Ax = b, for which the right hand side b also takes its values in Z. If the left hand side is a square matrix, A ∈ Z n×n , and is regular, the Cramer's expression of the solution allows to conclude that it lies in the field of the quotient space Q, and that the solution of the scaled problem Ax = (det A)b, where the determinant of the left hand side is (det A) ∈ Z, also lies in Z n .

An efficient direct solving procedure, especially with multiple right hand sides, is to provide a factorization of the matrix A, relying on the Gaussian elimination technique. If the factorization could be performed without computing fraction, all the intermediate terms will also lie in Z and the solution of the scaled problem is exact; the initial solution is then (det A) -1 x ∈ Q n . Such a factorization is named as fraction-free [START_REF] Bareiss | Sylvester's identity and multistep integer-preserving Gaussian elimination[END_REF][START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF][START_REF] Turner | A simplified fraction-free integer Gauss elimination algorithm[END_REF][START_REF] Lee | Fraction-free gaussian elimination for sparse matrices[END_REF][START_REF] Zhou | Fraction-free forms of LU matrix factoring[END_REF][START_REF] Corless | The Turing factorization of a rectangular matrix[END_REF][START_REF] Mulders | A generalized Sylvester identity and fraction-free random Gaussian elimination[END_REF] and recently, a complete fraction-free algorithm has been proposed [START_REF] Turner | A simplified fraction-free integer Gauss elimination algorithm[END_REF][START_REF] Zhou | Fraction-free matrix factors: new forms for LU and QR factors[END_REF] for factorization, forward and backward substitutions.

Section 2 briefly recalls the original complete fraction-free algorithms, Section 3 deals with the case of singular matrices, providing a suited regularization, and the kernel computations, as well as the treatment of the surjection case with a rectangular matrix. Finally, Section 4 proposes two test examples.

2 Complete fraction-free LD -1 U factorization of a full rank matrix

The algorithms provided in [START_REF] Zhou | Fraction-free matrix factors: new forms for LU and QR factors[END_REF] are partially coded in in fflas library [START_REF] Dumas | Dense linear algebra over word-size prime fields: the FFLAS and FFPACK packages[END_REF] (BLAS for matrices over finite fields), and in the SymPy library (Python library for symbolic mathematics) [START_REF]SymPy: Python library for symbolic mathematics[END_REF]. They factorize the regular matrix P A, where A ∈ Z n×m (m ≥ n) and P is a n-by-n permutation matrix, into P A = LD -1 U , with a n-by-n lower triangular matrix L, a n-by-n diagonal matrix D, a n-by-m upper trapezoidal matrix U , all having their entries in Z. The matrices L and U are of the form:

L =        p 1 L 21 p 2 . . . . . . . . . L n-1,1 L n-1,2 . . . p n-1 L n,1 L n,2 . . . L n,n-1 1        (1)
and

U =        p 1 U 12 . . . U 1,n-1 U 1,n . . . U 1,m p 2 . . . U 2,n-1 U 2,n . . . U 2,m . . . . . . . . . . . . p n-1 U n-1,n . . . U n-1,m p n . . . U n,m        (2)
The p i are the pivots, and

D = diag(p 1 , p 1 p 2 , . . . , p n-2 p n-1 , p n-1 ) (3)
The pseudo-codes for dense matrices are recalled in A.

If n = m, p n = det A and p i is the principal minor of the top left i-by-i sub-block of A (therefore the positive definiteness can be tested with ∀i, p i > 0).

Concerning the forward substitution, since all the divisions are exact in Z [START_REF] Turner | A simplified fraction-free integer Gauss elimination algorithm[END_REF], one has the property:

∀b ∈ Z n , DL -1 b ∈ Z n (4)
therefore DL -1 gets all its entries in Z.

For the backward substitution, only the first n-by-n sub-block Û of U is used. The property related to the backward substitution is:

∀y = DL -1 b, b ∈ Z n , Û -1 (det Â)y ∈ Z n ( 5
)
where  is the first n-by-n sub-block of A, and det  = p n = Ûn,n .

Case of a singular matrix

The case of a full rank surjective matrix A is trivial and has been considered in [START_REF] Dumas | Dense linear algebra over word-size prime fields: the FFLAS and FFPACK packages[END_REF][START_REF] Ibarra | A generalization of the fast LUP matrix decomposition algorithm and applications[END_REF]. We are concerned herein with the case of a singular square operator (m = n), and the case with more unknowns than equations (m > n) with a possible rank deficiency.

Standard LU factorization of a square singular matrix

When det A = 0, a solvability condition should be satisfied for the existence of a solution for the system Ax = b. If the kernel of the adjoint of A (or the transpose A T of A according to a scalar product on Z n ) is generated by the column vectors stored in matrix S, this solvability condition is

S T b = 0 (6)
Once satisfied, this condition allows to define the set of solutions, up to a vector in the kernel of A:

x = A † b + Ru (7)
R stores independent column vectors generating the kernel of A, and u is the column vector of the coordinates of a solution in the kernel of A. A † is an arbitrary generalized inverse of A [START_REF] Noble | A method for computing the generalized inverse of a matrix[END_REF][START_REF] Farhat | Encyclopedia of Vibration, chapter Linear Algebra[END_REF][START_REF] Golub | Matrix Computations[END_REF] (any generalized inverse generates the same family of solutions).

If S T b = 0, there is no solution.

For sake of simplicity, we assume that a standard LU factorization of A [START_REF] Dopico | Multiple LU factorizations of a singular matrix[END_REF], with a unitary lower triangular matrix L, can be performed without pivoting. Null diagonal entries of U correspond to null pivots; they are as numerous as the size of the kernel, which will be denoted with r [START_REF] Hwang | Rank revealing LU factorizations[END_REF]. The LU factorization can easily be modified to return the factorization of a regularized matrix Ā and to allow the extraction of the kernel vector sets R and S. The corresponding generalized inverse is then A † = Ā-1 . This is described in the following developments.

Consider first an arbitrary splitting of the unknowns in two sets, denoted with subscripts 1 and 2. The matrix A can be split accordingly in several blocks:

A 11 A 12 A 21 A 22 x 1 x 2 = b 1 b 2 (8)
where A 11 and A 22 are diagonal square blocks. If A 11 is regular (which is assumed here), this system can be condensed onto the second set of unknowns:

A 11 A 12 0 A 22 x 1 x 2 = b 1 b 2 (9)
with the Schur complement

A 22 = A 22 -A 21 A -1 11 A 12 and the condensed right hand side b 2 = b 2 - A 21 A -1
11 b 1 . The initial system reads:

I 1 0 A 21 A -1 11 I 2 A 11 A 12 0 A 22 A x 1 x 2 = I 1 0 A 21 A -1 11 I 2 b 1 b 2 b ( 10
)
where I k is the identity matrix (with a size consistent to its use).

The previous condensation operation can also be performed recursively, line per line of the left upper block A 11 of the system (8) to give:

U 11 U 12 0 A 22 x 1 x 2 = c 1 b 2 (11)
and the initial system reads:

L 11 0 L 21 I 2 L (1) U 11 U 12 0 A 22
U (1) A x 1 x 2 = L 11 0 L 21 I 2 c 1 b 2 b ( 12
)
for which the same Schur complement has been computed recursively (this can be easily proved using the sub-block factorization A 11 = L 11 U 11 , A 12 = L 11 U 12 , A 21 = L 21 U 11 and the uniqueness of the factorization without pivoting). Assume now that the second diagonal block, of size r ×r, corresponds to all the null pivots (this would be obtained with a total pivoting strategy): A 22 = 0. In this case, the factorization is completed and L = L (1) , U = U (1) . Existence of a solution is ensured if and only if b 2 = 0, which is the aforementioned solvability condition. This condition reads: b 2 = B T L -1 b = 0, where B ∈ Z n×r is a Boolean matrix containing non null terms only on its second block, which is an identity block of size r × r: B T = 0 I 2 . By identification with (6), the kernel of A T is therefore

S = L -T B ( 13
)
One way to resume the solution of the system is to 'clamp' the undetermined unknowns (the last ones) to an arbitrary fixed value, for instance identity, and to backsolve the remaining regular left upper block. Null pivots in the partial factorization are therefore replaced with identity, and dedicated partially condensed right hand sides (as many as the kernel size) are settled:

U 11 U 12 0 I 2 Ū R = B = 0 I 2 (14)
Ū is then invertible and the solution of this system is a generating set of independent vectors of the kernel of A:

R = Ū -1 B (15)
This is the factorization of a regularized version of A:

Ā = L Ū = L U + 0 0 0 I 2 = LU + L 11 0 L 12 I 2 0 0 0 I 2 = A + 0 0 0 I 2 (16)
A and Ā only differ with the kernel sub-space.

Completely fraction-free factorization of a singular square matrix

The goal of this section is to extend the previous factorization, regularization and kernel extraction to fraction-free computations. The proposed factorization requires a total pivoting strategy [START_REF] Peña | LDU decompositions with L and U well conditioned[END_REF][START_REF] Azulay | Which pivot to solve linear systems[END_REF] (though only the first non null pivot is searched for, not the pivot with maximal amplitude). In this case, the null pivots are postponed to the last degree-of-freedom positions, and a null block is built at the stage where no pivoting is possible anymore. Assuming at a first step that pivoting was not necessary for sake of simplicity, the system to solve is identical to [START_REF] Dureisseix | Information transfer between incompatible finite element meshes: Application to coupled thermo-viscoelasticity[END_REF] with a null block A

(1) 22 = 0, and a regular block U 11 :

U 11 U 12 0 0 x 1 x 2 = c 1 b (1) 2 (17)
and the initial system reads:

L 11 0 L 21 I 2 L D 11 0 0 I 2 -1 D -1 U 11 U 12 0 0 U A x 1 x 2 = LD -1 c 1 b (1) 2 b (18)
The blocks U 11 , U 12 , L 11 , L 21 and D 11 are given from expressions in (1), ( 3) and ( 2). The solvability condition is still b

(1) 2 = 0. To provide a fraction-free factorization and kernel extraction, two approaches can be derived: (i) using a fraction-free regularization Ā of A that allows to fulfill the standard fraction-free factorization, (ii) deal directly with the equivalent rectangular system (once the solvability condition is satisfied):

U 11 U 12 x 1 x 2 = c 1 (19)
Choosing one solution or the other depends on implementation issues. The first solution is described in the following, while the second one is discussed in Section 3.3.

As previously, to resume the factorization of a regularized matrix Ā ∈ Z n×n , the replacement of the null block

A (1)
22 should be performed in a suited way: Referring to the recursion formula on line 20 in algorithm 1, the null pivots should be replaced with the last non null pivot found in the previous factorization of the first block, i.e. p n-r (therefore, for k ≥ n -r + 1, A i,k = 0, and the new block A is not modified by the factorization). Finally, one gets the modified blocks:

L22 = p n-r I r-1 0 0 I r , D22 = p 2 n-r I r-1 0 0 p n-r and Ū22 = p n-r I r (20)
the modified factorization:

L = L 11 0 L 21 L22 , D = D 11 0 0 D22 and Ū = U 11 U 12 0 Ū22 (21)
and one recovers [START_REF] Gotsman | On the computation of null spaces of sparse rectangular matrices[END_REF]:

det Ā = p n-r , Ā = L D-1 Ū = A + 0 0 0 L22 D-1 22 Ū22 = A + 0 0 0 I 2 (22)
with the regularized problem:

Ū x = D L-1 b (23)
The fraction-free forward and backward substitution algorithms are preserved as previously; their usage on the factorization of a regularized matrix ensures that the solutions are still taking their entries in Z.

For kernel construction, one needs some matrices B and C in Z n×r with a regular second block of size r × r:

B = 0 B 22 and C = 0 C 22 (24)
Then, S = L-T DB and R = Ū -1 C. Suitable choices for B 22 and C 22 can be made to preserve fraction-free computations. In a first step, one can choose B 22 = I 2 , Indeed, since D L-1 ∈ Z n×n , see [START_REF] Bareiss | Computational solutions of matrix problems over an integral domain[END_REF], its transpose has also its entries in Z which allows to conclude that S ∈ Z n×r . In a second step, we propose to select

C 22 = (det Ā) D22 . Indeed, L D-1 0 D22 T ∈ Z n×r and since L D-1 Ū R = (det Ā) L D-1 0 D22
T , the solving fractionfree algorithm will provide R ∈ Z n×r by design [START_REF] Bini | Polynomial and matrix computations[END_REF].

If pivoting was necessary, the previous proofs are still valid since the pivoting does not alter the intrinsic properties. The factorization of A reads A = P L D-1 U Q and the factorization of the regularized matrix is Ā = P L D-1 Ū Q. P and Q are permutation matrices. They are orthogonal matrices and the reverse permutation operation is P T and Q T . This is used for notation purpose only since the permutation matrices are not explicitly stored in the pseudo-codes of this article in A and B, conforming to the BLAS implementation standards [START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF]. The associated kernels are:

S = P L-T DB and R = Q T Ū -1 (det Ā)C (25)
In practice, these expressions are used as follows. For the kernel S, once the multiple r right hand sides are stored in B, a call to the forward substitution routine that allows to solve for a transpose left hand side (algorithm 5) is performed. The result is then permuted with P to produce the kernel S. Concerning the kernel R, the multiple r right hand sides stored in C are used. This time, the call to the backward substitution routine (algorithme 3) is performed, and a reverse permutation with Q on the result allows to obtain R.

Case of a rectangular matrix

The extension of the previous algorithm to matrix A ∈ Z n×m with m > n is easy to derive for the standard case [START_REF] Gotsman | On the computation of null spaces of sparse rectangular matrices[END_REF].

For the fraction-free case, we will first consider the case where rank A = n. In this case, the factorization up to the stage [START_REF] Hwang | Rank revealing LU factorizations[END_REF] leads to a regular block U 11 . For an original rectangular system, it reads:

A = A 11 A 12 = LD -1 U 11 U 12 and b = LD -1 c 1 (26)
For any case, the problem is:

U 11 x 1 = c 1 -U 12 x 2 (27)
With a regular block U 11 , the kernel R is such that U 11 R is proportional to U 12 = DL -1 A 12 . With the same argument as before [START_REF] Bini | Polynomial and matrix computations[END_REF],

to get R ∈ Z n×(m-n) , it is proposed to choose R = U -1 11 (det A 11 )U 12 . Indeed, in such a case, LD -1 U 11 R = (det A 11 )A 12 .
With rank A = n, no solvability condition apply and S does not exist.

With a rank-deficient matrix A (rank A = n -r) and a total pivoting strategy, using the same regularization as before, the stage (19) reads:

Ū11 U 12 x 1 x 2 = D L-1 P T b with x 1 x 2 = Q x 1 x 2 (28) or Ū11 x 1 = D L-1 P T b -U 12 x 2 (29)
This problem has the same structure as [START_REF] Noble | A method for computing the generalized inverse of a matrix[END_REF], therefore, the kernel of A T has the same expression:

S = P L-T DB ∈ Z n×r (30)
and the kernel of A is:

R = Q T Ū -1 11 0 0 I d C U 12 0 I ∈ Z m×(r+m-n) ( 31
)
where d is the (n -r, n -r) entry of Ū11 . Once the matrix A has been factorized, and the permutations built, this kernel is obtained by: (i) calling the backward substitution routine with the multiple (r+m-n) right hand sides C U 12 , (ii) the completion of the results with the (m -n) lines d 0 I , and (iii) the reverse permutation of the result.

Implementation issues

Since the regularization procedure is quite easy to code, and does not require memory management for intermediate storage, this is the preferred solution we chose for dealing with the singular case. The treatment of the rectangular case also embeds this regularization procedure when rank A < n. The pseudo-codes in B describe this implementation.

The main drawback of the fraction-free algorithms is the large growing rate of the terms during the factorization and the substitutions (though this is reduced when compared to a division-free algorithm since exact divisions are taken into account here). For instance, for integer computations, Z = Z, the growth in the dynamic range of the integers that must be represented is recalled in [START_REF] Zhou | Fraction-free matrix factors: new forms for LU and QR factors[END_REF]: if the length of a ∈ Z, a = 0, is defined as λ(a) = log |a| + 1 where • is the rounding down to the nearest integer, and if the lengths of the elements of A are bounded by , then the elements of L and U have their lengths bounded by n( + log n).

A classical storage strategy for the integers is to use a double precision declaration and to certify that no cancellation error will occur. With IEEE standard for binary floating-point arithmetic (ANSI-IEEE 754-1985) on 32-bit architectures, to avoid cancellation error, integers stored in REAL*8 should be less than 2 53 -1, i.e. their length should be less than 16. For integer storage the value is limited to 2 31 -1, i.e. their length should be less than 10.

For the case of large-size systems, alternative solutions consist in using arbitrary adaptive precision algebra libraries [START_REF] Granlund | GMP, the GNU multiple precision arithmetic library[END_REF][START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF], modular methods such as a residue number system (RNS) [START_REF] Turner | Fraction-free RNS algorithms for solving linear systems[END_REF], or high-order lifting techniques [START_REF] Storjohann | The shifted number system for fast linear algebra on integer matrices[END_REF].

Though additional gain in dynamic range can be reached by using greatest common divisor searches per completed line of the matrix, this may not improve the worst case.

Test cases

For illustration and eventual checking purposes, we propose two examples for Z = Z. The first one arise from computational geometry: To match incompatible 3D finite element meshes, the core problem for mesh intersection of linear elements is to intersect a line (defined by a segment A 1 , A 2 ) with a plane (defined by a triangle M 1 , M 2 , M 3 ). Though using barycentric coordinates may not be the most efficient approach [START_REF] O'rourke | Computational Geometry in C[END_REF], these coordinates are required for the associated information transfer problem from one mesh to another [START_REF] Dureisseix | Information transfer between incompatible finite element meshes: Application to coupled thermo-viscoelasticity[END_REF]. Therefore, the problem is to find the point M = 3 i=1 λ i M i = 2 j=1 µ j A j with the partition of unity

3 i=1 λ i = 2 j=1 µ j = 1. The problem therefore reduces to find x = λ 2 λ 3 µ 2 T such that Ax = b with A = M 1 M 2 M 1 M 3 -A 1 A 2 and b = [M 1 A 1 ].
For robustness purposes, a snap rounding of nodal coordinates is used; this leads to integer coordinates. Checking singularities is mandatory due to singular positions that occur for instance on planar boundaries of the triangulated domain where segments and triangles are parallel. This is for instance the case when (see Figure 1) One can finally check that any combination of the form ( 7) is a solution, for instance: Ax -d b = 0. In this case, the solution of the intersection problem is λ 1 = 1/2, λ 2 = (48 -64u)/d, λ 3 = (-16 + 64u)/d, µ 1 = (d -64u)/d, µ 2 = 64u/d, with an arbitrary u; this is indeed the whole line (A 1 , A 2 ). Extremal solutions (for which at least one of the barycentric coordinate is in {0, 1} and all are in [0, 1]) can therefore easily be checked. They are the intersections with the triangle edges: (λ

M 1 =   0 0 0   , M 2 =   16 8 0   , M 3 =   0 0 4   , A 1 =   12 6 -1   , A 2 =   -4 -2 3
1 = λ 3 = 1/2, λ 2 = 0, µ 1 = 1/4, µ 2 = 3/4) and (λ 1 = λ 2 = 1/2, λ 3 = 0, µ 1 = 3/4, µ 2 = 1/4).
The next example is extracted from [START_REF] Chang | Error-free computer solution of certain systems of linear equations[END_REF], with a singular matrix obtained by replacing two columns by linear combination of others. It is designed to illustrate the growth in the dynamical range of the entries in Z. Consider the following problem:

A =      

Conclusions

This article proposes an extension of complete fraction-free algorithms, first to square singular matrices in an integral domain, A ∈ Z n×n with det A = 0, providing: (i) the factorization of a regularized matrix Ā ∈ Z n×n whose inverse is a pseudo-inverse of A, and (ii) the fraction-free determination of the kernels of A and of A T . Second, the case with rectangular matrices A ∈ Z n×m with m > n, and exhibiting a rank deficiency, is derived as an other extension. As for the original algorithm, all computations are exact in any integral domain and singular linear systems can be solved within their input domain.

A Pseudo-code for complete fraction-free factorization

The pseudo-codes for dense matrices are recalled in algorithms 1, 2 and 3. Supplementary information from reference [START_REF] Zhou | Fraction-free matrix factors: new forms for LU and QR factors[END_REF] are: the matrices L and U are overwritten on A during the factorization (the diagonal of L is not stored), the vector y is overwritten on b during the forward substitution (where b has been permuted if needed), and the vector x is overwritten on y during the backward substitution; the permutations are stored in a list P giving indexes for row permutations (consistent with the standard LAPACK [START_REF] Anderson | LAPACK Users' Guide[END_REF] factorization routine dgetf2, and BLAS [START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF] substitution routine dtrsm); and a typo correction has been fixed with respect to the original article (line 8 in algorithm 2). All the involved divisions are proved to be exact in Z, and the tests are strict [START_REF] Turner | A simplified fraction-free integer Gauss elimination algorithm[END_REF][START_REF] Zhou | Fraction-free matrix factors: new forms for LU and QR factors[END_REF].

B Pseudo-code for singular complete fraction-free factorization with kernel extraction

The additional codes for singular cases are given in algorithms 4 and 6. Note that the last one, for kernel extraction, deals only with a square matrix, which can be the first n-by-n block of a permuted factorization of a n-by-m (m > n) rectangular matrix. The algorithm 5 is a modified forward substitution with a transpose left hand side. All the divisions are exact in Z due to the property (DL -1 ) T ∈ Z n .

P (k) ← k 6: if A(k, k) = 0 then 7:
Search for the first non null pivot A(kpivot, k), for kpivot

∈ [k + 1, n] 8:
if No pivot found then 9:

Error: Matrix is rank deficient {Integer exact division} 20:

A(i, j) ← (pivot * A(i, j) -Aik * A(k, j))/oldpivot if No pivot found then 10:

{Matrix is rank deficient} {No more singularity here: Apply CFF transformation}

22:

pivot ← A(k, k)

for i = k + 1, . . . , n do

24:

Aik ← A(i, k)

25:

for j = k + 1, . . . , m do 26:

{Integer exact division} 27:

A(i, j) ← (pivot * A(i, j) -Aik * A(k, j))/oldpivot

 1

 1 The in-place factorization of A gives the matrix (the two last rows have been permuted: P is stored as [1 3 3] and Q is stored as[1 2 3] in the pseudo-code provided 1 , conforming to the BLAS standard[START_REF] Blackford | An updated set of basic linear algebra subprograms (BLAS)[END_REF]) some possible duplicates in P or Q arise from the choice of using an incremental permutation strategy for efficiency purposes, see associated pseudo-codes.

Figure 1 :LD - 1 U

 11 Figure 1: Test case of a 3D segment-triangle intersection

  The in-place factorization of Ā gives the matrix (with the permutation of rows corresponding to [1 2 4 4 5] and no column permutation, but with final pivots ln p = [4 5]): is obtained as: x = -14110 108710 -154840 0 0T with d = 11006.

25 :Algorithm 2 1 : 12 :Algorithm 3 1 :Algorithm 4 7 : 8 :

 25211231478 Return: modified A (storing L, D and U ), P Forward substitution: solve y from LD -1 y = b for a permuted b Inputs: A (storing L and D), b 2: n ← size(b) 3: oldpivot ← 1 4: for k = 1, 2, . . . , n -1 do 5: pivot ← A(k, k) 6: for i = k + 1, . . . , n do 7: {Integer exact division} 8: b(i) ← (pivot * b(i) -A(i, k) * b(k))/oldpivot Return: modified b (storing y) Backward substitution: solve x from U x = d y with d = det A Inputs: A (storing U ), y 2: n ← size(y) 3: d ← A(n, n) 4: for i = n, n -1,..., 1 do 5: {Integer exact division} 6: y(i) ← (d * y(i) -n k=i+1 A(i, k) * y(k))/A(i, i) 7: end for 8: Return: modified y (storing x), scaling factor d Completely fraction-free factorization of a dense regularization of A with full pivoting 1: Inputs: A 2: (n, m) ← size(A) 3: n p ← 0 4: oldpivot ← 1 5: for k = 1, 2, . . . , n do 6:P (k) ← k, Q(k) ← k if A(k, k) = 0 thenSearch for the first non null pivot A(kpivot, lpivot), for kpivot ∈ [k + 1, n] and lpivot ∈ [k + 1, n] 9:

32 :Algorithm 5 1 :

 3251 Return: modified A (storing L, D and U of the regularized matrix Ā), P , Q and ln p Modified forward substitution: solve y from L T D -1 y = b Inputs: A (storing L and D), b 2: n ← size(b) 3:for i = n, n -1,..., 1 do ) ← oldpivot * b(i) -( n k=i+1 A(k, i) * b(k))/A(i, i) 11: end for 12: Return: modified b (storing y) Algorithm 6 Find kernels from a fraction-free factorization of a regularized square matrix A 1: Inputs: A (storing L, D, U ), ln p (for original nul pivot locations) 2: n ← size(A, 1) 3: n p ← length(ln p ) 4: d ← A(n, n) 5: {Initialize right-hand-sides for kernel extractions} 6: for i = 1, 2, . . . , n p do 7: R(ln p (i), i) ← d 8: S(ln p (i), i) ← 1 9: end for 10: Solve U R = d b where the n p right-hand-sides b are stored in R on input, with fraction-free backward substitution (algorithm 3) 11: Solve L T D -1 S = b where the n p right-hand-sides b are stored in S on input, with fraction-free forward modified substitution (algorithm 5) 12: Return: kernels R and S

  Algorithm 1 Completely fraction-free factorization of dense A 1: Inputs: A 2: (n, m) ← size(A) 3: oldpivot ← 1 4: for k = 1, 2, . . . , n do

	5:

  p ← n p + 1, ln p (n p ) ← k

	11:
	12:	{Apply a suited regularization}
	13:	A(k, k) ← oldpivot
	14:	else
	15:	{Row interchange}
	16:
	17:	{Column interchange}
	18:
	19:	end if
	20:	end if
	21:

n A(k, 1 : m) ↔ A(kpivot, 1 : m), P (k) ← kpivot A(1 : n, k) ↔ A(1 : n, lpivot), Q(k) ← kpivot


Pierre-Olivier Goffard

Claude Lefèvre

Boundary Crossing of Order Statistics Point Processes

Keywords: 60G40, 12E10 Order statistic property, lower and upper boundaries, Appell and Abel-Gontcharoff polynomials

This paper is concerned with the first crossing of an order statistic point process through general moving boundaries. Our purpose is to determine exact boundary crossing probabilities in both one and two boundary cases. Simple recursive methods are obtained that exploit an underlying algebraic structure of polynomial type. This structure is a direct consequence of the order statistic property. The proposed approach is easy to implement and efficient. Perspectives in statistics are also announced.

MSC 2010: 60G55,

Introduction

Many problems in stochastic modeling come down to study the crossing time of a certain stochastic process through a given boundary, lower or upper. Typical fields of application are in risk theory, epidemic modeling, queueing, reliability and sequential analysis. The present paper is concerned with the well-known class of point processes with the order statistic property (OSPP: Order Statistic Point Processes). Definition 1.1. A point process {N (t), t ≥ 0} with N (0) = 0 is an OSPP if for every n ≥ 1, provided P (N (t) = n) > 0, then conditioned upon (N (t) = n), the successive jump times (T 1 , T 2 , . . . , T n ) are distributed as the order statistics of n i.i.d. random variables with distribution function F t (x) supported on [0, t]. So, for any fixed time t, the only knowledge of the distribution of N (t) and the distribution function F t (s), for 0 ≤ s ≤ t, supplies a complete description of the process over [0, t]. Several authors have characterized the class of OSPP under certain conditions. The reader is referred to Crump [START_REF] Crump | On point processes having an order statistic structure[END_REF], Feigin [START_REF] Feigin | On the characterization of point processes with the order statistic property[END_REF] and Puri [START_REF] Puri | On the characterization of point processes with the order statistic property without the moment condition[END_REF]. We focus our interest on the first-crossing time of an OSPP through a one-sided boundary, lower h α (t) or upper h β (t), and the first-exit time from a two-sided region delimited by two boundaries h α (t) and h β (t) that are parallel. The boundaries are allowed to be moving, linearly or not. They are non-decreasing, which is not restrictive. We write them under the form

h α (t) = h(t) -α, h β (t) = h(t) + β, t ≥ 0,
where α and β are two real numbers such that α > 0 and β ≥ 0 and h(t) is a non-decreasing function with h(0) = 0.

In the one boundary case, we will obtain the exact distribution of the first meeting time τ α in the lower boundary h α :

τ α = inf{t ≥ 0 : N (t) = h α (t)}, and the first crossing time τ β through the upper boundary h β :

τ β = inf{t ≥ 0 : N (t) ≥ h β (t)}.
For a region delimited by two such boundaries, we assume that h α and h β are parallel (i.e. share the same function h), and we derive the distribution of the first exit time from the region:

τ α,β = min(τ α , τ β ),
and the first exit time by the lower boundary h α without crossing the upper boundary h β :

τ * α,β = inf{t ≥ 0 : N (t) = h α (t) and N (s) < h β (s) for s < t}.

Crossing problems for Poisson and compound Poisson processes have received a considerable attention, especially for linear boundaries. We mention e.g. Zacks [START_REF] Zacks | Distributions of stopping times for Poisson processes with linear boundaries[END_REF], Gallot [START_REF] Gallot | Absorption and first-passage times for a compound poisson process in a general upper boundary[END_REF], Picard and Lefèvre [START_REF] Lefèvre | First crossing of basic counting processes with lower non-linear boundaries: a unified approach through pseudopolynomials[END_REF], Lehmann [START_REF] Lehmann | Boundary crossing probabilities of Poisson counting processes with general boundaries[END_REF], Perry et al. [START_REF] Perry | Contributions to the theory of firstexit times of some compound processes in queueing theory[END_REF], Stadje and Zacks [START_REF] Stadje | Upper first-exit times of compound Poisson processes revisited[END_REF], Ignatov and Kaishev [START_REF] Ignatov | A finite-time ruin probability formula for continuous claim severities[END_REF], Lefèvre [START_REF] Lefèvre | Discrete compound Poisson process with curved boundaries: polynomial structures and recursions[END_REF], Lefèvre and Loisel [START_REF] Lefèvre | Finite-time ruin probabilities for discrete, possibly dependent, claim severities[END_REF], Lefèvre and Picard [START_REF] Lefèvre | A new look at the homogeneous risk model[END_REF], Xu [START_REF] Xu | First exit times of compound Poisson processes with parallel boundaries[END_REF]. These works were motivated by various topics including the ruin in insurance, the final size in epidemics, the busy period in queues or the detection of changes.

To obtain crossing probabilities for these processes, a method frequently used is to work with Laplace transforms. An alternative approach, proposed intially by Picard and Lefèvre [START_REF] Lefèvre | First crossing of basic counting processes with lower non-linear boundaries: a unified approach through pseudopolynomials[END_REF], consists in applying simple recursive relations to calculate the probabilities. This is such a method that we choose to follow here. The recursions derived for the computations will rely on the existence of an underlying polynomial structure in the probabilities. The existence of this algebraic structure is a consequence of the order statistic property.

The formulas have been implemented in Mathematica. The source code is provided online at [START_REF] Goffard | Online accompaniment for "Boundary Crossing of Order Statistics Point Process[END_REF]. The recursive relations contribute to an efficient evaluation in terms of computing time. We invite the reader to appreciate the impact on the crossing probabilities when tuning the parameters of the OSPP or changing the shape of the boundary. A Monte Carlo estimator of the crossing probabilities is defined from our formulas and compared to a classical Monte Carlo estimator. The new Monte Carlo estimator beats the classical one in terms of variance. This paves the way for potential applications in a statistical framework.

The paper is organized as follows. Section 2 presents the class of OSPP under consideration and the families of polynomials used as a key tool in the study. Section 3 deals with the first meeting problem of an OSPP in a lower boundary. Section 4 examines the first crossing problem of an OSPP through an upper boundary. Section 5 considers the first exit problem for an OSPP that is trapped inside two boundaries. Section 6 discusses a Monte Carlo estimator constructed from the results derived before.

Order statistic point processes

A complete characterization of the OSPP's was obtained by Puri [START_REF] Puri | On the characterization of point processes with the order statistic property without the moment condition[END_REF], following on earlier work. His main result is stated below (see [START_REF] Puri | On the characterization of point processes with the order statistic property without the moment condition[END_REF]Theorem 5]). Proposition 2.1. Let {N (t), t ≥ 0} be an OSPP where µ(t) = E[N (t)] is finite for all t.

1. If lim t→∞ µ(t) = ∞, there exists a Poisson process {P(t), t ≥ 0} with rate 1, an independent nonnegative random variable W , both on the same space probability as is the OSPP, and a time transformation ν(t) such that

N (t) = P [W ν(t)] , t ≥ 0, a.s. (2.1)
In other words, {N (t), t ≥ 0} is a mixed Poisson process up to a time-scale transformation.

2. If lim t→∞ µ(t) = γ < ∞, there exists a process {D Z (t), t ≥ 0} that counts the deaths during (0, t] in a process with initially Z individuals and i.i.d. lifetimes of distribution function µ(t)/γ, t ≥ 0, and an independent non-negative integer-valued random variable Z, both on the same probability space as is the OSPP, such that

N (t) = D Z (t), t ≥ 0, a.s. (2.
2)

The process {D Z (t), t ≥ 0} is named a mixed sample process.

In both cases, the order statistic property holds over (0, t] with

F t (s) = µ(s)/µ(t), 0 ≤ s ≤ t. (2.3)
Note that E[N (t)] = µ(t) implies that E(Z) = γ in (2.2). As an illustration, a few particular cases of OSPP of special interest are highlighted.

Special cases.

(i) A Poisson process of parameter λ > 0. N (t) has a Poisson distribution of intensity λt and

F t (s) = s/t, 0 ≤ s ≤ t.
Formula (2.1) holds with ν(t) = t and W = λ.

(ii) A Pólya-Lundberg process of parameters λ > 0 and b ≥ 0, i.e. a non-homogeneous linear birth process of rate λ n (t) = λ(1 + bn)/(1 + λt). N (t) has a negative binomial distribution:

P [N (t) = n] = n -1 + 1/b n λbt 1 + λbt n 1 1 + λbt 1/b , n ≥ 0.
and, here too, F t (s) = s/t, for 0 ≤ s ≤ t. Formula (2.1) holds with ν(t) = t and W = d Γ(1/b, λb), where Γ(r, m) denotes the gamma distribution with shape parameter r > 0 and mean parameter m > 0.

(iii) A linear birth process with immigration of birth rate b > 0 and immigration rate λ ≥ 0. N (t) has a negative binomial distribution:

P [N (t) = n] = e -λt λ/b + n -1 n (1 -e -bt
) n , n ≥ 0, and F t (s) = (e bs -1)/(e bt -1), 0 ≤ s ≤ t.

Formula (2.1) holds with ν(t) = e bt -1 and W = d Γ(λ/b, 1).

(iv) A death counting process in a linear death process of rate b > 0 and initial size z ≥ 1. N (t) has a binomial distribution:

P [N (t) = n] = z n (1 -e -bt ) n e -bt(z-n) , 0 ≤ n ≤ z, and
F t (s) = (1 -e -bs )/(1 -e -bt ), 0 ≤ s ≤ t.
Formula (2.2) holds with µ(t) = z(1 -e -bt ), hence γ = z, and Z = z almost surely.

Order statistics and polynomial structures

The first-crossing or exit problem of an OSPP will be formulated in terms of the order statistics (U 1:n , . . . , U n:n ) for a sample of n uniform random variables on (0, 1). It is known that the joint distributions of (U 1:n , . . . , U n:n ) rely on an underlying polynomial structure (see Lefèvre and Picard [START_REF] Lefèvre | Risk models in insurance and epidemics: A bridge through randomized polynomials[END_REF] and the references therein). A few key points are recalled below.

One-sided joint distributions

Let U = {u i , i ≥ 1} be a sequence of reals, non-decreasing in our context. To U is attached a (unique) family of Appell polynomials of degree n in x, {A n (x|U ), n ≥ 0} defined as follows.

Definition 2.2. Starting with A 0 (x|U ) = 1, the A n (x|U )'s satisfy the differential equations

A (1) n (x|U ) = nA n-1 (x|U ), (2.4)
with the border conditions

A n (u n |U ) = 0, n ≥ 1. (2.5)
So, each A n has the integral representation

A n (x|U ) = n! x un yn u n-1
dy n-1 . . .

y 1 u 1 dy 1 dy n , n ≥ 1. (2.6)
In parallel, to U is attached a (unique) family of Abel-Gontcharov (A-G) polynomials of degree n in x, {G n (x|U ), n ≥ 0}.

G (1) n (x|U ) = nG n-1 (x|EU ), (2.7)
where EU is the shifted family {u i+1 , i ≥ 1}, and with the border conditions

G n (u 1 |U ) = 0, n ≥ 1. (2.8)
So, each G n has the integral representation

G n (x|U ) = n! x u 1 y 1 u 2 dy 2 . . . y n-1 un dy n dy 1 , n ≥ 1.
(2.9)

The Appell and A-G polynomials are closely related through the identity

G n (x|u 1 , . . . , u n ) = A n (x|u n , . . . , u 1 ), n ≥ 1.
(2.10)

The two families (i.e. for all n ≥ 0), however, are distinct and enjoy different properties. From (2.6) and (2.9), it is clear that the polynomials A n and G n can be interpreted in terms of the joint distribution of the vector (U 1:n , . . . , U n:n ).

Proposition 2.4. For 0 ≤ u 1 ≤ . . . ≤ u n ≤ x ≤ 1, P [U 1:n ≥ u 1 , . . . , U n:n ≥ u n and U n:n ≤ x] = A n (x|u 1 , . . . , u n ), n ≥ 1.
(2.11)

For 0 ≤ x ≤ u 1 ≤ . . . ≤ u n ≤ 1, P [U 1:n ≤ u 1 , . . . , U n:n ≤ u n and U 1:n ≥ x] = (-1) n G n (x|u 1 , . . . , u n ), n ≥ 1.
(2.12)

The representations (2.11) and (2.12) will play a key role for the first-crossing of an OSPP through an upper and the first meeting of an OSPP and a lower boundary respectively. Numerically, it will be necessary to evaluate some special values of the polynomials. To this end, it is convenient to use the following recusive relations.

Proposition 2.5. The Appell polynomials are computed through the expansion

A n (x|U ) = n k=0 n k A n-k (0|U )x k , n ≥ 1, (2.13)
where the A n (0|U )'s are obtained recursively from

A n (0|U ) = - n k=1 n k A n-k (0|U )u k n , n ≥ 1. ( 2

.14)

The A-G polynomials are computed through the recursion

G n (x|U ) = x n - n-1 k=0 n k u n-k k+1 G k (x|U ), n ≥ 1. (2.15)
Formulas (2.13), and (2.14) follow from the Taylor's expansion of A n , using also (2.4), and (2.5). Formula (2.15) follows from an Abelian expansion of x n based on (2.7), and (2.8). Details are omitted here. Of course, the computing time increases with the degree of the polynomials.

Note that

A n (x|a + bU ) = b n A n ((x -a)/b |U ) , n ≥ 1, (2.16)
with the same identity for G n . A simple special case is when the sequence U is affine.

Proposition 2.6.

If U = {u i = a + b(i -1), i ≥ 1} for some reals a, b, A n (x|U ) = (x -a -b(n -1))(x -a + b) n-1 , (2.17)
G n (x|U ) = (x -a)(x -a -bn) n-1 , n ≥ 1. (2.18)
Formula (2.17) is proven by induction from (2.6). Formula (2.18) is then a consequence of (2.10). In addition to speeding up the computations, the identities (2.17) and (2.18) will allow us to link our results to some classical ones in the literature.

Two-sided joint distributions

Let U = {u i , i ≥ 1} and V = {v i , i ≥ 1} be two sequences of non-decreasing reals such that 0 ≤ u 1 ≤ . . . ≤ u n ≤ 1, 0 ≤ v 1 ≤ . . . ≤ v n ≤ 1 and u 1 ≤ v 1 , . . . , u n ≤ v n .
Consider again the vector (U 1:n , . . . , U n:n ). The joint distributions of interest are rectangular probabilities of the form

d n (U, V ) = P [u 1 ≤ U 1:n ≤ v 1 , . . . , u n ≤ U n:n ≤ v n ] (2.19) = n! vn un v n-1 ∧yn u n-1 dy n-1 . . . v 1 ∧y 2 u 1 dy 1 dy n .
A polynomial structure can be exhibited in this case too. The trick consists in varying the sample size, k say, of uniforms on (0, 1) from 1 to n. Specifically, the following rectangular probabilities are defined

d k (U, V ) = P [u 1 ≤ U 1:k ≤ v 1 , . . . , u k ≤ U k:k ≤ v k ], 1 ≤ k ≤ n. (2.20) Proposition 2.7. For u n ≤ x ≤ 1, P [u 1 ≤ U 1:n ≤ v 1 , . . . , u n ≤ U n:n ≤ v n ∧ x] = n k=0 n k (v k+1 -x) n-k + (-1) n-k d k . (2.21) For instance, suppose that u n ≤ v n-m for some m ≥ 1. Consider the intervals [v n-i , v n-i+1 ] for i = 0, . . . , m, with v n+1 ≡ 1. For v n-i ≤ x ≤ v n-i+1 , (2.21) then reduces to P [u 1 ≤ U 1:n ≤ v 1 , . . . , u n ≤ U n:n ≤ v n ∧ x] = n k=n-i n k (x -v k+1 ) n-k d k ,
which is a polynomial in x of degree i. So, for u n ≤ x ≤ 1, (2.21) shows that the sought probabilty is a piecewise polynomial function. It has been named of Sheffer type by Niederhausen [START_REF] Niederhausen | Sheffer polynomials for computing exact Kolmogorov-Smirnov and Rényi distributions[END_REF]. Now, let us take x = u n in (2.21). Evidently, the probability in the l.h.s. becomes 0. This provides us with the simple formula (2.22) below.

Proposition 2.8. Let d 0 (U, V ) = 1. The probabilities d n (U, V ) are computed through the recursion

d n (U, V ) = - n-1 k=0 n k (v k+1 -u n ) n-k + (-1) n-k d k (U, V ), n ≥ 1. (2.22)
3 First-meeting with a lower boundary

Let {N (t), t ≥ 0} be an OSPP that jumps at times T n , n ≥ 1 (see Definition 1.1). We recall that given (N (t) = n), the vector (T 1 , . . . , T n ) is distributed as the order statistics of n i.i.d. random variables with distribution function F t (s), for 0 ≤ s ≤ t. Consider for this process a lower boundary of the form h α (t) = h(t) -α, where α is a non-negative real and h is a non-decreasing function with h(0) = 0. We want to determine the distribution of the first-meeting time τ α of the OSPP and this boundary. Let

h -1 α (x) = inf{t ≥ 0 : h α (t) ≥ x} denote the generalized inverse of h α .
As the jumps of the process are of unitary height, the possible meeting levels are the naturals 0, 1, 2, . . .. The corresponding meeting times are the instants α n = h -1 α (n), n ≥ 0, for which the boundary is integer-valued. This is illustrated in Figure 1 where the meeting happens at τ α = α 5 .

P (τ α = α n ) = (-1) n P [N (α n ) = n] G n [0|{F αn (α 0 ), . . . , F αn (α n-1 )}] , n ≥ 0, (3.1)
where G n [0|{. . .}] is an A-G polynomial such as defined in Subsection 2.2.1.

Proof. Obviously, P (τ α = α 0 ) = P [N (τ α ) = 0]. For n ≥ 1, the event {τ α = α n } can be expressed

{τ α = α n } = {T 1 ≤ α 0 } ∩ {T 2 ≤ α 1 } ∩ . . . ∩ {T n ≤ α n-1 } ∩ {T n+1 ≥ α n } = n i=1 {T i ≤ α i-1 } ∩ {N (α n ) = n}.
Conditioning on {N (α n ) = n}, we then have

P (τ α = α n ) = P [N (α n ) = n] P n i=1 {T i ≤ α i-1 }|N (α n ) = n . (3.2)
By the order statistic property, we know that

P n i=1 {T i ≤ α i-1 }|N (α n ) = n = P n i=1 {V i:n (α n ) ≤ α i-1 } , (3.3)
where (V 1:n (α n ), . . . , V n:n (α n )) are the order statistics for n i.i.d. random variables of distribution function F αn defined on (0, α n ). Note that

F αn (V i:n (α n )) = d U i:n , 1 ≤ i ≤ n,
where (U 1:n , . . . , U n:n ) are the order statistics for n independent uniforms on (0, 1). Thus, (3.3) may be rewritten as

P n i=1 {T i ≤ α i-1 }|N (α n ) = n = P n i=1 {F αn (V i:n (α n )) ≤ F αn (α i-1 )} = P n i=1 {U i:n ≤ F αn (α i-1 )} = (-1) n G n [0|{F αn (α 0 ), . . . , F αn (α n-1 )}] . (3.4)
The last equality follows from formula (2.12) for the A-G polynomials. Combining (3.2), (3.3), (3.4) yields the announced result (3.1).

Formula (3.1) has been implemented via Mathematica. The code is available online at [START_REF] Goffard | Online accompaniment for "Boundary Crossing of Order Statistics Point Process[END_REF]. The values of G n [0|{. . .}] are computed using the recursion (2.15). The characterization of an OSPP (see Proposition 2.1) enables us to provide a more explicit version of (3.1).

Proposition 3.2. 1. If lim t→∞ µ(t) = ∞, P (τ α = α n ) = (-1) n n! G n [0|{ν(α 0 ), . . . , ν(α n-1 )}] E[W n e -W ν(αn) ], n ≥ 0. (3.5) 2. If lim t→∞ µ(t) = γ < ∞, P (τ α = α n ) = (-1) n γ n G n [0|µ(α 0 ), . . . , µ(α n-1 )] E Z n 1 - µ(α n ) γ Z-n
, n ≥ 0.

P (τ α = α n ) = (-λ) n n! e -λαn G n (0|α 0 , . . . , α n-1 ), n ≥ 0,
(ii) a Pólya-Lundberg process,

P (τ α = α n ) = (-λb) n n -1 + 1/b n 1 1 + λbα n n+1/b
G n (0|α 0 , . . . , α n-1 ), n ≥ 0, (iii) a linear birth process with immigration,

P (τ α = α n ) = (-1) n e -(λ+bn)αn λ/b + n -1 n G n [1|{e bα 0 , . . . , e bα n-1 }], n ≥ 0,
(iv) a linear death counting process,

P (τ α = α n ) = e -(z-n)bαn z n G n [1|{e -bα 0 , . . . , e -bα n-1 })], 0 ≤ n ≤ z.
A simple situation is when {N (t), t ≥ 0} is a Poisson, or mixed Poisson, process and the lower boundary is linear.

Corollary 3.4. For a (mixed) Poisson process, if h α = ct -α with c > 0, P (τ α = α n ) = α 0 α n P [N (α n ) = n], n ≥ 0. (3.7)
Proof. If the OSPP is a mixed Poisson process, then F t (s) = s/t, 0 ≤ s ≤ t, so that (3.1) becomes

P (τ α = α n ) = (-1) n P [N (α n ) = n] (1/α n ) n G n [0|{α 0 , . . . , α n-1 }] , n ≥ 0. (3.8)
Now, the function h α being linear, its inverse is h -1 α (x) = (α + x)/c. Thus, in (3.8), α n = (α + n)/c, n ≥ 0, i.e. α n is an affine function of n. In such a case, G n is given by formula (2.18), which yields

G n [0|{α 0 , . . . , α n-1 }] = G n [0|{(α + i)/c, 0 ≤ i ≤ n -1}] = (-1) n (α/c)[(α + n)/c] n-1 = (-1) n α 0 α n-1 n . (3.9)
Inserting (3.9) in (3.8) then gives (3.7).

Formula (3.7) corresponds to a Kendall type identity (see e.g. Borovkov and Burq [START_REF] Borovkov | Kendall's identity for the first crossing time revisited[END_REF]). In fact, our main formula (3.1) can be viewed as a generalization of this identity to the first meeting problem of an OSPP through a moving lower boundary.

First-crossing through an upper boundary

Consider now the case of an upper boundary of the form h β (t) = h(t)+β, where β is a positive real and h is a non-decreasing function with h(0) = 0. We denote by β n = h -1 β (n), n ≥ 0, the instant at which the upper barrier reaches the (integer-valued) level n, where h -1 β is the generalized inverse of h β . We want to determine the distribution of the first-crossing time τ β of the OSPP through this boundary. Note that such a crossing will not occur through a meeting but by a jump over the boundary. This is illustrated in Figure 2 where the crossing occurs at τ β = T 5 . As the crossing is not a meeting, it is not convenient to work directly with the density function of τ β . We will derive a formula for the survival function instead. Denote by 1 A the indicator of an event A and by x the integer part of a real x.

P (τ β > t) = E A N (t) 1 {F t (β 1 ), . . . , F t (β N (t) )} 1 {N (t)≤ h β (t) } , t ≥ 0, (4.1)
where A n (1|{. . .}) is an Appell polynomial such as defined in Subsection 2.2.1.

Proof. On the set {τ β > t}, the number of jumps of the process up to time t cannot exceed the level h β (t) . Thus, we have

{τ β > t} = h β (t) n=0 {N (t) = n} ∩ {τ β > t} = h β (t) n=0 {N (t) = n} ∩ {∩ n i=1 {T i > β k }} .
Conditionning on {N (t) = n}, we then get

P (τ β > t) = h β (t) n=0 P n i=1 {T i > β i }|N (t) = n P [N (t) = n], (4.2)
and by the order statistic property,

P n i=1 {T i > β i }|N (t) = n = P n i=1 {V i:n (t) > t i } , (4.3)
where (V 1:n (t), . . . , V n:n (t)) are the order statistics for n i.i.d. random variables of distribution function

F t . Since F t [V i:n (t)] = d U i:n , 1 ≤ i ≤ n,
where (U 1:n , . . . , U n:n ) are the order statistics for n independent uniforms on (0, 1), formula (4.3) may be rewritten as

P n i=1 {T i > β i }|N (t) = n = P n i=1 {F t [V i:n (t)] > F t (β i )} = P n {U i:n > F t (β i )} = A n 1 {F t (β 1 ), . . . , F t (β n )} . (4.4)
Reinjecting (4.4) into (4.2) yields

P (τ β > t) = h β (t) n=0 A n [1 {F t (β 1 ), . . . , F t (β n )}] P [N (t) = n],
which leads to the result (4.1).

Here too, formula (4.1) has been implemented via Mathematica. The code is available online at [START_REF] Goffard | Online accompaniment for "Boundary Crossing of Order Statistics Point Process[END_REF]. Moreover, it is important to underline that the closed-form expression of (4.1) is quite suitable for an estimation by Monte Carlo simulation. In Section 6 we will show numerically that the Monte Carlo simulation of (4.1) leads to a less volatile estimator than the classical Monte Carlo estimator.

From Proposition 2.1 and using the identity (2.16), we can then deduce a more explicit version of formula (4.1). Corollary 4.2.

1. If lim t→∞ µ(t) = ∞,

P (τ β > t) = h β (t) n=0 A n [ν(t)|{ν(β 1 ), . . . , ν(β n )}] E W n e -W ν(t) n! , t ≥ 0. (4.5) 2. If lim t→∞ µ(t) = γ < ∞, P (τ β > t) = h β (t) n=0 1 γ n A n [µ(t)|{µ(β 1 ), . . . , µ(β n )}] E Z n 1 - µ(t) γ Z-n , t ≥ 0. (4.6)
Let us examine the special cases (i)-(iv) presented in Section 2.1.

Corollary 4.3. (i) For a Poisson process,

P (τ β > t) = h β (t) n=0 λ n e -λt n! A n [t {β 1 , . . . , β n }], t ≥ 0,
(ii) a Pólya-Lundberg process,

P (τ β > t) = h β (t) n=0 1/b + n -1 n λb 1 + λbt n 1 1 + λbt 1/b A n [t {β 1 , . . . , β n }], t ≥ 0,
(iii) a linear birth process with immigration, λ+bn) A n [e bt -2 {e bβ 1 , . . . , e bβn }], t ≥ 0, (iv) a linear death counting process,

P (τ β > t) = h β (t) n=0 (-1) n λ/b + n -1 n e -t(
P (τ β > t) = h β (t) n=0 (-1) n e -(z-n)bt z n A n [e -bt -2 {e -bβ 1 , . . . , e -bβn }], t ≥ 0.
Formula (4.1) becomes extremely simple in the particular case where {N (t) ; t ≥ 0} is a Poisson, or mixed Poisson, process and the upper boundary is a diagonal line. Let x + denote the positive part of a real x.

Corollary 4.4. For a (mixed) Poisson process, if h β (t) ≡ h 0 (t) = ct with c > 0,

P (τ β > t) = E [1 -N (t)/ct] + , t ≥ 0., (4.7)
Proof. The OSPP being a mixed Poisson process, we have F t (s) = s/t, 0 ≤ s ≤ t, and since h β (t) ≡ h 0 (t) = ct, then β n ≡ 0 n = n/c, n ≥ 0. Thus, in (4.1),

A N (t) 1 {F t (β 1 ), . . . , F t (β N (t) )} = A N (t) [1| {i/ct, i ≥ 1}] = 1 -N (t)/ct, (4.8)
by virtue of formula (2.17). Inserting (4.8) in (4.1) then yields formula (4.7).

A result of this form has a long history and different proofs are available. In particular, it has has been derived in the framework of ballot type theorems (e.g. Takács [START_REF] Takaćs | Combinatorial Methods in the Theory of Stochastic Processes[END_REF]). In risk theory, it corresponds to the finite time non-ruin probability in the compound Poisson risk model with no initial reserves and claim sizes equal to 1 (e.g. Asmussen and Albrecher [1, Theorem 2.1]).

First-exit with two parallel boundaries

Finally, we suppose that the OSPP is trapped into a region bounded from below by a lower boundary h α (t) = h(t) -α and from above by an upper boundary h β (t) = h(t) + β, where α > 0, β ≥ 0 and the function h is non-decreasing with h(0) = 0. Observe that h α and h β are built with the same function h, i.e. the two boundaries are parallel. As before, we use the notation

α n = h -1 α (n), n ≥ 0, and β n = h -1 β (n), n ≥ 1.
Our aim is to derive the distribution of the first exit-time from this region, τ α,β say, and of the first-exit time through the lower boundary, τ * α,β say. This is illustrated in Figure 3 where τ α,β = τ * α,β = α 5 . The assumption of parallel boundaries will allow us to deal with the class of OSPP processes. Let us mention that Xu [START_REF] Xu | First exit times of compound Poisson processes with parallel boundaries[END_REF] has examined the first-exit problem of a compound Poisson process with parallel linear boundaries. Applications of stopping times with parallel boundaries are also provided in that paper. Proposition 5.1. For n ≥ 1,

P τ * α,β = α n = d n [{F αn (β 1 ), . . . , F αn (β n )}, {F αn (α 0 ), . . . , F αn (α n-1 )}] P [N (α n ) = n], (5.1)
and for t ≥ 0,

P (τ α,β > t) = E d N (t) {F t (β 1 ), . . . , F t (β N (t) )}, {F t (α 0 ), . . . , F t (α N (t)-1 )} 1 { hα(t) ≤N (t)≤ h β (t) } , (5.2
) where d n [{. . .}, {. . .}] is a rectangular probability such as defined in Subsection 2.2.2.

Proof. We proceed as for Propositions 3.1 and 4.1. First, consider the event {τ * α,β = α n }. For n = 0, it means {N (α 0 ) = 0}. For n ≥ 1, it is equivalent to

{τ α = α n } ∩ {τ β > α n } = {β 1 ≤ T 1 ≤ α 0 } ∩ . . . ∩ {β n ≤ T n ≤ α n-1 } ∩ {T n+1 ≥ α n } = n i=1 {β i ≤ T i ≤ α i-1 } ∩ {N (α n ) = n}.
By conditioning on {N (α n ) = n} and using the order statistic property, we then have

P τ * α,β = α n = P [N (α n ) = n] P n i=1 {β i ≤ T i ≤ α i-1 }|N (α n ) = n = P [N (α n ) = n] P n i=1 {β i ≤ V i:n (α n ) ≤ α i-1 } , (5.3)
where (V 1:n (α n ), . . . , V n:n (α n )) are the order statistics for n i.i.d. random variables of distribution function F αn . Let us apply the transform F αn to get the order statistics (U 1:n , . . . , U n:n ) for n independent uniforms on (0, 1). Thus, we can write

P n i=1 {β i ≤ V i:n (α n ) ≤ α i-1 } = P n i=1 {F αn (β i ) ≤ U i:n ≤ F αn (α i-1 )} = d n [{F αn (β 1 ), . . . , F αn (β n )}, {F αn (α 0 ), . . . , F αn (α n-1 )}] , (5.4)
by virtue of definition (2.20). Combining (5.3) and (5.4) gives the result (5.1). Now, consider the event {τ α,β > t} = {τ α > t} ∩ {τ β > t}. We see that it can be expressed as

{τ α,β > t} = h β (t)
n= hα(t)

{N (t) = n} ∩ {τ α > t} ∩ {τ β > t} = h β (t)
n= hα(t)

{N (t) = n} ∩ n i=1 {β i ≤ T i ≤ α i-1 } .
Arguing as above, we then obtain

P (τ α,β > t) = h β (t) n= hα(t) P n i=1 {β i ≤ T i ≤ α i-1 }|N (t) = n P [N (t) = n] = h β (t) n= hα(t) d N (t) {F t (β 1 ), . . . , F t (β N (t) )}, {F t (α 0 ), . . . , F t (α N (t)-1 )} P [N (t) = n],
which is formula (5.2).

First-exit problems with two boundaries have been less studied in the literature. In the case of the (compound) Poisson process, some interesting papers are e.g. Xu [START_REF] Xu | First exit times of compound Poisson processes with parallel boundaries[END_REF] cited before, Perry et al. [START_REF] Perry | A two-sided first-exit problem for a compound Poisson process with a random upper boundary[END_REF] (in which the upper boundary is random) and Lehmann [START_REF] Lehmann | Boundary crossing probabilities of Poisson counting processes with general boundaries[END_REF] (in which the boundaries are not parallel). The approach in these works is different from that followed here and relies on the specific structure of the Poisson process. For the general class of OSPP, the question of two arbitrary boundaries requires further investigation.

An induced simulation study

Our aim in this Section is to discuss an evaluation by simulation of the survival function of τ β , the first-crossing time through an upper boundary. We recall that by Proposition 4.1,

P (τ β > t) = E A N (t) 1 {F t (β 1 ), . . . , F t (β N (t) )} 1 {N (t)≤ h β (t) } ≡ E p N (t) say, t ≥ 0. (6.1)
The survival function being expressed as an expectation, a Monte Carlo simulation is suitable for its estimation. We name it an Appell Polynomial Monte Carlo (APMC) procedure. We will show below the superiority of an APMC simulation over a Crude Monte Carlo (CMC) procedure. Note that a priori, an approximation by simulation could be viewed of limited relevance as formula (6.1) can be evaluated by recursion. Nevertheless, we will point out that our study paves the way for potential statistical applications.

For our problem, the CMC technique consists in writing the survival function of τ β as

P (τ β > t) = E 1 τ β >t , t ≥ 0. (6.
2)

The expectation in (6.2) is then approximated by the common Monte Carlo estimator based on trajectories drawn for the OSPP up to time t. For that, replications of N (t) are simulated and the jump times are then given by the order statistics of a sample of size N (t) random variables with distribution function F t . The CMC evaluation of (6.2) falls in the rare event simulation framework where the quantity of interest is of the same magnitude as the variance of the CMC estimator, which makes it inefficient from a point of view of the committed relative error. For a detailed account on the CMC procedure and rare event simulations, see e.g. the book of Asmussen and Glynn [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]. Denote p t = P (τ β > t). The variance of the CMC estimator is given by

σ 2 CMC = p t (1 -p t ). (6.3)
By comparison, the variance of the APMC estimator is

σ 2 AP M C = V ar A N (t) 1 {F t (β 1 ), . . . , F t (β N (t) )} 1 {N (t)≤ h β (t) } = h(t)+β n=0 A n [1 {F t (β 1 ), . . . , F t (β n )}] 2 P [N (t) = n] -p 2 t ≤ p t (1 -p t ) = σ 2 CM C .
Thus, the APMC procedure enables us to reduce the variance of the estimator. The magnitude of the variance reduction is not easy to capture and depends on the situation. For illustration, consider a Pólya-Lundberg process where W = d Γ(2, 1) (see (ii) in Subsection 2.1) and a polynomial upper boundary given by

h β (t) = h(t) + β = t 2 + 3/2, t ≥ 0,
The probability of interest is p 2 = P (τ β > 2) = 0.568265. The variance of the APMC estimator is 0.184338 while the variance of the CMC estimator is 0.24534. So, the variance reduction when choosing the APMC procedure is arround 25%.

To explain a variance reduction, let us suppose that 5 trajectories are drawn for the process; they are displayed in Figure 4 (with the upper boundary). In general, the trajectories fall in three possible categories:

1. those ending below β for which both p N (t) = 1 τ β >t = 1 since there is no crossing, 2. those ending above h(t) + β for which both p N (t) = 1 τ β >t = 0 since crossing is a.s., 3. those ending between β and h(t) + β for which p N (t) ≤ 1 and 1 τ β >t = 0 or 1. For the 5 trajectories in Figure 4, the values of 1 τ β >t and p N (t) are given in Table 1. Intuitively, if many trajectories finish their run either below β or above h(t) + β, then p N (t) = 1 τ β >t very often and the variance reduction is not significant. On the contrary, when N (t) is between β and h(t) + β, then the time at which the jumps occur plays an important role and the APMC simulation decreases the variance significantly. In practice, if N (t) ∈ (β, h(t) + β), getting an approximation of p N (t) will imply repeated simulations of the jump times. This requires nested simulations that are known to be intense from a computational viewpoint. The variance reduction is expected to be larger when the likelihood for N (t) of being between β and h(t) + β is high. For the present example, the probability P [N (t) ∈ (β, h(t) + β)] is increasing in t, as shown in Figure 5.

Moreover, in Figure 6, the variance of the CMC and APMC estimators is represented as function of t, as well as the relative difference, expressed as a percentage,

∆σ 2 = σ 2 AP M C -σ 2 CM C σ 2 CM C
. We observe that for t < 1, the variances are of similar size because the two estimators are rather close. For t > 1, the decay of the variance of the APMC estimator becomes clear. So, a variance reduction of more than 30% is reached when t = 4. Table 2 gives some numerical values for p t , σ However, a direct link between the probability P [N (t) ∈ (β, h(t) + β)] and the variance reduction is not established so far. Changing the parameter values, the type of OSPP and the boundary shape allows us to appreciate the advantage of the APMC simulation. Further experiences can be conducted through the online accompaniment [START_REF] Goffard | Online accompaniment for "Boundary Crossing of Order Statistics Point Process[END_REF].

Let us add a remark on the use of formula (6.1) itself. This requires to determine the values of the Appell polynomials up to the degree h(t) + β , which is done by applying formula (2.13) with the recursion (2.14). These computations will be time consuming when h(t) + β is large. However, if it is unlikely that the OSPP reaches the h(t) + β , the APMC estimator provides a way to reduce the computational effort as the Appell polynomials will be computed up to the degree given by min { h(t) + β , N * (t)}, where N * (t) is the maximum of the OSPP at time t over all the replications.

Statistical perspectives. To close, it is worth to emphasize potential applications in statistics. Suppose that a sample of trajectories up to a given time horizon is available. Each trajectory can be summarized by the jump times and the number of jumps. As the APMC estimator is less volatile than the CMC estimator, it leads to a more reliable estimation. What makes the APMC estimator even more appealing is that we do not need to keep track of the jump times to evaluate it, only the number of jumps is necessary. This is quite convenient because it is not always possible in practice to spot and report exactly when does an event occurs. In such a situation of incomplete information upon the observed trajectories, we can still propose an estimation of the probability under study. The cost is, of course, the assumption that the observed phenomenon is governed by an OSPP for which we have to identify the type and the parameters. One could argue that the mixed Poisson process emcompasses a broad class of stochastic processes and is well suited in many modelling problems (see the monograph by Grandell [START_REF] Grandell | Mixed Poisson Process[END_REF]).

Definition 2 . 3 .

 23 Starting with G 0 (x|U ) = 1, the G n (x|U )'s satisfy the differential equations

Figure 1 :

 1 Figure 1: (blue, dashed) trajectory of the OSPP {N (t) ; t ≥ 0}; (black, plain) lower boundary.

(3. 6 )

 6 To obtain (3.5) and(3.6), it suffices to use in (3.1) the identity (2.16) for the A-G polynomials. For illustration, let us examine the special cases (i)-(iv) presented in section 2.1. Then, (3.1) and (2.16) give the following distributions for τ α . Corollary 3.3. (i) For a Poisson process,

Figure 2 :

 2 Figure 2: (blue, dashed) trajectory of the OSPP {N (t) ; t ≥ 0}; (black, plain) upper boundary.

Figure 3 :

 3 Figure 3: (blue, dashed) trajectory of the OSPP {N (t) ; t ≥ 0}; (black, plain) upper boundary; (red, plain) lower boundary.

Figure 4 :Table 1 :

 41 Figure 4: (blue, dashed) trajectories of the OSPP {N (t) ; t ≥ 0}; (red, plain) upper boundary. Trajectory 1 2 3 4 5 p N (t) 0.875 0 1 0 0.0821616 1 τ β >t 0 0 1 0 0Table 1: Values of 1 τ β >t and p N (t) for the 5 trajectories.

Figure 5 :

 5 Figure 5: The probability P [N (t) ∈ (β, h(t) + β)] as a function of t.

Figure 6 :

 6 Figure 6: In funcrion of t, on the left side the variances σ 2 CM C (dashed) and σ 2 AP M C (plain), on the right side the relative difference of variance ∆σ 2 .

Table 2 :

 2 2 AP M C , σ 2 CM C , and ∆σ 2 at different time horizons. Values of p t , σ 2 AP M C , σ 2 CM C , and ∆σ 2 at different time horizons.

	t	p t	σ 2 AP M C	σ 2 CM C	∆σ 2
	1.	0.62963 0.196159 0.233196 -15.8824
	2. 0.568265 0.184338 0.24534 -24.8643
	3. 0.562963 0.167682 0.246036 -31.8464
	4. 0.562619 0.158262 0.246079 -35.6867

Acknowledgement

We received support from the ARC project IAPAS of the Fédération Wallonie-Bruxelles.