Eigenvalue perturbation theory via Feynman diagrams Suppose I have a matrix given by a sum
$$A=D+\epsilon B$$
where $D$ is diagonal and $\epsilon$ is small, and I want the eigenvalues of $A$ as a power series in $\epsilon$. The first two orders in perturbation theory are well known. Third and higher orders are briefly discussed here. However, the equations become horrible. 
I hear that Feynman diagrams are an efficient way to formulate perturbation theory, but I can't find an accessible exposition of this approach. Note that I have in mind the simple matrix setting. I don't want vacuum states, quantum field theory, Dirac equation, etc. Can someone help?
 A: Feynman diagrams work really well when considering Green's functions which encapsulate information about eigenvectors and eigenvalues but in an indirect way. If you want these objects themselves, then you need an explicit and therefore necessarily combinatorial/diagrammatic form of the Rayleigh-Schrödinger perturbation series which, in general, is more complicated than the Feynman diagram expansions for Green's functions. Such a formalism is developed in the article "The Rayleigh-Schrödinger perturbation series of quasi-degenerate systems" by Brouder, Duchamp, Patras and Toth. A similar philosophy is followed in the recent work by Imbrie on Anderson localization.
A: I think the simplest way is to use the very simple and very useful resolvant formula
$$ (A+B)^{-1}=A^{-1}-A^{-1}B(A+B)^{-1}.$$
For perturbative theory, we just iterate this formula 
$$ (A+B)^{-1}=A^{-1}-A^{-1}BA^{-1}+A^{-1}BA^{-1}B(A+B)^{-1}$$
and obtain a power series in $B$.
$$ (A+B)^{-1}=A^{-1}-A^{-1}BA^{-1}+A^{-1}BA^{-1}BA^{-1}-A^{-1}BA^{-1}BA^{-1}BA^{-1} +\cdots$$
With your notation it gives 
$$ A^{-1}=D^{-1}-\epsilon D^{-1}BD^{-1}+\epsilon^2 D^{-1}BD^{-1}BD^{-1}-\epsilon^3 D^{-1}BD^{-1}BD^{-1}BD^{-1} +\cdots$$
For the eigenvalues there is then the Cauchy formula
$$ z_0 = \frac{1}{2i\pi}\oint \frac{z}{z-z_0} dz$$
which for matrices gives
$$\lambda_0 = \frac{1}{2i\pi}\oint \lambda Tr(A-\lambda)^{-1}d\lambda$$
where the integral is a closed loop around an eigenvalue of $D$.
And we obtain the perturbative expansion:
$$ \lambda_0 =\frac{1}{2i\pi} \oint \lambda Tr(D-\lambda)^{-1} d\lambda \\-\epsilon \frac{1}{2i\pi} \oint \lambda Tr (D-\lambda)^{-1}B(D-\lambda)^{-1} \\+\epsilon^2 \frac{1}{2i\pi} \oint \lambda Tr (D-\lambda)^{-1}B(D-\lambda)^{-1}B(D-\lambda)^{-1}\\-\epsilon^3 \cdots$$
Now take a loop sufficiently close to the initial eigenvalue of $D$ and it will give what you wanted. 
