Inverse of symmetric matrix plus identity matrix

Let $p(x)=\sum_{k=0}^na_kx^k$ (with $a_n=1$) and $q(x)$ be the characteristic polynomials of $A$ and $B=I+A$ respectively. By Cayley-Hamilton theorem, $B^{-1}=g(B)$ where \begin{align} g(x) &=\frac{q(0)-q(x)}{xq(0)}\\ &= \frac{p(-1)-p(x-1)}{xp(-1)}\\ &=\frac{-\sum_{j=0}^n a_j[(x-1)^j - (-1)^j]}{xp(-1)}\\ &=\frac{-\sum_{j=1}^n a_j \sum_{k=0}^{j-1}(-1)^{j-1-k}(x-1)^k} {\sum_{j=0}^n (-1)^ja_j}\\ &=\frac{-\sum_{j=0}^{n-1} a_{j+1} \sum_{k=0}^j(-1)^{j-k}(x-1)^k}{\sum_{j=0}^n (-1)^ja_j}. \end{align} Therefore, $(I+A)^{-1}$ can be expressed in terms of $A$ by $f(A)$, where \begin{align} f(x) &=\frac{-\sum_{j=0}^{n-1} a_{j+1} \sum_{k=0}^j(-1)^{j-k}x^k} {\sum_{j=0}^n (-1)^ja_j}\\ &=\frac{-\sum_{k=0}^{n-1} x^k \sum_{j=k}^{n-1} (-1)^{j-k} a_{j+1}} {\sum_{j=0}^n (-1)^ja_j}\\ &=\frac{-\sum_{k=0}^{n-1} x^k \sum_{j=0}^{n-1-k} (-1)^j a_{j+k+1}} {\sum_{j=0}^n (-1)^ja_j}. \end{align} That is, $$ (I+A)^{-1}=\frac{-\sum_{k=0}^{n-1}\left[\sum_{j=0}^{n-1-k} (-1)^j a_{j+k+1}\right] A^k} {\sum_{j=0}^n (-1)^ja_j}.\tag{1} $$ It is well-known that each coefficient $a_i$ can be expressed in terms of the traces of the powers of $A$. More specifically, by Girard-Waring formula, if we define $s_k=\operatorname{tr}(A^k)$, then $$ a_{n-m} = \sum_{\substack{k_1+2k_2+\cdots+mk_m=m\\ k_1,k_2,\ldots,k_m\ge0}}(-1)^{k_1+k_2+\cdots+k_m}\frac{1}{k_1!k_2!\cdots k_m!}\left(\frac{s_1}1\right)^{k_1}\left(\frac{s_2}2\right)^{k_2}\cdots\left(\frac{s_m}m\right)^{k_m}.\tag{2} $$ In particular, when $n=5$, we have $a_5=1,\ a_0=-\det(A)$ and \begin{align} a_1 &=-\frac{s_4}4 + \frac{s_2^2}8 + \frac{s_1s_3}3 - \frac{s_1^2s_2}4 + \frac{s_1^4}{24},\\ a_2&= -\frac{s_3}3 + \frac{s_1s_2}2 - \frac{s_1^3}6,\\ a_3&= \frac{s_1^2-s_2}2,\\ a_4&= -s_1, \end{align} but surely, if $A^{-1}$ is also known, we would express $a_1$ as $\det(A)\operatorname{tr}(A^{-1})$ instead. Plug these into $(1)$, you can express $(I+A)^{-1}$ as a weighted sum of the powers of $A$ when $n=5$.

Let $A\in M_n(\mathbb{C})$, $spectrum(A)=(\lambda_i)$ and $B=I+A$ (we assume that $-1\notin spectrum(A)$); note that $spectrum(B)=(1+\lambda_i)$.

Proposition. $(I+A)^{-1}$ is a polynomial of degree $n-1$ in $A$, the coefficients of which, are rational functions of the $(trace(A^k))_{k\leq n}$.

Proof (for $n=5$). Let $(\sigma_k)$ be the elementary symmetric functions of the $(\lambda_i)$ and $(\tau_k)$ be the elementary symmetric functions of the $(1+\lambda_i)$. According to Hamilton-Cayley:

$B^{-1}=\dfrac{1}{\tau_5}(B^4-\tau_1B^3+\tau_2B^2-\tau_3B+\tau_4I_5)$.

Each $B^k$ is a polynomial of degree $k$ in $A$.

On the other hand, each $\tau_k$ is a symmetric function of the $(\lambda_i)$, then a polynomial in $\sigma_1,\cdots,\sigma_k$. Since each $\sigma_k$ is a polynomial in $trace(A),\cdots,trace(A^k)$ (use the Newton's formulae), $\tau_k$ is a polynomial in $trace(A),\cdots,trace(A^k)$.

Finally, $(I+A)^{-1}$ is a polynomial in $A$, the coefficients of which, are rational functions of the $(trace(A^k))_{k\leq 5}$.

Remark. Of course, we may replace the characteristic polynomial of $A$ with its minimal polynomial.

EDIT. Answer to @ancientmathematician ; cf. below a procedure (in Maple) that works for every dimension $m$ (I use the Bell function).

***Decomposition of $(I+A)^{-1}$; the result is in "res"; $a$ denotes $A\in M_m$ and $t[i]=t_i$ denotes $Trace(A^i)$.

restart;

with(LinearAlgebra);

m:=5:

TAU:=Vector(m): for n from 1 to m do U:=Matrix(n,n): for i to n-1 do U[i+1, i] := -1 end do; S := Vector(n, symbol = s); X := Vector(n, symbol = x); for i to n do x[i] := -factorial(i-1)*s[i] end do;

for i to n-1 do for k to i do U[i+1-k, i] := x[k]*binomial(n-i-1+k, k-1) end do end do;

for i to n do U[i, n] := x[n-i+1] end do;

TAU[n] := (-1)^n*Determinant(U)/factorial(n);

T:=Vector(n,symbol=t):

for k from 1 to n do

rr:= (1+y)^(k): ss:=m:

for i from 1 to k do ss:=ss+coeff(rr,y,i)*t[i]: od: s[k]:=ss: od:

TAU[n]:=expand(TAU[n]): x:='x':s:='s':t:='t': od: ;

***$(I+A)^{-1}$ as a polynomial in $A$

res := (1+a)^(m-1);

for i from m-2 by -1 to 0 do res := res+(-1)^(i-m+1)TAU[m-i-1](1+a)^i end do;

res := (-1)^(m-1)*collect(expand(res), a)/TAU[m];

***VERIFICATION

A := RandomMatrix(m,m);

B := 1/(IdentityMatrix(m)+A);

a := A; for i to m do t[i] := Trace(A^i) end do;

simplify(res-B);