Generating function for $(A+B)^n$ in terms of generating functions for $A^n$, $B^n$

Question

I need to express generating function of sequence $(A+B)^n$ in terms of generating functions of sequences $A^n$ or $B^n$.

For instance, consider the following matrix sequences.

$$ \begin{array} \ A^0,A^1,A^2,\ldots\\ B^0,B^1,B^2\ldots\\ \end{array} $$

We can consider a new sequence $(A+B)^n$ $$\ (A+B)^0,(A+B)^1,(A+B)^2,\ldots$$

Its generating function was derived here, in terms of generating function of $A^n$. Is there's functionality in Mathematica I could've reused to help obtain this expression?

$$ G_{A+B}=(I-x G_A B)^{-1} G_A $$

ClearAll["Global`*"];
SeedRandom[1];
d = 2;

ii = IdentityMatrix[d];
A = RandomInteger[{-10, 10}, {d, d}];
B = RandomInteger[{-10, 10}, {d, d}];
gA = GeneratingFunction[MatrixPower[A, n], n, 
   x]; (* same as Inverse[ii-x A] *)
gAB = Function @@ {{x}, 
    GeneratingFunction[MatrixPower[A + B, n], n, x]};
formula = Function @@ {{x}, Inverse[ii - x gA . B] . gA};

x = 2;
Print["Formula works: ", formula[x] == gAB[x] // Reduce] (* True *)

Answer 1

Your formula reminds the Sherman-Morrison formula: $$(A+b c^T)^{-1}=A^{-1}-\frac{A^{-1}b c^T A^{-1}}{1+c^TA^{-1}b},$$ but there are some important details. I am quite sure that one cannot derive from here a generating function for $(A+B)^n$, where $A$ and $B$ are general matrices. In fact, by looking at your code we see that you consider a special matrix $B=h h^T$ (I use here a convention that $h$ is a vector). $A$ is also diagonal in your case.

Even though it is elemental on the paper, I would be pleasantly surprised if you can derive the Sherman-Morrison equation (or more general Woodbury equation) with Mathematica. The problem is that MA knows

Inverse[Inverse[b]] // TensorExpand
(* b *)

But it is not capable of simplifying

Inverse[Inverse[a].Inverse[b]] // TensorReduce
(* MatrixPower[MatrixPower[a, -1].MatrixPower[b, -1], -1] *)

Should be just

b.a