On the Total Positivity and Accurate Computations of r -Bell Polynomial Bases

: A new class of matrices deﬁned in terms of r -Stirling numbers is introduced. These r - Stirling matrices are totally positive and determine the linear transformation between monomial and r -Bell polynomial bases. An efﬁcient algorithm for the computation to high relative accuracy of the bidiagonal factorization of r -Stirling matrices is provided and used to achieve computations to high relative accuracy for the resolution of relevant algebraic problems with collocation, Wronskian, and Gramian matrices of r -Bell bases. The numerical experimentation conﬁrms the accuracy of the proposed procedure.


Introduction
The Stirling numbers form combinatorial sequences with remarkable properties and many applications in Combinatorics, Number Theory, or Special Functions, among others.The Stirling numbers were so named by Nielsen (cf.[1]) in honor of James Stirling, who computed them in 1730 (cf.[2]).Interested readers are referred to [3][4][5] to find a nice introduction and many nice properties of Stirling numbers.In the literature, many kinds of generalizations of Stirling numbers have been introduced and studied ( [3,6]).In this paper, we shall focus on the r-Stirling numbers whose combinatorial and algebraic properties, in most cases, generalize those of the regular Stirling numbers.
The r-Bell polynomials (see [7]) are considered as an efficient mathematical tool for combinatorial analysis (cf.[5]) and can be applied in many different contexts: in the evaluation of some integrals and alternating sums [8,9], for the analysis of the internal relations for the orthogonal invariants of a positive compact operator [10], in the Blissard problem [5], for the Newton sum rules for the zeros of polynomials [11], in the recurrence relations for a class of Freud-type polynomials [12] and in many other subjects.
In this paper, a class of r-Stirling matrices will be introduced.It will be shown that these matrices determine the linear transformations between monomial and r-Bell polynomial bases and have the property of being totally positive and, therefore, having all of their minors as non-negative.Several applications of the class of totally positive matrices can be found in [13][14][15].Totally positive matrices can be expressed as a particular product of bidiagonal matrices (see [16,17]).This bidiagonal decomposition provides a representation of the totally positive matrices that exploits their total positivity property and allows us to derive algorithms to high relative accuracy for the resolution of relevant algebraic problems, such as, the computation of their inverses, eigenvalues, singular values, or the solution of some systems of linear equations.
A major issue in Numerical Linear Algebra is to obtain algorithms to high relative accuracy, because the relative errors in their computations will have the order of the machine precision and will not be drastically affected by the dimension or conditioning of the considered matrices.Consequently, the design of algorithms adapted to the structure of totally positive matrices and achieving computations to high relative accuracy has attracted the interest of many researchers (see [18][19][20][21][22][23][24]).
Let us recall that many interpolation, numerical quadrature, and approximation problems require algebra computations with collocation matrices of the considered bases.When solving Taylor interpolation problems, Wronskian matrices have to be considered.On the other hand, Gramian matrices define linear transformations to obtain orthogonal from nonorthogonal bases.Moreover, the inversion of Gramian matrices is also required when approximating, in the least-squares sense, curves by linear combinations of control points and the basis functions.The large range of applications of r-Bell polynomials has motivated us to exploit the total positivity property of r-Stirling matrices to design fast and accurate algorithms for the resolution of algebraic problems with collocation, Wronskian, and Gramian matrices of r-Bell polynomial bases.
The layout of this paper is as follows.Section 2 summarizes some notations and auxiliary results.Section 3 focuses on the class of r-Stirling matrices.Their bidiagonal decomposition is provided and their total positivity property is deduced.In Section 4, r-Bell polynomials are considered.Collocation, Wronskian, and Gramian matrices of r-Bell polynomial bases are shown to be totally positive.Then, using the proposed factorization for r-Stirling matrices, the resolution of the above mentioned algebraic problems can be achieved to high relative accuracy.Finally, Section 5 presents numerical experiments confirming the accuracy of the proposed methods for the computation of eigenvalues, singular values, inverses, or the solution of some linear systems related to collocation, Wronskian, and Gramian matrices of r-Bell polynomial bases.

Notations and Auxiliary Results
Let us recall that a matrix is totally positive (respectively, strictly totally positive) if all its minors are nonnegative (respectively, positive).
The Neville elimination is an alternative procedure to Gaussian elimination (see [16,17,25]).Given a nonsingular matrix A = (a i,j ) 1≤i,j≤n+1 , its Neville elimination calculates a sequence of (n + 1) × (n + 1) matrices A (k) , k = 1, . . ., n + 1, so that their elements below the main diagonal of the k − 1 first columns are zeros and then A (n+1) is upper triangular. From ) 1≤i,j≤n+1 as follows, with A (1) := A. The element is called the (i, j) pivot and we say that p i,i is the i-th diagonal pivot of the Neville elimination of A. Whenever all pivots are nonzero, no row exchanges are needed in the Neville elimination procedure.The value for 1 ≤ j < i ≤ n + 1, is the (i, j) multiplier of the Neville elimination of A.
By Theorem 4.2 and the arguments of p. 116 of [17], a nonsingular totally positive A ∈ R (n+1)×(n+1) can be written as follows, where n+1) and G i ∈ R (n+1)×(n+1) , i = 1, . . ., n, are the totally positive, lower and upper triangular bidiagonal matrices with the unit diagonal described by and D ∈ R (n+1)×(n+1) is a diagonal matrix with positive diagonal entries.In fact, the elements in the diagonal of D are the diagonal pivots of the Neville elimination of A. On the other hand, the off-diagonal elements m i,j and m i,j are the multipliers of the Neville elimination of A and A T , respectively.Using the results in [16,17,25] and Theorem 2.2 of [22], the Neville elimination of A can be used to derive the following bidiagonal factorization (5) of its inverse matrix, where F i ∈ R (n+1)×(n+1) and G i ∈ R (n+1)×(n+1) are the lower and upper triangular bidiagonal matrices with the form described in (6), which are obtained by replacing the off-diagonal entries {m i+1,1 , . . ., m n+1,n+1−i } and { m i+1,1 , . . ., m n+1,n+1−i } by {−m i+1,i , . . ., −m n+1,i } and {− m i+1,i , . . ., − m n+1,i }, respectively.The structure of the bidiagonal matrix factors in ( 7) is described as follows, Let us also observe that, if a matrix A ∈ R (n+1)×(n+1) is nonsingular and totally positive, its transpose A T is also nonsingular and totally positive, and where Fi = G T i , Gi = F T i and F i , G i i = 1, . . ., n, are the lower and upper triangular bidiagonal matrices in (5).
The total positivity property of a given matrix can be characterized by analyzing the sign of the diagonal pivots and multipliers of its Neville elimination, as shown in Theorem 4.1, Corollary 5.5 of [25] and the arguments of p. 116 of [17].
Theorem 1.A nonsingular matrix A is totally positive (respectively, strictly totally positive) if and only if the Neville elimination of A and A T can be performed without row exchanges, all the multipliers of the Neville elimination of A and A T are nonnegative (respectively, positive) and the diagonal pivots of the Neville elimination of A are all positive.
In [19], the bidiagonal factorization (5) of a nonsingular and totally positive A ∈ R (n+1)×(n+1) is represented by defining a matrix BD(A) = (BD(A) i,j ) 1≤i,j≤n+1 such that This representation will allow us to define algorithms adapted to the totally positive structure and provide accurate computations with the matrix.
Let us observe that, for the matrix A in (4), the bidiagonal decomposition (5) can be written as follows: Let us also recall that a real value x is computed to high relative accuracy whenever the computed value x satisfies where u is the unit round-off and K > 0 is a constant, which is independent of the arithmetic precision.High relative accuracy implies a great accuracy in the computations, since the relative errors have the same order of the machine precision and the accuracy is not affected by the dimension or the conditioning of the problem we are solving.
A sufficient condition to assure that an algorithm can be computed to high relative accuracy is the non-inaccurate cancellation condition, sometimes denoted as the NIC condition, which is satisfied if the algorithm does not require inaccurate subtractions and only evaluates products, quotients, sums of numbers of the same sign, subtractions of numbers of the opposite sign, or subtraction of initial data (cf.[18,19]).
If the bidiagonal factorization (5) of a nonsingular and totally positive matrix A can be computed to high relative accuracy, the computation of its eigenvalues and singular values, the computation of A −1 , and even the resolution of systems of linear equations Ax = b, for vectors b with alternating signs, can be also computed to high relative accuracy, using the algorithms provided in [26].
Let (u 0 , . . ., u n ) be a basis of a space U of functions defined on I ⊆ R. Given a sequence of parameters x 1 < . . .< x n+1 on I, the corresponding collocation matrix is defined by The system (u 0 , . . ., u n ) of functions defined on I ⊆ R is totally positive if all of its collocation matrices M(x 1 , . . ., x n+1 ) are totally positive.
If the space U is formed by n-times continuously differentiable functions and x ∈ I, the Wronskian matrix at x is defined by: where u (i) (x), i ≤ n, denotes the i-th derivative of u at x. Now, let us suppose that U is a Hilbert space of functions under a given inner product u, v defined for any u, v ∈ U.Then, given linearly independent functions u 0 , . . ., u n in U, the corresponding Gram matrix is the symmetric matrix defined by:

Bidiagonal Factorization of r-Stirling Matrices
Let us recall that given n, m ∈ N with m < n, the (signless) Stirling number of the first kind, usually denoted by n m , can be defined combinatorially as the number of permutations of the set {1, . . ., n} having m cycles.
On the other hand, the Stirling number of the second kind, usually denoted by n m , coincides with the number of partitions of the set {1, . . ., n} into m non-empty disjoint blocks (cf.[3]).Moreover, the Bell number B n gives the total number of partitions of the set {1, . . ., n}; that is, The first and second kinds of Stirling numbers can be recursively computed as follows, with the initial conditions The r-Stirling numbers are polynomials in r generalizing the regular Stirling numbers and they count some restricted permutations and partitions of sets with a given number of elements.
For r ∈ N, the first kind r-Stirling number, denoted by n m r , counts the number of permutations of the elements {1, . . ., n} having m cycles, such that 1, 2, . . ., r are in distinct cycles.In addition, the second kind r-Stirling number, denoted by n m r , is defined as the number of set partitions of {1, . . ., n + r} into m + r blocks such that the first r elements are in distinct blocks.Thus the r-Bell number B r n gives the total number of partitions of the set {1, . . ., n + r}, such that the first r elements are in distinct blocks; that is, Clearly, for r = 0, the r-Stirling numbers coincide with the corresponding regular Stirling numbers.Theorems 1 and 2 of [3] prove that the r-Stirling numbers satisfy the recurrence relations verified by the regular Stirling numbers (see ( 14)); that is, but with different initial conditions, Moreover, for some particular values, the corresponding r-Stirling numbers can be easily obtained: Finally, let us recall a "cross" recurrence relating r-Stirling numbers of the second kind with a different r that is going to be used in this paper (see Theorem 4 of [3]).The r-Stirling numbers of the second kind satisfy Now, we are going to define triangular matrices whose entries are given in terms of r-Stirling numbers.We shall analyze their Neville elimination and derive their bidiagonal factorization (5).Definition 1.For r, n ∈ N with 0 ≤ r ≤ n, the (n + 1) × (n + 1) the first kind of r-Stirling matrix is the lower triangular matrix L r n = (l r i,j ) 1≤i,j≤n+1 with l r i,j := On the other hand, the (n + 1) × (n + 1) second kind of r-Stirling matrix, is the lower triangular matrix Lr n = ( lr i,j ) 1≤i,j≤n+1 with The following results provide the expression of the pivots and multipliers of the Neville elimination of r-Stirling matrices and their inverses.Then, their bidiagonal factorization (5) will be deduced and their total positivity property analyzed.
Theorem 2. For r, n ∈ N with r ≤ n, the first kind of r-Stirling matrix L r n ∈ R (n+1)×(n+1) is totally positive and admits the following factorization where F i , i = 1, . . ., n, are the lower triangular bidiagonal matrices of the form (6), whose offdiagonal entries are given by Moreover, (L r n where F i , i = 1, . . ., n, are the lower triangular bidiagonal matrices of the form (6), whose offdiagonal entries are given by m i,j = −(r Proof.Let us define L (1) := L r n and . ., n, the matrix obtained after k − 1 steps of the Neville elimination of L r n .We shall deduce by induction on k ∈ {1, . . ., n + 1} that For k = 1, identities (28) follow from (22).Now, suppose that ( 28) is verified for some k ∈ {1, . . ., n}.Then, taking into account (19), we can write By (1), we have l l−1,j and, using ( 16), (28), and (29), we derive corresponding to the identity (35) for k + 1.By ( 2), (19), and (28), we deduce that the pivots of the Neville elimination of L r n satisfy and so the diagonal pivots are p i,i = r r r = 1 for i = 1, . . ., n + 1 (see (17)).Moreover, the following expression for the multipliers can be deduced Clearly, p i,i > 0, for i = 1, . . ., n + 1, and m i,j > 0, for 1 ≤ j < i ≤ n + 1.Then, using Theorem 1, we deduce that L r n is totally positive.Finally, taking into account the bidiagonal factorization (7) for the inverse of a nonsingular totally positive matrix, the factorization ( 26) for (L r n ) −1 can be deduced easily.
The provided bidiagonal factorization (5) of the matrix Lr n can be stored in a compact form through BD( Lr elsewhere. (37) Remark 1.Given a lower triangular (n + 1) × (n + 1) matrix, for any k ≤ n + 1, the determinants of the submatrices using row i 1 , . . ., i k and columns j 1 , . . ., j k with i h ≥ j h for all h ≤ k are called nontrivial minors because all other minors are zero.A lower triangular matrix is called ∆STP if all its nontrivial minors are positive.Since Theorem 2 (respectively, Theorem 3) proves that the lower triangular matrix with unit diagonal L r n (respectively, Lr n ) has all multipliers of its Neville elimination positive, by Theorem 4.3 of [16], we conclude that L r n (respectiveley, Lr n ) is in fact ∆STP.
Now, taking into account Theorems 2 and 3, we provide the pseudocode of Algorithms 1 and 2, respectively.Specifically, Algorithm 1 computes the matrix form BD(L r n ) in (30), for the bidiagonal decomposition of the first kind r-Stirling matrix L r n in (24).Moreover, Algorithm 2 computes the matrix form BD( Lr n ) in (37) for the bidiagonal decomposition of the r-Stirling matrix of the second kind Lr n in (31).The computational cost of both algorithms is O(n 2 ).Let us observe that none of the algorithms requires inaccurate subtractions, and so the provided matrices are computed to high relative accuracy.In fact, all of their entries are natural numbers.
Algorithm 1: Computation to high relative accuracy of BD(L r n ) (see (30)) Require: n, r Ensure: BDL r n bidiagonal decomposition of L r n to high relative accuracy BDL r n = eye(n + 1)

Accurate Computations with r-Bell Bases
Let P n (I) be the (n+1)-dimensional linear space formed by all polynomials in the variable x defined on a real interval I and whose degree is not greater than n; that is, For n ∈ N, the n degree r-Bell polynomial of the first kind is defined as Then we can define a system (B r 0 , . . ., B r n ) of r-Bell polynomials of the first kind that satisfy where m i (x) := x i , i = 0, . . ., n and L r n is the (n + 1) × (n + 1) first kind of r-Stirling matrix defined by (22).By (24), det L r n = 1 and we can guarantee that (B r 0 , . . ., B r n ) T is a basis of P n (R).
Definition 2. We say that (B r 0 , . . ., B r n ) is the r-Bell basis of the first kind of the polynomial space P n (R).
On the other hand, the r-Bell polynomials of the second kind { B r n } n≥0 in [3,27] are called r-Bell polynomials, which are defined by their generating function They can also be written in terms of the monomial basis as follows, and then ( Br 0 , . . ., B r n where Lr n is the (n + 1) × (n + 1) second kind of r-Stirling matrix defined by (23).By (31), det Lr n = 1 and we can guarantee that ( Br 0 , . . ., Br n ) T is a basis of P n (R).
Definition 3. We say that ( Br 0 , . . ., B r n ) is the r-Bell basis of the second kind of P n (R).
Let us recall that the monomial basis (m 0 , . . ., m n ), with m i (x) = x i for i = 0, . . ., n, is a strictly totally positive basis of P n (0, +∞), and so the collocation matrix is strictly totally positive for any increasing sequence of positive values 0 < x 1 < . . .< x n+1 (see Section 3 of [19]).In fact, V is the Vandermonde matrix at the considered nodes.Let us recall that Vandermonde matrices have relevant applications in linear interpolation and numerical quadrature (see for example [28,29]).As for BD(V), we have and it can be easily checked that the computation of BD(V) does not require inaccurate cancellations and can be performed to high relative accuracy.Taking into account the total positivity of the Vandermonde matrices and the total positivity of r-Stirling matrices of the first and second kinds, we shall derive the total positivity of r-Bell bases, as well as factorizations providing computations to high relative accuracy when considering their collocation matrices.
Theorem 4. The r-Bell basis of first and second kind are totally positive bases of P n (0, +∞).Moreover, given 0 < x 1 < . . .< x n+1 , the collocation matrices and their bidiagonal factorization (5) can be computed to high relative accuracy.
Proof.Given 0 < x 1 < . . .< x n+1 , by formulae (39) and ( 41), the collocation matrices in (44) satisfy where V is the Vandermonde matrix (42), and L r n and Lr n are r-Stirling matrices of the first and second kind, respectively.
It is well known that V is strictly totally positive at 0 < x 1 < . . .< x n+1 and its decomposition (5) can be computed to high relative accuracy (see [19]).By Theorem 2 (respectively, Theorem 3), L r n (respectively, Lr n ) is a nonsingular totally positive matrix and its decomposition (5) can be computed to high relative accuracy.Taking into account these facts, (L r n ) T (respectively, ( Lr n ) T ) is a nonsingular totally positive matrix and its bidiagonal decomposition can be also computed to high relative accuracy (see (10)).Therefore, we can deduce that B (respectively, B) is a totally positive matrix since it is the product of totally positive matrices (see Theorem 3.1 of [13]).Moreover, using Algorithm 5.1 of [26], if the decomposition (5) of two nonsingular totally positive matrices is provided to high relative accuracy, then the decomposition of the product can be obtained to high relative accuracy.Consequently, B (respectively, B) and its decomposition (5) can be obtained to high relative accuracy.
Corollary 1 of [20] provides the factorization (5) of W := W(m 0 , . . ., m n )(x), the Wronskian matrix of the monomial basis (m 0 , . . ., m n ) at x ∈ R. For the matrix representation BD(W), we have Taking into account the sign of the entries of BD(W) and Theorem 1, one can derive that the Wronskian matrix of the monomial basis is totally positive for any x > 0.Moreover, the computation of (46) satisfies the NIC condition and W can be computed to high relative accuracy.Using (46), computations to high relative accuracy when solving algebraic problems related to W have been achieved in [20] for x > 0.
Using formula (39), it can be checked that So, with the reasoning of the proof of Theorem 4, the next result follows and we can also guarantee computations to high relative accuracy when solving algebraic problems related to Wronkian matrices of r-Bell polynomial bases at positive values.
Theorem 5. Let (B r 0 , . . ., B r n ) and ( Br 0 , . . ., Br n ) be the r-Bell bases of the first and second kind, respectively.For any x > 0, the Wronskian matrices W r n := W(B r 0 , . . ., B r n )(x) and Wr n := Let us recall that if BD(A) can be obtained to high relative accuracy, then the MATLAB functions TNEigenValues, TNSingularValues, TNInverseExpand and TNSolve available in the software library TNTools in [30], take as input argument BD(A), and compute to high relative accuracy the eigenvalues and singular values of A, the inverse matrix A −1 (using the algorithm presented in [22]) and even the solution of linear systems Ax = b, for vectors b with alternating signs.
In order to check the accuracy of our algorithms, we have performed several matrix computations using the mentioned routines available in [30], with the matrix form (11) of the bidiagonal factorization (5) as an input argument.The obtained approximations have been compared with the respective approximations obtained by traditional methods provided in Matlab R2022b.In this context, the values provided by Wolfram Mathematica 13.1 with 100-digit arithmetic have been taken as the exact solution of the considered algebraic problem.For the sake of brevity, only a few of these experiments will be shown.
The relative error of each approximation has also been computed in Mathematica with 100-digit arithmetic as e := |y − y|/|y|, where y denotes the exact solution and y the computed approximation.
Computation of eigenvalues and singular values.For all considered matrices, we have compared the smallest eigenvalue and singular value obtained using the proposed bidiagonal decompositions provided by Algorithms 3-5 with the functions TNEigenvalues and TNSingularvalues, and the smallest eigenvalue and singular value computed with the Matlab commands eig and svd, respectively.Note that the eigenvalues of the Wronskian matrices are exact.Furthermore, the singular values of the Gramian matrices considered coincide with their eigenvalues (since Gramian matrices are symmetric).
The relative errors are shown in Figure 2. Note that our approach accurately computes the smallest eigenvalue and singular value regardless of the 2-norm condition number of the considered matrices.In contrast, the Matlab commands eig and svd return results that are not accurate at all.
Computation of inverses.Further to this, for all considered matrices, we have compared the inverse obtained using the proposed bidiagonal decompositions provided by Algorithms 3-5 with the function TNInverseExpand and the inverse computed with the Matlab command inv.As shown in Figure 3, our procedure provides very accurate results.On the contrary, the results obtained with Matlab reflect poor accuracy.
Resolution of linear systems.

Conclusions
We have introduced a new class of matrices called r-Stirling matrices.These matrices are defined in terms of r-Stirling numbers and have interesting properties, such as being totally positive.They also play a significant role in determining the linear transformation between monomial and r-Bell polynomial bases.
The paper further discusses an efficient algorithm for computing the bidiagonal factorization of r-Stirling matrices to high relative accuracy.This factorization is used for solving relevant algebraic problems involving collocation, Wronskian, and Gramian matrices associated with r-Bell polynomial bases.To validate the proposed procedure and its usefulness, the paper presents numerical experiments that confirm its accuracy.
the diagonal pivots and multipliers of the Neville elimination of A and A T are all positive, we conclude that A is strictly totally positive.The above factorization can be represented in the following matrix form:

Figure 1 .
Figure 1.The 2-norm condition number of collocation, Wronskian, and Gramian matrices of r-Bell bases of the first and second kind.
Finally, for all considered matrices, we have compared the solution of the linear systems B r n c = d, Br n c = d, W r n c = d, Wr n c = d, G r n c = d and Ḡr n c = d where d = ((−1) i+1 d i ) 1≤i≤n+1 and d i , i = 1, . . ., n + 1, are random nonnegative integer values, obtained using the proposals for bidiagonal decompositions provided by Algorithms 3-5 with the function TNSolve and the solutions obtained within the Matlab command \.As opposed to the results obtained with the command \, the proposed procedure preserves the accuracy for all of the dimensions which have been taken into account.Figure 4 illustrates the relative errors.

Figure 3 .
Figure 3. Relative error of the approximations to the inverse of W r n with r = 5 and t = 60 and Ḡrn