What can we say about $(I-AD)^{-1}$ if $D$ is a diagonal matrix? Assume we know that square  matrices $A$ and $(I-AD)^{-1}$ are invertible and also $D$ is a diagonal matrix. Also assume that $A$ is a symmetric matrix. My question is when we can express $(I-AD)^{-1}$ as a function of $A, D, A^{-1}, D^{-1}, (I-A)^{-1}$ (without terms $(A-D)^{-1}$ or $(A^{-1}-D)^{-1}$) ? 
For example when matrix A is rank 1, then we have:
$(I-AD)^{-1}=I+\frac{1}{1-tr(AD)} AD$. As you can see if A is rank 1 then we can do this easily. 
The only related paper I found is a paper by Kenneth S. Miller, but it is not useful for higher rank matrices.
I know it might be very hard for general matrix $A$ but can it be done for special cases where for instance matrix A is positive semidefinite? Any comment is highly appreciated.
 A: I have an answer to the following question : Is there a four-variable (not necessarily commutative) polynomial $f$ such that the identity
$$
(I-AD)^{-1}=f((I-A)^{-1},A,D,A^{-1},D^{-1}) \tag{1}
$$
holds, whenever $A$ is symmetric positive definite, and $D$ is invertible and diagonal ?
The answer is NO. Indeed, this is already impossible when $n=2$ and 
$$D=\left(\begin{matrix} 2 & 0 \\ 0 & 3 \end{matrix}\right).$$ If we write
$$
A=\left(\begin{matrix} a & b \\ b & a \end{matrix}\right)
$$
with $a\gt 0$ and $a\gt b$, then
$$
\det(A)=a^2-b^2, \det(I-A)=a^2-2a-b^2+1, \det(I-AD)=6a^2-5a-6b^2+1
$$
Thus the RHS in (1) will always have a denominator of the form
$$
(a^2-b^2)^p (a^2-2a-b^2+1)^q,
$$
where $p$ and $q$ are integers. Now the LHS in (1) will always a denominator of the form 
$(6a^2-5a-6b^2+1)^r$. Since we have three distinct irreducible polynomials in ${\mathbb Q}[a,b]$
here, the denominators will never coincide.
A: In his book Some Eclectic Matrix Theory, Kenneth S. Miller presents the following generalization of your rank-$1$ result. This result is handy if your diagonal matrix $D$ is low rank.

Theorem (Miller, p. 14): Let $H$ be an $n\times n$ real matrix or rank $r$, and let
  $\sigma_1, \ldots, \sigma_n$ be the elementary symmetric functions of
  its eigenvalues. If $I+H$ is nonsingular then $a_r \ne 0$ and
  $$(I+H)^{-1} = I-\frac{1}{a_r}\big(a_{r-1}H - a_{r-2}H^2 + \cdots + (-1)^{r-1}a_0 H^r\big)$$ where $$a_m = 1 + \sum_{1 \le k \le m} \sigma_k, \quad m = 0,1,\ldots,r.$$

As for the elementary symmetric functions of the eigenvalues: They can be computed (via the Newton–Girard identities) from the power sums of the eigenvalues, which themselves are simply the traces of the powers of $H$ (as can be seen from writing $H$ in Jordan form):
$$\operatorname{tr}(H^m) = \sum_{i=1}^n \lambda_i^m, \quad m \ge 0.$$
The upshot is that $(I+H)^{-1}$ is a degree-$\operatorname{rk}(H)$ polynomial in $H$ whose coefficients are rational functions in the traces of the powers of $H$ (Miller, Some Eclectic Matrix Theory, pp. 14-15).
