Thanks to Higham I know that $A f(BA) = f(AB) A$ for any two matrices whose sizes are compatible.


Now I believe that $A (BA)^D = (AB)^D A$, even though the Drazin inverse is not the same function (polynomial?) for $AB$ as for $BA$. 


I have validated this relationship via numerical experiments with random matrices, I just can't $prove$ it. 


Can you prove (or disprove) it?



An identity involving the Drazin inverse
[These notes](http://benisrael.net/GI-LECTURE-7.pdf) say that the Drazin Inverse is the matrix function corresponding to $f(z) = 1/z$, defined on the *nonzero* eigenvalues. Thus, by the theorem that you cite, the said equality should hold.


