This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/BellTransform
The Bell Transform
Summary: We introduce the partial Bell polynomials in two flavors: in the traditional way based on the variables x_{1}, x_{2}, ... and based on the variables x_{0}, x_{1}, x_{2}, .... The coefficients of the x_{0}based univariate polynomials give the new and unexplored A257563.
Next we base the computation of the Bell transform on partitions which gives a surprising fast algorithm for small n which is similar to the algorithm used in SageMath for Bell polynomials.
Iterating the Bell transform using the row sums of the triangles as input for the next level gives us the higher order Bell numbers. In particular we recover Joerg Arndt's award winning A187761 as the secondorder Bell numbers.
We then apply the Bell transform and its inverse to the power functions, the multifactorials, the rising factorials and the falling factorials showing more than 50 examples of sequences in the database. We found that almost none of them was described as the result of a Bell transformation.
Contents
 1 ♦ Bell polynomials
 2 ♦ The Bell transform based on partitions
 3 ♦ The Bell matrix
 4 ♦ The inverse Bell transform
 5 ♦ The Bell inverse of a sequence
 6 ♦ The coefficients of the Bell polynomials
 7 ♦ Bell numbers of higher order
 8 ♦ Associated Stirling subset numbers
 9 ♦ The Bell transform of multifactorials
 10 ♦ The Bell transform of powers
 11 ♦ The Bell transform of rising factorials
 12 ♦ The Bell transform of falling factorials
 13 ♦ The Bell transform of monotonic factorials
 14 ♦ An implementation for Maple
 15 ♦ Index of triangles generated by the Bell transform
 16 ♦ References
♦ Bell polynomials
The Bell transform is a transformation which maps a sequence
a = a_0, a_1, a_2,...
to a triangle
T(0,0) T(1,0) T(1,1) T(2,0) T(2,1) T(2,2) T(3,0) T(3,1) T(3,2) T(3,3)
defined by
T(n,k) = Sum_{j=0..nk+1} binomial(n1,j1)*a_j*T(nj,k1) if k != 0 else T(n, k) = (a_0)^n.
With Sage this can be implemented as
def bell_trans(n,k,a): @cached_function def T(n, k): if k==0: return a[0]^n return sum(binomial(n1,j1)*a[j]*T(nj,k1) for j in (0..nk+1)) return T(n,k)
Let us look at the first few rows of the generated triangle.
(*) a_0^0 a_0^1, a_1 a_0^2, a_1*a_0+a_2, a_1^2 a_0^3, a_1*a_0^2+2*a_2*a_0+a_3, a_1^2*a_0+3*a_2*a_1, a_1^3 a_0^4, a_1*a_0^3+3*a_2*a_0^2+3*a_3*a_0+a_4, a_1^2*a_0^2+5*a_1*a_2*a_0 +4*a_3*a_1+3*a_2^2, a_1^3*a_0+6*a_2*a_1^2, a_1^4
For example to compute the unsigned Stirling numbers of the first kind A132393 one can write
s = [0]+[factorial(i) for i in (0..5)] for n in (0..5): [bell_trans(n,k,s) for k in (0..n)] [1] [0, 1] [0, 1, 1] [0, 2, 3, 1] [0, 6, 11, 6, 1] [0, 24, 50, 35, 10, 1]
A second example shows that there are other important numbers which can be generated this way which do not always reduce to integers.
N = 6 t = [2^(1j) if is_odd(j) else 0 for j in (0..N)] for n in (0..N): [bell_trans(n,k,t) for k in (0..n)] [1] [0, 1] [0, 0, 1] [0, 1/4, 0, 1] [0, 0, 1, 0, 1] [0, 1/16, 0, 5/2, 0, 1] [0, 0, 1, 0, 5, 0, 1]
These are the coefficients of the central factorials [1]
Another way to look at the symbolic sums in the triangle (*) above is to interpret the a_{n} as variables x_{n} of a multivariate polynomial. Then one arrives at the partial Bell polynomials:
def partial_bell_polynomial(n,k): v = var(['x_'+str(i) for i in (0..n+1)]) return bell_trans(n,k,v).expand() for n in (0..4): [partial_bell_polynomial(n,k) for k in (0..n)] [1] [x_0, x_1] [x_0^2, x_0*x_1 + x_2, x_1^2] [x_0^3, x_0^2*x_1 + 2*x_0*x_2 + x_3, x_0*x_1^2 + 3*x_1*x_2, x_1^3] [x_0^4, x_0^3*x_1 + 3*x_0^2*x_2 + 3*x_0*x_3 + x_4, x_0^2*x_1^2 + 5*x_0*x_1*x_2 + 3*x_2^2 + 4*x_1*x_3, x_0*x_1^3 + 6*x_1^2*x_2, x_1^4]
Note that these polynomials depend on the variables x_{0}, x_{1}, x_{2}, ... For some but probably not for a good reason the partial Bell polynomials are often defined to depend only on the variables x_{1}, x_{2}, ... In this case the first column of the above triangle is just ignored (except for the case (0, 0) which is added somewhat unmotivated and unsystematic). To get this in the older literature common form (see for example L. Comtet, Advanced Combinatorics, page 135) from our form is easy: we just have to set x_{0} = 0.
With Sage we can write:
for n in (0..4): [partial_bell_polynomial(n,k).subs(x_0=0) for k in (0..n)] [1] [0, x_1] [0, x_2, x_1^2] [0, x_3, 3*x_1*x_2, x_1^3] [0, x_4, 3*x_2^2 + 4*x_1*x_3, 6*x_1^2*x_2, x_1^4]
So far we have only looked at the partial Bell polynomials. The Bell polynomials are defined as
B_n = Sum_{k=0..n} B_{n,k}. def bell_polynomial(n): return sum(partial_bell_polynomial(n,k) for k in (0..n)) for n in (0..3): bell_polynomial(n) 1 x_0 + x_1 x_0^2 + x_0*x_1 + x_1^2 + x_2 x_0^3 + x_0^2*x_1 + x_0*x_1^2 + x_1^3 + 2*x_0*x_2 + 3*x_1*x_2 + x_3
Two things remain to be mentioned:
First, it is sometimes convenient to reduce the infinite number of variables x_0, x_1, x_2, ... to one variable x only.
Second, in the considerations above we worked with 'symbolic sums' (or in the Sage parlance in the Symbolic Ring SR). But of course we can also work with 'real' polynomials.
Putting these two requirements together we arrive at polynomials in the univariate polynomial ring in x over the rational field.
def univariate_bell_polynomial(n): p = bell_polynomial(n).subs(x_0=0) q = p({p.variables()[i]:x for i in range(len(p.variables()))}) R = PolynomialRing(QQ, 'x') return R(q) for n in (0..6): univariate_bell_polynomial(n) 1 x x^2 + x x^3 + 3*x^2 + x x^4 + 6*x^3 + 7*x^2 + x x^5 + 10*x^4 + 25*x^3 + 15*x^2 + x x^6 + 15*x^5 + 65*x^4 + 90*x^3 + 31*x^2 + x
If we extract the coefficients of these polynomials we find:
for n in (0..6): univariate_bell_polynomial(n).list() [1] [0, 1] [0, 1, 1] [0, 1, 3, 1] [0, 1, 7, 6, 1] [0, 1, 15, 25, 10, 1] [0, 1, 31, 90, 65, 15, 1]
Lo and behold, these are the Stirling subset numbers A048993! This must have been observed before ;)
Note the first substitution subs(x_0=0) which indicates the 'classical case' in the sense discussed above. If instead we use the uniform substitution subs(x_n=x) for all n we get:
def x0_based_univariate_bell_polynomial(n): p = bell_polynomial(n) q = p({p.variables()[i]:x for i in range(len(p.variables()))}) R = PolynomialRing(QQ,'x') return R(q) for n in (0..6): x0_based_univariate_bell_polynomial(n) 1 2*x 3*x^2 + x 4*x^3 + 5*x^2 + x 5*x^4 + 14*x^3 + 10*x^2 + x 6*x^5 + 30*x^4 + 48*x^3 + 19*x^2 + x 7*x^6 + 55*x^5 + 158*x^4 + 149*x^3 + 36*x^2 + x for n in (0..6): x0_based_univariate_bell_polynomial(n).list() [1] [0, 2] [0, 1, 3] [0, 1, 5, 4] [0, 1, 10, 14, 5] [0, 1, 19, 48, 30, 6] [0, 1, 36, 149, 158, 55, 7]
The row sums
1, 2, 4, 10, 30, 104, 406, 1754, 8280, 42294, 231950, ...
are Paul Barry's A186021, Bell(n)*(20^n). Geoffrey Critzer commented: These are "the number of collections of subsets of {1,2,...,n1} that are pairwise disjoint."
The coefficients of the x_{0}based univariate polynomials are now in A257563. What is their combinatorial meaning?
♦ The Bell transform based on partitions
We did not yet consider the question of efficiency of our implementation, to which we turn now.
If we delete the intermediate steps which we introduced for the sake of exposition above we can sum up our implementation as:
def bell_polynomial(n): X = var(['x_'+str(i) for i in (0..n+1)]) @cached_function def T(n, k): if k==0: return k^n return sum(binomial(n1,j1)*T(nj,k1)*X[j] for j in (0..nk+1)).expand() return sum(T(n,k) for k in (0..n))
The nice thing is that we can speed up things by explicitly summing over partitions of n instead of using the recursive form. This is somewhat surprising since the relevant formula looks quite formidable; the reader will find it as the first formula on Wikipedia's entry on Bell polynomials.
(Note that Wikipedia calls 'complete Bell polynomials' what we call 'Bell polynomials' and Mathworld calls 'Bell polynomials' what we call 'univariate Bell polynomials'.)
We start with the simpler case of the univariate Bell polynomials and look at an implementation which is essentially the implementation of Sage (GPL licensed), written by Blair Sutton.
def partition_based_univariate_bell_polynomial(n): polynom = x^n fn = factorial(n) for k in (0..n1): result = 0 for p in Partitions(n, length=k): factorial_product = 1 power_factorial_product = 1 for part, count in p.to_exp_dict().iteritems(): factorial_product *= factorial(count) power_factorial_product *= factorial(part)**count coefficient = fn//(factorial_product*power_factorial_product) result += coefficient polynom += result*x^k return polynom
So let's benchmark!
%timeit univariate_bell_polynomial(10) 5 loops, best of 3: 105 ms per loop %timeit partition_based_univariate_bell_polynomial(10) 5 loops, best of 3: 14.1 ms per loop %timeit univariate_bell_polynomial(20) 5 loops, best of 3: 632 ms per loop %timeit partition_based_univariate_bell_polynomial(20) 5 loops, best of 3: 79.4 ms per loop %timeit univariate_bell_polynomial(30) 5 loops, best of 3: 2.15 s per loop %timeit partition_based_univariate_bell_polynomial(30) 5 loops, best of 3: 1.16 s per loop
Thus for small n we observe a speedup between 2 and 8.
Back to math. Our main goal is to explore the Bell transform. Therefore we skip the details concerning the implementation of the multivariate Bell polynomials (as polynomials) and give directly the formalism in its most useful form: as a sequence to triangle transformation. To make our setup more generic we now change from lists as input to functions.
def bell_transform(f, n): # partition_based row = [] fn = factorial(n) for k in (0..n): result = 0 for p in Partitions(n, length=k): factorial_product = 1 power_factorial_product = 1 for part, count in p.to_exp_dict().iteritems(): factorial_product *= factorial(count) power_factorial_product *= factorial(part)**count coefficient = fn//(factorial_product*power_factorial_product) result += coefficient*prod([f(i1) for i in p]) row.append(result) return row
♦ The Bell matrix
For example let's look at the family of Stirlingcycle/Lah numbers, which are the Bell transforms of the functions "n → factorial(n+j)" (for some fixed j).
for n in (0..7): bell_transform(factorial, n) [1] [0, 1] [0, 1, 1] [0, 2, 3, 1] [0, 6, 11, 6, 1] [0, 24, 50, 35, 10, 1] [0, 120, 274, 225, 85, 15, 1] [0, 720, 1764, 1624, 735, 175, 21, 1] Unsigned Stirling numbers of the first kind, A132393. for n in (0..7): bell_transform(lambda x: factorial(x+1), n) [1] [0, 1] [0, 2, 1] [0, 6, 6, 1] [0, 24, 36, 12, 1] [0, 120, 240, 120, 20, 1] [0, 720, 1800, 1200, 300, 30, 1] [0, 5040, 15120, 12600, 4200, 630, 42, 1] Unsigned Lah numbers, cf. A111596. for n in (0..7): bell_transform(lambda x: factorial(x+2), n) [1] [0, 2] [0, 6, 4] [0, 24, 36, 8] [0, 120, 300, 144, 16] [0, 720, 2640, 2040, 480, 32] [0, 5040, 25200, 27720, 10320, 1440, 64] [0, 40320, 262080, 383040, 199920, 43680, 4032, 128] Unsigned coefficients for rewriting generalized falling factorials into ordinary falling factorials, A136656.
The next step is a trival one: We collect these lists in a matrix, the Bell matrix, and call the function given as input the generator of the Bell matrix. With SageMath we can write this:
def bell_matrix(generator, dim): row = lambda n: bell_transform(generator, n) return matrix(ZZ, [row(n)+[0]*(dimn1) for n in srange(dim)])
Alternatively, if we want to avoid the matrix constructor and base the Bell transform on its recursive form we can compute the Bell matrix with this function:
def Bell_Matrix(generator, dim): A = [[0]*(n+1) for n in range(dim)] for n in range(dim): A[n][0] = 1 if n==0 else 0 if n>0: A[n][1] = generator(n1) for k in (2..n): A[n][k] = sum(binomial(n1,j1)*A[nj][k1]*A[j][1] for j in (1..nk+1)) return A
♦ The inverse Bell transform
The inverse Bell transform is the map generator → inverse(bell_transform(generator)). Here we consider 'bell_transform(generator)' as an infinite lower matrix of which we calculate the matrix inverse.
For example using the factorial function as the generator we can write:
bell_matrix(factorial, 8) [ 1 0 0 0 0 0 0 0] [ 0 1 0 0 0 0 0 0] [ 0 1 1 0 0 0 0 0] [ 0 2 3 1 0 0 0 0] [ 0 6 11 6 1 0 0 0] [ 0 24 50 35 10 1 0 0] [ 0 120 274 225 85 15 1 0] [ 0 720 1764 1624 735 175 21 1]
Then the inverse Bell transform can be computed as:
bell_matrix(factorial, 8).inverse() [ 1 0 0 0 0 0 0 0] [ 0 1 0 0 0 0 0 0] [ 0 1 1 0 0 0 0 0] [ 0 1 3 1 0 0 0 0] [ 0 1 7 6 1 0 0 0] [ 0 1 15 25 10 1 0 0] [ 0 1 31 90 65 15 1 0] [ 0 1 63 301 350 140 21 1]
Another example are the coefficients of the Abel polynomials A137452 which have the enumeration of the natural numbers as their generator:
bell_matrix(lambda n: n+1, 8).inverse() [ 1 0 0 0 0 0 0 0] [ 0 1 0 0 0 0 0 0] [ 0 2 1 0 0 0 0 0] [ 0 9 6 1 0 0 0 0] [ 0 64 48 12 1 0 0 0] [ 0 625 500 150 20 1 0 0] [ 0 7776 6480 2160 360 30 1 0] [ 0 117649 100842 36015 6860 735 42 1]
This is the inverse of the triangle of idempotent numbers A059297 :
bell_matrix(lambda n: n+1, 8) [ 1 0 0 0 0 0 0 0] [ 0 1 0 0 0 0 0 0] [ 0 2 1 0 0 0 0 0] [ 0 3 6 1 0 0 0 0] [ 0 4 24 12 1 0 0 0] [ 0 5 80 90 20 1 0 0] [ 0 6 240 540 240 30 1 0] [ 0 7 672 2835 2240 525 42 1]
Alternatively, if we want to avoid the matrix constructor and base the inverse Bell transform on its recursive form we can use the following function:
def Inverse_Bell_Matrix(f, dim): A = Bell_Matrix(f, dim) M = [[0]*(n+1) for n in range(dim)] for n in range(dim): M[n][n] = 1/A[n][n] for k in range(n1,0,1): M[n][k] = sum(A[i][k]*M[n][i] for i in range(n,k,1))/A[k][k] return M Inverse_Bell_Matrix(factorial, 7)
♦ The Bell inverse of a sequence
Given an invertible Bell matrix B with generator g then we call the generator of the inverse of B the Bell inverse of g. Thus we are looking here at a sequence to sequence transformation.
With Sage we can implement this as:
# Given a sequence f returns the inverse Bell sequence of f. def Inverse_Bell_Sequence(f, dim): M = Inverse_Bell_Matrix(f, dim) return [M[n][1] for n in (1..dim1)]
Some examples:
# [1, 1, 1, 1, 1, 1, 1, 1, ...] maps to # [1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, ...] print Inverse_Bell_Sequence(lambda n: (1)^n, 9) # [1, 2, 2, 0, 0, 0, 0, 0, 0, ...] maps to # [1, 2, 10, 80, 880, 12320, 209440, 4188800, ...] print Inverse_Bell_Sequence(lambda n: [1,2,2][n] if n<3 else 0, 9) # [1, 2, 1, 0, 0, 0, 0, 0, 0, ...] maps to # [1, 2, 11, 100, 1270, 20720, 413000, 9726640, ...] print Inverse_Bell_Sequence(lambda n: [1,2,1][n] if n<3 else 0, 9)
More general: Let g be a 0based sequence then we require g(0) = 0 and g(1) != 0 and write a = g(1), b = g(2), c = g(3), d = g(4) to simplify the notation. The Bell inverse of g then starts
1/a, b/a^3, (c*a+3*b^2)/a^5, (d*a^210*b*c*a+15*b^3)/a^7, ...
If we assume a = 1, which will be often the case, this reduces to:
1, b, c+3*b^2, d+10*b*c15*b^3, ...
1 
b 
3*b^2  c 
15*b^3 + 10*b*c  d 
105*b^4  105*b^2*c + 10*c^2 + 15*b*d  e 
Therefore the skeleton of the coefficients is
1 
1 
3, 1 
15, 10, 1 
105, 105, 10, 15, 1 
Looking up this table in the OEIS we find Wolfdieter Lang's A176740 called 'Inversion of e.g.f. formal power series' and his paper E.g.f. Lagrange inversion.
There are two differences: Lang's table misses the leading '1' and the order of coefficients is different. For our purposes the latter is of no importance.
♦ The coefficients of the Bell polynomials
But why stop half way? So let's look at all the coefficients of the inverse Bell polynomials. The coefficients in the table above appear here in column 1.
n\k  0  1  2  3  4  5 
0  [1]  
1  [0]  [1]  
2  [0]  [1]  [1]  
3  [0]  [3, 1]  [3]  [1]  
4  [0]  [15, 10, 1]  [15, 4]  [6]  [1]  
5  [0]  [105, 105, 10, 15, 1]  [105, 60, 5]  [45, 10]  [10]  [1] 
It seems somewhat surprising that this table was not yet in the database. Now it is A268442. Let's also add the table of the coefficients of the Bell polynomials A268441 (also A036040 with different order and column 0 missing).
n\k  0  1  2  3  4  5 
0  [1]  
1  [0]  [1]  
2  [0]  [1]  [1]  
3  [0]  [1]  [3]  [1]  
4  [0]  [1]  [3, 4]  [6]  [1]  
5  [0]  [1]  [10, 5]  [15, 10]  [10]  [1] 
There is a simple check of the tables: replacing the lists by their sums reduces the triangles to the Stirling number tables: to the Stirling numbers of first kind and to the Stirling numbers of the second kind respectively.
The algorithm for building these tables is given below, implemented with Sage.
def bell_polynomial_matrix(dim, inverse): def bell_polynomial(n): X = var(['x'+str(i) for i in (0..dim)]) @cached_function def T(n, k): if k==0: return k^n return sum(binomial(n1,j1)*T(nj,k1)*X[j1] for j in (0..nk+1)).expand() return [T(n,k) for k in (0..n)] A = [[f for f in bell_polynomial(n)] for n in range(dim)] if not inverse: return A M = [[0 for k in (0..n)] for n in range(dim)] for n in range(dim): M[n][n] = 1/A[n][n] for k in range(n1,1,1): M[n][k] = expand(sum(A[i][k]*M[n][i] for i in range(n,k,1))/A[k][k]) return M def coefficient_matrix(M): def coefficient(p): c = SR(p).fraction(ZZ).numerator().coefficients() return [0] if not c else c return [[coefficient(p) for p in M[n]] for n in range(len(M))]
The two tables above are then generated by the function call
M = bell_polynomial_matrix(8, inverse=true) coefficient_matrix(M)
respectively by the call
L = bell_polynomial_matrix(8, inverse=false) coefficient_matrix(L)
♦ Bell numbers of higher order
Consider the sequence
S_{0} → T_{0} → S_{1} → T_{1} → S_{2} → T_{2} → ...
Here S_{n} → T_{n} indicates the Bell transform mapping a sequence S_{n} to the triangle T_{n} and Tn → S_{n+1} the operator associating a triangle with the sequence of its row sums. For example if we start this iteration with the sequence <1,1,1,...> we get
S0 = A000012 = <1,1,1,...> T0 = A048993 # Stirling subset numbers, S1 = A000110 # Bell numbers, T1 = A264428 # Bell transform of Bell numbers, S2 = A187761 # secondorder Bell numbers, T2 = A264430 # Bell transform of secondorder Bell numbers, S3 = A264432 # thirdorder Bell numbers.
Similarly if we start the iteration with the sequence <1,1,1,1,...> we get the higher order complementary Bell numbers.
Sequence  Bell transform 
Inverse Bell transform 

[1,1,1,1,...]  A000012  A048993  A048994 
Bell numbers  A000110  A264428  A264429 
Second order Bell numbers  A187761  A264430  A264431 
Third order Bell numbers  A264432  A264433  A264434 
[1,1,1,1,...]  A033999  A080417  A132393 
Complementary Bell numbers  A000587  A264435  A264436 
Second order complementary Bell numbers 
A265023 
A nice observation is that Joerg Arndt's award winning A187761 (number of maps f: [n] → [n] with f(x) ≤ x and f(f(x)) = f(f(f(x)))) here turns up as the secondorder Bell numbers.
def bell_second_order(generator, n): G = [generator(k) for k in range(n)] row = lambda n: bell_transform(n, G) S = [sum(row(k)) for k in range(n)] return bell_transform(n, S)
Arndt's sequence can now be easily computed from the constant sequence 1,1,1,... as:
[sum(s) for s in [bell_second_order(lambda k: 1, n) for n in range(10)]]
The sequence is also column 1 in the triangle A264430.
[bell_second_order(bell_number, n) for n in range(8)] [[1], [0, 1], [0, 1, 1], [0, 2, 3, 1], [0, 6, 11, 6, 1], [0, 23, 50, 35, 10, 1], [0, 106, 268, 225, 85, 15, 1], [0, 568, 1645, 1603, 735, 175, 21, 1]]
♦ Associated Stirling subset numbers
As we have seen in the table above the Stirling subset numbers are generated by applying the Bell transform to the simplest sequence of positive numbers: the all 1's sequence.
On the other hand one of the simplest methods to transform a sequence is the Hilbert's Hotel Transform (HHT) also known as the 'shift right operation' in CS communities: move simultaneously cell n to cell n+1 and 0 to cell 0.
So what happens if we apply this idea to 1,1,1,... ? We define in Sage:
hht = lambda k: lambda n: 1 if n>=k else 0 for n in (0..5): print bell_matrix(hht(n), 12)
In terms of Anumbers we get for n = 0, 1, 2,... the triangles:
A048993, A008299, A059022, A059023, A059024, A059025.
which are called the (n+1)associated Stirling subset numbers following Comtet.
However compared to our construction in the last section the construction of this sequence of sequences has a disadvantage: apart from the case n = 0 the resulting matrices are singular; in other words, the inverse Bell transform cannot be applied to these generators.
♦ The Bell transform of multifactorials
An overview about multifactorials can be found in my July 2011 blog post. In Sage parlance the general definition is:
multifactorial = lambda a,b: lambda n: prod(a*k + b for k in (0..n1))
Here the multifactorials depends on two parameter, a and b. For example multifactorial(2, 2) are the double factorials of even numbers (2*n)!! = 2^n*n!, A000165. Note that multifactorial(2, 2) is a function, evaluated at n=4 gives multifactorial(2, 2)(4) = 384.
Now we can study the Bell transform of multifactorials.
for a in (1..4): for b in (1..a): print "Bell matrix generated by multifactorial", (a,b) bell_matrix(multifactorial(a,b), 6) print "Inverse Bell matrix generated by multifactorial", (a,b) inverse_bell_matrix(multifactorial(a,b), 6)
The table below shows the output in terms of Anumbers. From the 20 basic cases 17 have been in the OEIS (the missing three were added by the author). Looking closer at the entries it becomes obvious that the Bell transform of multifactorials was never studied systematically before.
Type  Multi factorial 
Bell transform 
Inverse Bell transform 
MF_{1,1}  A000142  A132393  A048993 
MF_{2,1}  A001147  A132062  A122848 
MF_{2,2}  A000165  A039683  A075497 
MF_{3,1}  A032031  A203412  A265605 
MF_{3,2}  A007559  A004747  A049404 
MF_{3,3}  A008544  A051141  A075498 
MF_{4,1}  A007696  A265606  A265604 
MF_{4,2}  A001813  A119274  A119275 
MF_{4,3}  A008545  A265606  A049410 
MF_{4,4}  A047053  A051142  A075499 
Some special cases of the Bellinverse of multifactorials are:
bell_inverse(MF_{n, n}) = n^k for k≥0.
bell_inverse(MF_{n, n1}) = falling_factorial(n  1, k) for k≥0.
bell_inverse(MF_{n, 1}) = MF_{n1, n2} (for n≥3).
bell_inverse(MF_{2*n, n}) = [1, n, 0, 0, 0,...].
bell_inverse(MF_{2*n, 2*n2}) = A161381_row(n1).
♦ The Bell transform of powers
The power functions can be seen as the degenerated cases of multifactorials: as those cases where a = 0 in the definition of multifactorials.
multifactorial = lambda a,b: lambda n: prod(a*k + b for k in (0..n1))
These cases have been studied by Wolfdieter Lang. Lang calls the Bell transform of powers "Stirling2 triangles with scaled diagonals".
for n in range(6): bell_matrix(lambda k: n^k, 7)
Lang calls the signed inverse Bell transform of powers "Generalized Stirling number triangles of first kind".
for n in (0..5): bell_matrix(lambda k: n^k, 7).inverse()
n  Powers of n 
Bell transform 
Inverse Bell transform 
0  A000007  A023531  A023531 
1  A000012  A048993  A132393 
2  A000079  A075497  A039683 
3  A000244  A075498  A051141 
4  A000302  A075499  A051142 
5  A000351  A075500  A051150 
6  A000400  A075501  A051151 
7  A000420  A075502  A051186 
8  A001018  A075503  A051187 
9  A001019  A075504  A051231 
The Bellinverse of the power functions are
bell_inverse(n^k) = n^k*k! for k≥0.
♦ The Bell transform of rising factorials
risingfactorial = lambda n: lambda k: rising_factorial(n, k) for n in (0..7): print bell_matrix(risingfactorial(n), 7) print inverse_bell_matrix(risingfactorial(n), 7)
The Bell transform of the rising factorial and its inverse lead to:
n  Rising factorial 
Bell transform 
Inverse Bell transform 
0  A000007  A023531  A023531 
1  A000142  A132393  A048993 
2  A000142*  A111596  A111596 
3  A001710  A046089  A035342 
4  A001715  A049352  A035469 
5  A001720  A049353  A049029 
6  A001725  A049374  A049385 
7  A001730  A134141  A092082 
The Bellinverse of rising factorials are multifactorials:
bell_inverse(risingfactorial(n)) = MF(n1,1)(k+1) for k≥0.
♦ The Bell transform of falling factorials
Similarly for the Bell transform of falling factorials:
fallingfactorial = lambda n: lambda k: falling_factorial(n, k) for n in (0..5): print bell_matrix(fallingfactorial(n), 7) print inverse_bell_matrix(fallingfactorial(n), 7)
n  Falling factorial A008279_row 
Bell transform 
Inverse Bell transform 
0  row(0)  A023531  A023531 
1  row(1)  A049403 A122848 
A132062 
2  row(2)  A049404  A004747 
3  row(3)  A049410  A000369 
4  row(4)  A049424  A011801 
5  row(5)  A049411  A013988 
The Bellinverse of falling factorials are multifactorials:
bell_inverse(fallingfactorial(n)) = MF(n+1,n)(k) for k≥0.
♦ The Bell transform of monotonic factorials
There seems to be no standard name for the function
monotonicfactorial = lambda r: lambda n: rising_factorial(r, n)/factorial(n)
I therefore dubbed it monotonic factorial for the purpose of this post (as it counts the number of monotone words with length n over an ralphabet; most often it is written as binomial(n+r1, n) ).
for n in (0..4): print bell_matrix(monotonicfactorial(n), 7) print inverse_bell_matrix(monotonicfactorial(n), 7)
Astonishing little is found in the OEIS about this class of sequences, with the notable exception of Roger L. Bagula's A137452 (Olivier Gérard's A061356) which is the inverse Bell transform of A000169. Certainly more triangles of this type should be entered into the OEIS.
♦ An implementation for Maple
The Maple implementation of the Bell matrix given below is compatible with the interface and conventions of Sloane's Maple transforms.
# BELL: Computes the Bell matrix of a sequence. # Given the list [a(0),...,a(dim)] returns the Bell matrix # with dimension dim+1 of the sequence a. BELL := proc(a) local M, dim, n, k: if whattype(a) <> list then RETURN([]) fi: dim := nops(a); M := Matrix(dim, shape=triangular[lower]); M[1,1] := 1; for n from 1 to dim1 do M[n+1,2] := a[n] od; for n from 1 to dim1 do for k from 1 to n do M[n+1,k+1] := add(binomial(n1,j1)*M[nj+1,k]*M[j+1,2], j=1..nk+1) od od; RETURN(M) end: BELLi:= proc(a) if whattype(a) <> list then RETURN([]) fi: RETURN(linalg[inverse](BELL(a))) end: # Examples: a := [seq(1,n=0..8)]: BELL(a); # A048993 b := [seq(n!,n=0..8)]: BELL(b); # A132393 c := [seq(n,n=1..9)]: BELL(c); # A059297 d := [seq((1)^n*(n+1)^n,n=0..8)]: BELL(d); # A137452 e := [seq(`if`(n::odd,0,2*n!),n=0..8)]: BELL(e); # A137513 f := [seq(`if`(n<2,(1)^n,0),n=0..8)]: BELL(f); # A104556 g := [seq(doublefactorial(2*n1),n=0..8)]: BELL(g); # A001497
♦ Index of triangles generated by the Bell transform
The format of the list is [An[Ak]] which means that the triangle An is generated by the sequence Ak. If Ak is not in the OEIS than Ak might point to a related sequence or to dummy sequence Axxxxxx.
[A000369[A008545]] [A001497[A001147]] [A004747[A008544]] [A008275[A133942]] [A008277[A000012]] [A008296[A265313]] [A008297[A133942]] [A008298[A038048]] [A011801[A008546]] [A013988[A008543]] [A023531[A000007]] [A035342[A001147]] [A035469[A007559]] [A038455[A006963]] [A039621[A177885]] [A039683[A000165]] [A039692[A039647]] [A039810[A000110]] [A039811[A000258]] [A039812[A000307]] [A039813[A000357]] [A039814[A003713]] [A039815[A000268]] [A039816[A000310]] [A039817[A000359]] [A046089[A001710]] [A048176[A051262]] [A048786[A034177]] [A048993[A000012]] [A048994[A133942]] [A049029[A007696]] [A049218[A005359]] [A049352[A001715]] [A049353[A001720]] [A049374[A001725]] [A049385[A008548]] [A049403[A000161]] [A049404[A000925]] [A049410[A004552]] [A049411[A008279]] [A049424[A265609]] [A051141[A032031]] [A051142[A047053]] [A051150[A052562]] [A051151[A047058]] [A051186[A051188]] [A051187[A051189]] [A051231[A051232]] [A059297[A000027]] [A059298[A000027]] [A059419[A009006]] [A060281[A001865]] [A061356[A000169]] [A075497[A000079]] [A075498[A000244]] [A075499[A000302]] [A075500[A000351]] [A075501[A000400]] [A075502[A000420]] [A075503[A001018]] [A075504[A001019]] [A075505[A011557]] [A075525[A265024]] [A078521[A038048]] [A079621[A002866]] [A079638[A029767]] [A079639[A006252]] [A079640[A007840]] [A079641[A000629]] [A079642[A089064]] [A086915[A052849]] [A088729[A000670]] [A088814[A000262]] [A092082[A008542]] [A104556[A003475]] [A105278[A000142]] [A105599[A000272]] [A105786[A000272]] [A105819[A055860]] [A106239[A057500]] [A111246[A130716]] [A111593[A155585]] [A111594[A005359]] [A111596[A155456]] [A119274[A001813]] [A119275[A130706]] [A121408[A177145]] [A122848[Axxxxxx]] [A122850[A001147]] [A125553[A208529]] [A129062[A000629]] [A130123[A000038]] [A130191[A000110]] [A130534[A000142]] [A131222[A029767]] [A132056[A045754]] [A132062[A001147]] [A132393[A000142]] [A134141[A001730]] [A135338[A155456]] [A135494[A153881]] [A136590[A136591]] [A136595[A048287]] [A136630[A000035]] [A136656[A155456]] [A137312[A052849]] [A137320[A066459]] [A137339[A052560]] [A137378[A052510]] [A137431[Axxxxxx]] [A137433[Axxxxxx]] [A137452[A213236]] [A137513[Axxxxxx]] [A143395[A000225]] [A143543[A001187]] [A144402[A000161]] [A144633[A144636]] [A144634[A144636]] [A144644[A001899]] [A145520[A000040]] [A147308[A027641]] [A147309[A022902]] [A147311[A122045]] [A147312[A122045]] [A147315[A000111]] [A151509[Axxxxxx]] [A151511[Axxxxxx]] [A166317[A002436]] [A166318[A002436]] [A171996[Axxxxxx]] [A171998[Axxxxxx]] [A174893[A000142]] [A184962[A000670]] [A185285[A004123]] [A185296[Axxxxxx]] [A185415[Axxxxxx]] [A185419[A143523]] [A185422[A080635]] [A185690[A056594]] [A185951[A193356]] [A186366[A000111]] [A187082[Axxxxxx]] [A187084[Axxxxxx]] [A188062[Axxxxxx]] [A188066[Axxxxxx]] [A188832[A009843]] [A189898[A003027]] [A191249[A062738]] [A194938[A039647]] [A195204[A076726]] [A195205[Axxxxxx]] [A202183[Axxxxxx]] [A202184[Axxxxxx]] [A202185[Axxxxxx]] [A202189[Axxxxxx]] [A202190[Axxxxxx]] [A203412[A007559]] [A209849[A155585]] [A215771[A001710]] [A215861[A215851]] [A217756[A129271]] [A223511[A045755]] [A223512[A045756]] [A223513[A045757]] [A223514[Axxxxxx]] [A223515[Axxxxxx]] [A223516[Axxxxxx]] [A223517[Axxxxxx]] [A223518[Axxxxxx]] [A223522[Axxxxxx]] [A225171[A225170]] [A227450[A007395]] [A228534[A009444]] [A228550[A033678]] [A228859[A001832]] [A228892[A002027]] [A247232[A001929]] [A256041[A005212]] [A256042[A256033]] [A256892[A000262]] [A256893[A000670]] [A259286[Axxxxxx]] [A264428[A000110]] [A264429[Axxxxxx]] [A264430[A125273]] [A264431[Axxxxxx]] [A264433[A047889]] [A264434[A011634]] [A264435[A000587]] [A264436[A007549]] [A265314[Axxxxxx]] [A265602[Axxxxxx]] [A265608[A000262]]
In particular this table shows that the coefficients of the Touchard or Bell polynomials (A048993), the Abel polynomials (A137452), the MittagLeffler polynomials (A137513), the modified Hermite polynomials (A104556) and the Bessel polynomials (A001497) all are Bell matrices generated by elementary sequences.
♦ References
 D. E. Knuth, Convolution polynomials, Mathematica J. 2.1 (1992), no. 4, 6778.
Accompanying, there is an Jupyter Notebook which was developed on SageMath Cloud. SageMath is a free open source alternative to Magma, Maple, Mathematica and Matlab. Get a free account on SageMath Cloud, download this notebook, upload it to SMC and start exploring!