Can we find the scalar $\lambda$ in the following quadratic vector-matrix form? Assume $\lambda$ and $\rho$ are scalers, $b \in \mathbb{R}^d$, $I$ is the identity $d \times d$ matrix and $D$ is a $d \times d$ diagonal positive semidefinite matrix. Is there a way to calculate $\lambda$ from the following quadratic equation?
$$\frac 14 b^T(\lambda I - D)^{-2}b = \rho^2$$
 A: $
\def\d{\delta}
\def\l{\lambda}
\def\o{{\tt1}}
\def\a{\approx}
\def\p{d}
\def\LR#1{\left(#1\right)}
\def\trace#1{\operatorname{trace}\left(#1\right)}
\def\qiq{\quad\implies\quad}
\def\grad#1#2{\frac{\p #1}{\p #2}}
\def\c#1{\color{red}{#1}}
\def\CLR#1{\c{\LR{#1}}}
$For typing convenience, define the matrix
$$\eqalign{
M &= \LR{\l I-D}^{-1} \qiq \c{dM=-M^2\,d\l} \\
}$$
Set up the nonlinear function and calculate its derivative
$$\eqalign{
f(\l) &= 4\rho^2 - b^TM^2\,b \\
df &= -b^T\LR{2M\,\c{dM}}b \\
  &= -b^T\LR{2M\CLR{-M^2\,d\l}}b \\
  &= 2\,b^TM^3b\:d\l  \\
f'=\grad{f}{\l} &= 2\,b^TM^3\,b \\
}$$
Now apply Newton's method to find the roots of $f(\l)=0,\,$ i.e.
$$\eqalign{
\l_+
 &= \l - \frac{f}{f'} \\
 &= \l - \frac{4\rho^2 - b^TM^2\,b}{2\,b^TM^3\,b} \\
 &= \l + \frac{b^T\LR{\l I-D}^{-2}b\:-\:4\rho^2}{2\,b^T\LR{\l I-D}^{-3}b} \\
}$$
This function isn't particularly tricky, so the iterations should converge rapidly. However, since there are lots of roots, different initial guesses will converge to different solutions.
Approximate Solutions
Replace $D$ by its isotropic approximation
$$D\a\d I\qiq \d=\frac{\trace D}{\trace I}$$
This reduces the matrix equation to a scalar equation which is easily solved
$$\eqalign{
4\rho^2 &\a \|b\|^2\:\LR{\l-\d}^{-2} \\
\l &\a \d + \frac{\|b\|}{2\rho} \\
}$$
One Newton iteration yields a better closed-form approximation
$$\eqalign{
\mu &= \d + \frac{\|b\|}{2\rho}, \qquad M = \LR{\mu I-D}^{-1} \\
\l &\a \frac{b^T\!\LR{M^2+2\mu M^3}b\:-\:4\rho^2}{2\,b^TM^3b} \\
}$$
A: Because $(\lambda I -D)^{-2}(\lambda I -D)^{-1} = I $
and   $(\lambda I -D)(\lambda I -D)^{-1} = I $,
Then  $(\lambda I -D)^{-2} = (\lambda I -D) $
The above equation becomes  $b^T(\lambda I -D)b=(2 \rho)^2$
Then,
$\Sigma (\lambda - d_i)b_i^2 = (2 \rho)^2$
$\Sigma \lambda b_i^2 - \Sigma d_i b_i^2 = (2 \rho)^2$
$\Sigma \lambda b_i^2   = (2 \rho)^2 + \Sigma d_i b_i^2$
$\lambda = \frac{(2 \rho)^2 + \Sigma d_i b_i^2}{\Sigma b_i^2}$
A: You have
$ b^T (\lambda I - D)^{-2} b = 4 \rho^2 $
The left hand side expands as follows
$ b^T (\lambda I - D)^{-2} b = \displaystyle \sum_{i = 1 }^d \dfrac{b_i^2}{ (\lambda - D_{ii})^2} $
Define the $2d$-th polynomial
$ f( \lambda) = \displaystyle \prod_{i=1}^d (\lambda - D_{ii})^2 $
And the $2(d-1)$-th polynomials
$ f_j (\lambda ) = \displaystyle \prod_{i=1, i \ne j}^d (\lambda - D_{ii})^2 $
Then
$\displaystyle \sum_{j = 1}^d b_j^2 f_j(\lambda) = 4 \rho^2 f(\lambda) $
And this a polynomial equation of degree $2 d$, whose roots are the required $\lambda$'s.
