Inverse of multivariate function power series

Question

Consider a vector function $\mathbf y = \mathbf{f}(\mathbf x)$, where $\mathbf y,\mathbf x \in \mathbb R^N$. To be more explicit,

$$y_1 = f_1(x_1,\dots,x_N)$$ $$y_2 = f_2(x_1,\dots,x_N)$$ $$...$$ $$y_N = f_N(x_1,\dots,x_N)$$

I will assume for simplicity that $\mathbf f(\mathbf 0) = \mathbf 0$, and that $\mathbf f$ is locally invertible.

Suppose I know the form of the function $\mathbf f$, and in particular the power series,

$$y_n = \sum_i A_i^n x_i + \sum_{ij} A_{ij}^n x_ix_j + \sum_{ijk} A_{ijk}^n x_ix_jx_k + \dots$$

I want to find the power series of the inverse function, that is,

$$x_n = \sum_i B_i^n y_i + \sum_{ij} B_{ij}^n y_iy_j + \sum_{ijk} B_{ijk}^n y_iy_jy_k + \dots$$

How can I compute the $B$'s from the $A$'s?

I would suspect there is some literature on the relations between the coefficients $A$ and $B$ here, but I have not been able to find it. Except for the one-dimensional case ($N=1$), where we have explicit expressions involving Bell numbers. Are there similar results / formalisms for the multivariate case?

For example, it is easy to see that $B_i^n$ is a matrix inverse of $A_i^n$ (because the Jacobian of an inverse function is the inverse of the Jacobian). What about the higher-order terms?

Note: In practice I need to compute the B's only for small-order.

Answer 1

Statement of problem. Consider the $n\times n$ nonlinear system $Y=\vec y = f(\vec x)$ that is assumed to satisfy $0= f(0)$. Assume the each component of $f$ is a multivariable power series in the components of $\vec x$. Assume that the Jacobian matrix $J= [df]_0=[dY/dX]|_{X=0}$ at the origin is invertible. The goal is to find the power series expansion of the inverse map $\vec x = P(\vec y)$. Note that $Y= f(X)=J.X + Q(X)$ where the series for the correction term $Q(X)$ is quadratic or higher. (Henceforth the period . means matrix product).

Reductions. After a constant-coefficient linear change of variables $U= J^{-1} X$ we can restate the problem as $Y= U + \tilde Q(U)$ for some new correction term $\tilde Q(U)=Q(JU)$.

Thus WLOG we assume that the Jacobian matrix is the identity matrix.

Caratheodory format.

Returning to our original notation, we now assume

$Y= X+ Q(X)$ and note that quadratic correction can be expressed in "Caratheodory format" as $Q(X)= [A(X)].X$ where the square matrix $A= [A(X)]$ contains terms that are known polynomials in $X$ that are of degree 1 or higher. (I use this phrasing because Caratheodory pointed out that this factorization can be used to give an elegant definition of differentiability. Note the $A$ is not unique.)

We seek the inverse map, also assumed to written in Caratheodory format, as $X= Y+ B(Y).Y$ where the square matrix of power series terms $B(Y)$ must be determined. (Each matrix entry of $B(Y)$ is an unknown polynomial of degree one or higher).

Derivation of a fixed-point equation for $B$

Substituting the series $X= Y+ B(Y).Y$ inside the series $Y= X+ A(X).X$ we require that we obtain the identity map:

$Y= I.Y= (I.Y+ B(Y).Y) + [A (Y+B(Y).Y)] .(I.Y+B(Y).Y)$ This can be satisfied if the matrix equation

$$0= B(Y) + [A(Y+[B(Y)]Y] (I+ [B(Y)]$$ holds true identically in $Y$. (This identity is what you get by erasing the column vector $Y$ that occurs everywhere on the right side of the former equation).

Equivalently $B(Y)= -A[Y+[B(Y)]Y] (I+ [B(Y)])$.

This is a fixed point equation for the matrix $B(Y)$.

A recursion that boosts the degree of accuracy

It is instructive to expand $B(Y) =\sum_k b_k(Y)$ as a sum of terms of progressively higher total degree.

(Trick: When implementing this in Mathematica, it is helpful to organize the calculation by replacing $Y-> tY$ throughout, and match corresponding powers of the scalar variable $t$.) The identity

$$(1)\qquad B(tY) = - A(tY+[B(tY).(tY)). (I + B(tY))$$ has the feature that the if we substitute the truncated partial sum expansion $$(2)\qquad B_n(t)= \sum_{k=1}^n b_k(Y) t^k$$

into (1), then the term that appears on the left provides one higher-degree term than the right.

So we define recursively $$ (3) \qquad B_{n+1}(t)= - A(tY+[B_n(tY).(tY)). (I + B_n(tY))$$ and this determines more and more terms in the power series expansion. (The initial value in the recursion is the zero matrix.)

Below we illustrate the algorithm in an example.

Consider $Y= f(X) = X+ A(X).X$ whose component formulas are

$[y1,y2,y3]= [x1,x2,x3]+ [3 \text{x1} \text{x3}-\text{x1} \text{x2},\text{x3} \left(\text{x1}+\text{x2}^2\right)+3 \text{x1} \text{x3}+\text{x2}^2,2 \text{x1} \text{x2}^2+\text{x3} (\text{x1}+\text{x2} \text{x3})+\text{x1} \text{x3}^2]$

This factors as $X+ A(X).X$ in many ways, one of which is $$A(X)=\begin{array}{ccc} -2 \text{x2} & \text{x1} & 3 \text{x1} \\ 3 \text{x3} & \text{x2} & \text{x1}+\text{x2}^2 \\ \text{x3}^2 & 2 \text{x1} \text{x2} & \text{x1}+\text{x2} \text{x3} \\ \end{array}$$

The Mathematica code that defines the perturbation matrix A (called p in the code) and the recursion function (3) is given below.

After the first three steps of recursion are completed, the truncated third-order answer $B_3$ is obtained. In the image below, we reveal the total answer in pieces, as linear $(b_1)$, quadratic $(b_2)$, and third order $b_3$:

Later I checked (by Mathematica) that the composition of the original series and the inverse series actually returned the identity map (modulo higher-order terms of fourth order).