[proofplan]
Set $B := A^{-1}E$, so the hypothesis becomes $\|B\|_{\mathcal L(X)} < 1$. We first prove the Neumann-series argument directly in the Banach algebra $\mathcal L(X)$: the series $\sum_{n=0}^{\infty}(-B)^n$ converges in [operator norm](/page/Operator%20Norm) and its sum is the two-sided inverse of $I_X+B$. We then factor $A+E$ as $A(I_X+B)$, preserving the order of factors, and multiply by the inverse of $I_X+B$ on the correct side to obtain the displayed formula.
[/proofplan]
[step:Introduce the perturbation inside the identity factor]
Let $I_X: X \to X$ denote the identity operator on $X$, and define the [bounded linear operator](/page/Bounded%20Linear%20Operator)
\begin{align*}
B: X &\to X
\end{align*}
by
\begin{align*}
B := A^{-1}E.
\end{align*}
Since $A^{-1} \in \mathcal L(X)$ and $E \in \mathcal L(X)$, their composition $B$ belongs to $\mathcal L(X)$. The hypothesis gives
\begin{align*}
\|B\|_{\mathcal L(X)} < 1.
\end{align*}
Using associativity of operator composition and the identity $AA^{-1}=I_X$, we have
\begin{align*}
A(I_X+B) = A + AA^{-1}E = A+E.
\end{align*}
Thus it remains to invert $I_X+B$ in $\mathcal L(X)$ and then multiply by $A^{-1}$ on the right.
[/step]
[step:Show that the Neumann series converges in operator norm]
For each $N \in \mathbb N$, define the partial-sum operator $S_N \in \mathcal L(X)$ by
\begin{align*}
S_N := \sum_{n=0}^{N}(-B)^n.
\end{align*}
We claim that $(S_N)_{N=1}^{\infty}$ is Cauchy in $\mathcal L(X)$ with respect to $\|\cdot\|_{\mathcal L(X)}$. If $M,N \in \mathbb N$ and $M>N$, then the triangle inequality and [submultiplicativity of the operator norm](/theorems/1054) give
\begin{align*}
\|S_M-S_N\|_{\mathcal L(X)} \le \sum_{n=N+1}^{M}\|(-B)^n\|_{\mathcal L(X)} \le \sum_{n=N+1}^{M}\|B\|_{\mathcal L(X)}^n.
\end{align*}
Since $0 \le \|B\|_{\mathcal L(X)} < 1$, the scalar geometric series $\sum_{n=0}^{\infty}\|B\|_{\mathcal L(X)}^n$ converges, so the right-hand side tends to $0$ as $N \to \infty$, uniformly in $M>N$. Hence $(S_N)$ is Cauchy.
Because $X$ is Banach, the bounded-operator space $\mathcal L(X)$ is complete under the operator norm. Therefore there exists $S \in \mathcal L(X)$ such that
\begin{align*}
S_N \to S
\end{align*}
in $\mathcal L(X)$ as $N \to \infty$. By definition,
\begin{align*}
S = \sum_{n=0}^{\infty}(-B)^n
\end{align*}
with convergence in operator norm.
[guided]
The purpose of this step is to construct a candidate inverse for $I_X+B$. The natural candidate is the geometric series
\begin{align*}
I_X - B + B^2 - B^3 + \cdots.
\end{align*}
To make this legitimate, we must prove that the series converges in the normed space where it is supposed to live, namely $\mathcal L(X)$ with the operator norm.
For each $N \in \mathbb N$, define
\begin{align*}
S_N := \sum_{n=0}^{N}(-B)^n \in \mathcal L(X).
\end{align*}
We estimate the distance between two partial sums. Let $M,N \in \mathbb N$ with $M>N$. Then
\begin{align*}
S_M-S_N = \sum_{n=N+1}^{M}(-B)^n.
\end{align*}
Taking the operator norm, applying the triangle inequality, and then applying submultiplicativity of the operator norm repeatedly, we obtain
\begin{align*}
\|S_M-S_N\|_{\mathcal L(X)} \le \sum_{n=N+1}^{M}\|(-B)^n\|_{\mathcal L(X)} \le \sum_{n=N+1}^{M}\|B\|_{\mathcal L(X)}^n.
\end{align*}
The assumption $\|B\|_{\mathcal L(X)}<1$ is used exactly here: it says that the scalar geometric series with ratio $\|B\|_{\mathcal L(X)}$ converges. Therefore the tail
\begin{align*}
\sum_{n=N+1}^{M}\|B\|_{\mathcal L(X)}^n
\end{align*}
tends to $0$ as $N \to \infty$, uniformly over all $M>N$. Hence $(S_N)$ is a [Cauchy sequence](/page/Cauchy%20Sequence) in $\mathcal L(X)$.
Now we use the Banach-space hypothesis on $X$. Since $X$ is complete, the space $\mathcal L(X)$ of bounded linear maps from $X$ to itself is complete under the operator norm. Thus every Cauchy sequence in $\mathcal L(X)$ converges in operator norm, so there exists an operator $S \in \mathcal L(X)$ such that
\begin{align*}
S_N \to S
\end{align*}
in $\mathcal L(X)$. This limit is precisely the operator-norm sum of the Neumann series:
\begin{align*}
S = \sum_{n=0}^{\infty}(-B)^n.
\end{align*}
[/guided]
[/step]
[step:Identify the Neumann series as the inverse of $I_X+B$]
For each $N \in \mathbb N$, the finite geometric product identity in the associative algebra $\mathcal L(X)$ gives
\begin{align*}
(I_X+B)S_N = I_X-(-B)^{N+1}.
\end{align*}
The same identity also gives
\begin{align*}
S_N(I_X+B) = I_X-(-B)^{N+1}.
\end{align*}
Moreover,
\begin{align*}
\|(-B)^{N+1}\|_{\mathcal L(X)} \le \|B\|_{\mathcal L(X)}^{N+1} \to 0.
\end{align*}
Since multiplication by the fixed bounded operator $I_X+B$ is continuous in the operator norm, passing to the limit in the two displayed identities yields
\begin{align*}
(I_X+B)S = I_X
\end{align*}
and
\begin{align*}
S(I_X+B) = I_X.
\end{align*}
Therefore $I_X+B$ is invertible in $\mathcal L(X)$ and
\begin{align*}
(I_X+B)^{-1}=S=\sum_{n=0}^{\infty}(-B)^n.
\end{align*}
[/step]
[step:Factor $A+E$ and preserve the order of multiplication]
Define $R \in \mathcal L(X)$ by
\begin{align*}
R := SA^{-1}.
\end{align*}
Since $S \in \mathcal L(X)$ and $A^{-1}\in\mathcal L(X)$, the composition $R$ belongs to $\mathcal L(X)$. Using the factorization $A+E=A(I_X+B)$, the identity $S(I_X+B)=I_X$, and the identity $AA^{-1}=I_X$, we obtain
\begin{align*}
R(A+E)=SA^{-1}A(I_X+B)=S(I_X+B)=I_X.
\end{align*}
Using instead $(I_X+B)S=I_X$, we also obtain
\begin{align*}
(A+E)R=A(I_X+B)SA^{-1}=AA^{-1}=I_X.
\end{align*}
Thus $R$ is a two-sided inverse of $A+E$ in $\mathcal L(X)$, so $A+E$ is invertible in $\mathcal L(X)$ and
\begin{align*}
(A+E)^{-1}=R=SA^{-1}.
\end{align*}
Substituting the operator-norm sum for $S$ gives
\begin{align*}
(A+E)^{-1}=\left(\sum_{n=0}^{\infty}(-B)^n\right)A^{-1}.
\end{align*}
Since $B=A^{-1}E$, this is exactly
\begin{align*}
(A+E)^{-1}=\sum_{n=0}^{\infty}(-A^{-1}E)^nA^{-1}.
\end{align*}
The convergence is in operator norm because multiplication on the right by the fixed operator $A^{-1}$ is operator-norm continuous:
\begin{align*}
\| (S_N-S)A^{-1}\|_{\mathcal L(X)} \le \|S_N-S\|_{\mathcal L(X)}\|A^{-1}\|_{\mathcal L(X)} \to 0.
\end{align*}
This proves both the invertibility assertion and the displayed inverse formula.
[/step]