On the Square Root of a Bell Matrix

In the context of Riordan arrays, the problem of determining the square root of a Bell matrix R=R(f(t)/t,f(t))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R={\mathcal {R}}(f(t)/t,\ f(t))$$\end{document} defined by a formal power series f(t)=∑k≥0fktk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(t)=\sum _{k \ge 0}f_kt^k$$\end{document} with f(0)=f0=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(0)=f_0=0$$\end{document} is presented. It is proved that if f′(0)=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f^\prime (0)=1$$\end{document} and f″(0)≠0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f^{\prime \prime }(0)\ne 0$$\end{document} then there exists another Bell matrix H=R(h(t)/t,h(t))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H={\mathcal {R}}(h(t)/t,\ h(t))$$\end{document} such that H∗H=R;\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H*H=R;$$\end{document} in particular, function h(t) is univocally determined by a symbolic computational method which in many situations allows to find the function in closed form. Moreover, it is shown that function h(t) is related to the solution of Schröder’s equation. We also compute a Riordan involution related to this kind of matrices.


Introduction
The concept of a (proper) Riordan array was introduced [18,19] as a generalization of the Pascal triangle. A Riordan array is defined as a pair of formal power series D = R(d(t), h(t)) with d(0) = 0, h(0) = 0, h (0) = 0. The usual way to represent the Riordan array R(d(t), h(t)) is by means of an infinite matrix (d n,k ), n, k ∈ N , its generic element being: The Pascal triangle is just the case: any sequence having s(t) = ∞ k=0 s k t k as its generating function, it is possible to prove that: n k=0 d n,k s k = [t n ]d(t)s(h(t)) (2) thus reducing the sum to the extraction of a coefficient from a formal power series.
Riordan arrays constitute a group with respect to the product defined by: which corresponds to the usual row-by-column product of two (infinite) matrices. The Riordan array I = R(1, t) acts as the identity and the inverse of D = R(d(t), h(t)) is the Riordan array: , h(t) (4) where h(t) is the compositional inverse of h(t), that is the power series such that h(h(t)) = h(h(t)) = t. For Pascal triangle we can easily find: Riordan arrays have another important characterization, due to Rogers [15]: every element d n+1,k+1 , n, k ∈ N , can be expressed as a linear combination of the elements in the preceding row, i.e.: d n+1,k+1 = a 0 d n,k + a 1 d n,k+1 + a 2 d n,k+2 + · · · = ∞ j=0 a j d n,k+j .
The sum is actually finite and the sequence A = (a k ) k∈N is called the Asequence of the Riordan array and we have h(t) = tA(h(t)), or A(t) = t h (t) . (6) For example, in the Pascal case we have A(t) = 1 + t. Sometimes, the Asequence is simpler than h(t) and in this case, going through this sequence can simplify the manipulation of the matrix. While computing the product D 1 * D 2 = D 3 between Riordan arrays, He and Sprugnoli [6] proved the following relation involving the corresponding A-sequences A 1 (t), A 2 (t) and A 3 (t) : There are several known subgroups of the Riordan group, in particular Bell subgroup contains elements of the form R(h(t)/t, h(t)); therefore, a Bell matrix is completely defined by its A-sequence. In this paper we investigate some questions about powers of Bell matrices. Let us start by considering the generalized Pascal triangle: By using the rules for the product between Riordan arrays we can easily compute: and by proceeding by mathematical induction we can prove the more general result: In classical applications, s ∈ N and P s r corresponds to a standard product among matrices.
However, let us suppose that rs = 1 with r, s ∈ Q; then in particular, is the square root of the Pascal triangle since P 1/2 * P 1/2 = P, as can be verified for the first 6 rows by performing the following product:  Now, let us consider the Bell matrix H = R(h(t)/t, h(t)) and consider the product of H by itself: therefore H can be seen as the square root of R(f (t)/t, f (t)).
Given h(t) we can compute f (t) : the question is, given f (t), is it possible to compute h(t) such that f (t) = h(h(t))? We point out that in these computations we are considering formal power series over R and that we need to have h(0) = f (0) = 0.
Obviously the problem can be generalized and one can ask the same question for the product of H by itself three times (or more) looking for a function h(t) such that g(t) = h(h(h(t))) for a given g(t): Studying fractional iteration of power series is a classical topic, addressed in several classical books such as [8,14,20] and the compositional square root is a particular case. However, as far as we know, this is the first attempt to study the problem in the context of Riordan arrays. Functions of Riordan arrays chunks, in particular the square root of Riordan arrays chunks, have been recently studied in [9] in terms of generalized Lagrange bases and Hermite polynomials. Here we are interested in studying the problem with a generating function approach and identifying the square root for the Riordan array in its entirety. In the next sections, we analyse the square root case and give some hints for general powers. Moreover, a result of Barry [1] is specialized to Bell matrices to find Riordan involutions associated with Bell matrices, another interesting problem related to powers of Riordan arrays.

The Square Root Case
We start by giving some example of square root matrices.

Self-composition of the Fibonacci Sequence
Let us consider the generating function of Fibonacci numbers, h(t) = t/(1 − t − t 2 ) = t + t 2 + 2t 3 + 3t 4 + 5t 5 + 8t 6 + · · · (sequence A000045 in OEIS), and the corresponding Riordan array: as can be verified on the first 6 rows of the involved matrices by computing the usual row-by-column product:

Self-composition of the Catalan Sequence
Let us consider the shifted generating function of Catalan numbers (sequence A000108 in OEIS), h(t) = (1− √ 1 − 4t)/2 = t+t 2 +2t 3 +5t 4 +14t 5 +42t 6 +· · · and the corresponding Riordan array: , matrix H is s the square root of the Riordan array as can be verified on the first 6 rows of the involved matrices by computing the usual row-by-column product:

The General Case
As already observed, given h(t) computing the self composition h(h(t)) is simple and proceeding in this way we can find square roots of many Bell matrices. The problem to be investigated is the converse, that is, given f (t) how to compute h(t) such that h(h(t)) = f (t). We study the problem by a Riordan array approach, as follows. We consider the Riordan array S = R(1, h(t)) and denote s n,k = [t n ]h(t) k ; the first 6 rows are which corresponds to a matrix-vector product, as illustrated by the first 6 involved rows: The matrix-vector product translates into the following system of equations: For simplifying, let us consider the quite common case in combinatorics f 1 = 1 and f 2 = 0, then, from the first equation we have h 1 = ±1 but from the previous observation and the additional condition on f 2 we need to choose the + sign and therefore we have: as expected and illustrated in Sect. Since, in general, a matrix can have several square roots (see, e.g., [7]), the previous result means that a Bell matrix has a unique square root which is a Bell matrix that can be computed by the computational method shown before; this is stated in the following theorem: We observe that different algorithms for computing the coefficients in the expansion of the rth root of a formal power series f (t) under composition can be found in [2].
Next questions is: is it possible to find a closed form expression h(t) from f (t)? For example, given Before we start discussing this problem, let's see some features of Bell matrices and their square root.

Theorem 2.2. Let F = R(f (t)/t, f (t)) be a Bell matrix with f (0) = 0, f (0) = 1 and f (0) = 0 and H = R(h(t)/t, h(t))
with h(0) = 0, h (0) = 1 and h (0) = 0 its square root; moreover let A f and A h be the A-sequences of F and H. Then, the following relations hold: For example, for the Pascal case f (t) = t/(1 − t) and A(t) = 1 + t, thus we obtain the following relations: However, these relations are equivalent and they do not allow to find h(t).

Guessing h(t)
In some situations, by looking at the very first coefficients of the function h(t) defined in Theorem 2.1, it is possible to guess the generating function, for example consulting the OEIS [12]. This is surely true for the examples just discussed, since they correspond to well-known sequences. Sometimes this is possible also for less well-known sequences. As an example, let us consider the generating function By using the relations between the coefficients of f (t) and h(t) we find the first values: but looking at them it is not easy to guess the corresponding generating function. However, we can make an hypotheses and imagine that it is a rational function, that is, sequence (h n ) n∈N satisfies a linear recurrence relation with

Finally, by computing h(h(t)) we can verify that it corresponds to f (t). This proves that
is the square root of the Riordan array: thus generalizing the Fibonacci case.
If the check f (t) = h(h(t)) had not given the desired result, we could have tried a higher order recurrence, which exists if the function h(t) is rational and is computable with the method of coefficient already mentioned. There is the classical GFUN Maple package [16] that allow you to guess the generating function with a similar approach, even in the algebraic case, that is, when the linear recurrence relation has polynomial coefficients. For example, giving the list 0, 1, 1, 2, 5, 14, 42 as input to the guessgf function of the package, will return the generating function of Catalan numbers 1− . As a more complex example, if we consider the series and give it as input to the seriestoalgeq function of the package, we obtain the equation:  h(h(t)) and applying the computational method illustrated in the previous section. As a final example, let us consider again the Catalan triangle, which is one of the most studied Riordan array due to the many combinatorial interpretations of the involved numbers, and look for its square root. The computational method gives for the first 6 rows the following decomposition: However, we did not succeed in finding the closed form or an equation for or, equivalently, for apart from those corresponding to Theorem 2.2 which correspond to A f (t) = 1 1−t . We will return to this example in the next section.

A Solution Through Schröder's Equation
The computation of h(t) in Theorem 2.1 seems to be connected to the solution of Schröder's equation [17], a functional equation with one independent variable: given the function f the aim is to find the function Ψ such that Upon iteration of the functional equation we have: and where n is not necessarily positive or integer. In particular, for the functional square root we have Hence, for computing f 1/2 (t) we need to find the Schröder's auxiliary function Ψ(t), and construct the inverse function Ψ(t). For example, as illustrated in [5], if f (t) = 2t(1 + t) then Schröder's equation is solved by Ψ(t) = 1 2 ln(1 + 2t) and Ψ(t) = 1 2 (e 2t − 1) which gives f 1/2 (t) = 1 2 ((1 + 2t) √ 2 − 1). However, as observed in [5], it is crucial to note that in the limit s → 1, all iterates and inverse iterates lose their distinction and degenerate to the identity map, f 0 (t) = t, and the method fails. For this reason, if f (0) = 1, as in our case, one augments f 1 (t) in Schröder's equation to sf 1 (t), and takes the s → 1 limit only at the very end of the calculations.
Let us show this approach for the Pascal case, by considering f (t) = st 1−t . We need to find a function Ψ(t) with Ψ(0) = 0 and Ψ (0) = 0 such that for some constant s = 1. By assuming that Ψ(t) = k≥1 Ψ k t k we have: The coefficients in the development of the previous expression should be equal to zero, and since Hence, by extracting the (n + 1)th coefficient we have: or, By varying n, we find: and, more generally, Ψ n =
Finally, the solution corresponding to s = 1 can be found by computing the following limit: f 1/r (t) = Ψ(s 1/r Ψ(t)) = s 1/r (1 − s)t 1 − s + s 1/r t − t and the s → 1 limit gives t 1− t r , as expected.
By using this approach, in principle, one could compute h(t) = f 1/2 (t) and find the square root of a Bell matrix explicitly. However, solving the Schröder's equation for a generic function f (t) is not simple: the crucial point, in fact, is to determine the function Ψ.
As another example, let us consider f (t) = s(1− √ 1−4t) 2 and the corresponding Schröder's equation: We proceed similarly to the previous example and begin by performing the following extraction, related to the generating function C(t) of Catalan numbers: where 1 + w = C(t) and we used Lagrange Inversion Formula (see, e.g, [10]). Therefore we have: Hence, by extracting the nth coefficient, we have: We can state this result in the following theorem: