Inverse of symmetric matrix plus identity matrix Consider the symmetric, positive definite matrix $\mathbf{A}$. I'd like to find a general form for
$$(\mathbf{I} + \mathbf{A})^{-1}$$
that only involves $\mathbf{A}^{-1}$, i.e., no other inverse appears in the solution (as, for instance, in the Woodbury matrix identity).
I've tried to derive the inverse by hand but I could only obtain a result up to he $4 \times 4$ case as follows.


*

*$2 \times 2$:


$$(\mathbf{I} + \mathbf{A})^{-1} = \frac{\mathrm{det}(\mathbf{A}) \mathbf{A}^{-1} + \mathbf{I}}{\mathrm{det}(\mathbf{A}) + \mathrm{tr}(\mathbf{A}) + 1}$$


*

*$3 \times 3$:


$$(\mathbf{I} + \mathbf{A})^{-1} = \frac{\mathrm{det}(\mathbf{A}) \mathbf{A}^{-1} - \mathbf{A} + \big( \mathrm{tr}(\mathbf{A}) + 1 \big) \mathbf{I}}{\mathrm{det}(\mathbf{A}) + \mathrm{det}(\mathbf{A}) \mathrm{tr}(\mathbf{A}^{-1}) + \mathrm{tr}(\mathbf{A}) + 1}$$


*

*$4 \times 4$:


$$(\mathbf{I} + \mathbf{A})^{-1} = \frac{\mathrm{det}(\mathbf{A}) \mathbf{A}^{-1} + \mathbf{A}^2 - \mathrm{tr}(\mathbf{A}) \mathbf{A} - \mathbf{A} + \big( \tfrac{1}{2}(\mathrm{tr}(\mathbf{A})^2-\mathrm{tr}(\mathbf{A}^2)) + \mathrm{tr}(\mathbf{A}) + 1 \big) \mathbf{I}}{\mathrm{det}(\mathbf{A}) + \tfrac{1}{2}(\mathrm{tr}(\mathbf{A})^2-\mathrm{tr}(\mathbf{A}^2)) + \mathrm{det}(\mathbf{A}) \mathrm{tr}(\mathbf{A}^{-1}) + \mathrm{tr}(\mathbf{A}) + 1}$$
Can we find a general expression for higher dimensions that builds on the ones above? Even obtaining the $5 \times 5$ case is prohibitive, and so far I haven't been able to spot the pattern.
 A: 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$.
A: 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);
