For given $n\times n$ real symmetric positive definite matrices $X_{1}, X_{2}$, let $Y := X_{1}^{-1}X_{2}$, and let $I$ be the identity matrix.
I would like to solve the following equation for the unknown $p>0$, in terms of functions of the matrices $X_{1}$ and $X_{2}$:
$$\text{trace}\left[\left(I + pY\right)^{-1} \left(I - p^{2}Y\right)\right] = 0.$$
I tried expanding the inverse in Neumann series but I am not sure how that helps to solve for $p$. This is not homework.
Write $Y=SDS^{-1}$, where $D$ is diagonal (since $X_1$ and $X_2$ are psd, $Y$ is diagonalizable). Then, observe that
\begin{equation*} f(p) = \text{tr}(S^{-1}(I+pSDS^{-1})^{-1}SS^{-1}(I-p^2SDS^{-1})S). \end{equation*} This simplifies into a simpler function of $p$ involving only the diagonal values in $D$. That particular equation you can solve using Newton's method.
EDIT: as noted in my comments below, the same idea also works more generally if we use the Schur decomposition to reduce matrices to an upper triangular form; we do not need to assume diagonalizability in this case.
Given $\mathrm A \in \mathbb R^{n \times n}$, we define $f : \mathbb R \to \mathbb R$ as follows
$$f (x) := \mbox{tr} \left( (\mathrm I_n + x \mathrm A)^{-1} (\mathrm I_n - x^2 \mathrm A) \right)$$
It is not necessary that $\mathrm A$ be diagonalizable. Using the Schur decomposition $\mathrm A = \mathrm Q \mathrm U \mathrm Q^{\top}$, where $\mathrm Q$ is orthogonal and $\mathrm U$ is upper triangular with the eigenvalues of $\mathrm A$ on its main diagonal, we have
$$\begin{array}{rl} f (x) &= \mbox{tr} \left( (\mathrm I_n + x \mathrm A)^{-1} (\mathrm I_n - x^2 \mathrm A) \right)\\ &= \mbox{tr} \left( (\mathrm Q \mathrm Q^{\top} + x \, \mathrm Q \mathrm U \mathrm Q^{\top})^{-1} (\mathrm Q \mathrm Q^{\top} - x^2 \, \mathrm Q \mathrm U \mathrm Q^{\top}) \right)\\ &= \mbox{tr} \left( \mathrm Q (\mathrm I_n + x \, \mathrm U )^{-1} \mathrm Q^{\top} \mathrm Q (\mathrm I_n - x^2 \, \mathrm U) \mathrm Q^{\top} \right)\\ &= \mbox{tr} \left( (\mathrm I_n + x \, \mathrm U)^{-1}(\mathrm I_n - x^2 \, \mathrm U) \right)\end{array}$$
where both $\mathrm I_n + x \, \mathrm U$ and $\mathrm I_n - x^2 \, \mathrm U$ are upper triangular. Note that
The inverse of an upper triangular matrix is also upper triangular. Since $\mathrm I_n + x \, \mathrm U$ is upper triangular, $(\mathrm I_n + x \, \mathrm U)^{-1}$ is also upper triangular.
The product of two upper triangular matrices is also upper triangular with the main diagonal being the entrywise product of the main diagonals of the factors. Hence, we can conclude that $(\mathrm I_n + x \, \mathrm U)^{-1}(\mathrm I_n - x^2 \, \mathrm U)$ is also upper triangular.
Suppose that $\mathrm A$ has $m$ distinct eigenvalues with multiplicities $n_1, n_2, \dots, n_m$. Hence,
$$f (x) = \sum_{k=1}^m n_k \left( \dfrac{1 - x^2 \, \lambda_k (\mathrm A)}{1 + x \, \lambda_k (\mathrm A)} \right)$$
Thus, the equation $f (x) = 0$ can be written as follows
$$n_1 \left( \dfrac{1 - x^2 \, \lambda_1 (\mathrm A)}{1 + x \, \lambda_1 (\mathrm A)} \right) + n_2 \left( \dfrac{1 - x^2 \, \lambda_2 (\mathrm A)}{1 + x \, \lambda_2 (\mathrm A)} \right) + \cdots + n_m \left( \dfrac{1 - x^2 \, \lambda_m (\mathrm A)}{1 + x \, \lambda_m (\mathrm A)} \right) = 0$$
Suppose we have
$$\mathrm A = \begin{bmatrix} 1 & 0 & 0\\ 0 & 2 & 1\\ 0 & 0 & 2\end{bmatrix}$$
which is in Jordan normal form and, thus, is already upper triangular. Matrix $\mathrm A$ has $m = 2$ distinct eigenvalues with multiplicities $n_1 = 1$ and $n_2 = 2$. Using SymPy:
>>> from sympy import *
>>> A = Matrix([[1,0,0],[0,2,1],[0,0,2]])
>>> x = Symbol('x')
>>> ((eye(3) + x * A)**-1 * (eye(3) - x**2 * A)).trace()
2*(-2*x**2 + 1)/(2*x + 1) + (-x**2 + 1)/(x + 1)
We obtain the equation
$$\frac{1 - x^{2}}{x + 1} + 2 \left(\frac{1 - 2 x^{2}}{2 x + 1}\right) = 0$$
Solving the equation,
>>> solve(2*(-2*x**2 + 1)/(2*x + 1) + (-x**2 + 1)/(x + 1))
[1/12 + sqrt(73)/12, -sqrt(73)/12 + 1/12]
