Combinatorics of multivariate Faà di Bruno formula

Question

This question is a followup on this one on stackoverflow where i implement in python this issue. I am applying the Faà di Bruno formula to obtain the $\mathbf{i}$th derivatives of the function $f = ln \circ g$ in terms of the derivatives of $g$, a function with $n$ variables. This yield the following formula, denoting by $f^{(\mathbf i)}$ the $\mathbf i$th derivative of a function f

$$f^{(\mathbf i)}(t) = \sum\limits_{\pi \in \Pi(\mathbf{i})} (\lvert \pi \rvert -1) (-1)^{\lvert \pi \rvert -1} g(t)^{- \lvert \pi \rvert} \prod\limits_{B \in \pi} g^{(\mathbf{i}(B))}$$

where the notation follows this section of the wikipedia page, plus a small variation that recover the multiindex $\mathbf i(B)$ from the set of dimension indexes $B$.

The problem is that, since the exponential function has all it's derivatives equal, as soon as we want to derivate twice on the same variable, i obtain several $\pi$ inside $\Pi(\mathbf i)$ that are equivalent.

For exemple, set $\mathbf i = (2,1)$ in a two-dimensional problem. Then, the partitions inside $\Pi(\mathbf i)$ are:

• $(0,0,1)$
• $(1), (0,0)$
• $(0), (0,1)$
• $(0, 1), (0)$
• $(0) (0) (1)$

There is here a partition that appears twice (without taking into account the order of blocks inside a partition), and which will have the same value into the sum.

This happends a lot (run my code from stackoverflow to find out). Is there a way i could 'factorise' the formula ?

Edit: For $n = 1$, in a one-dimensional problem, the formula indeed factorises to the Bell polynomials, and the number of occurences of similar partitions are the Bell numbers. Is there somewhere a multivariate equivalent ?

Edit: This problem is equivalent to finding all multiset-partitions of a multiset. The solution in term of multivariate partial Bell polynomials exists (see article in a comment), but the polynomials in question are still hard to implement.

Answer 1

The Faa-di-Bruno formulae are really a clumsy ways of expressing the fact that (in the appropriate quotient algebra) the Taylor polynomial of a composition is the (formal) composition of the relevant Taylor polynomials. Working with formal Taylor polynomials should (I think) give you cleaner code and let you delegate most of the fiddly algebra/combinatorics to something like SymPy.

You should be able to get the relevant composition rule by doing a Google-search for "composition rule for Taylor polynimials" or similar, but you can find it Chapter 1 of Malgrange's Ideals of differentiable functions if you need.

Edit: I thought I'd provide the following rough sketch just to illustrate what I mean... and because I couldn't resist having a go!

import sympy
import math

def get_indeterminates(n_indeterminates):
    return [sympy.Poly('X{}'.format(i)) for i in range(1, n_indeterminates+1)]

def monomial_from_multi_index(multi_index, indeterminates):
    monomial = indeterminates[0]**0
    if type(multi_index) is int:
        multi_index = (multi_index,)
    for x, power in zip(indeterminates, multi_index):
        monomial *= x**power/math.factorial(power)
    return monomial

def taylor_polynomial_from_partials(partials, indeterminates):
    """
    'partials' should be a dictionary mapping integer tuples to floats
    """
    Tf = 0
    for k, fk in partials.items():
        Tf += fk * monomial_from_multi_index(k, indeterminates)
    return Tf

def get_composite_taylor_polynomial(derivatives_of_f, partials_of_g, n_variables, truncation_order):
    x = get_indeterminates(n_variables)
    y = get_indeterminates(1)
    Tf = taylor_polynomial_from_partials(derivatives_of_f, y)
    Tg = taylor_polynomial_from_partials(partials_of_g, x)
    Tfog = sympy.polys.polytools.compose(Tf, Tg)
    return