I'm looking for a closed form solution to the following equation: $$ \DeclareMathOperator*{\argmax}{arg\,max} \begin{align*} E^* = \argmax_E \ r^T(I - \gamma (A + BE))^{-1}s \end{align*} $$
for fixed matrices $A, B$, vectors $r, s$ and scalar $\gamma$. $E$ is also a matrix, but is the focus of our computation and isn't "fixed" in the same way. As for shapes, $A$ is square, and $BE$ is square, but $B$ and $E$ separately are not square, $r$ and $s$ have the same dimensionality. By construction, $\lVert \gamma (A + BE) \rVert < 1$ needs to be true which I think points the quadratic in the right direction.
A first question is, of course, if $E^*$ even exists. I think it does, but I'd love help proving that.
My initial thoughts were solve for the derivative, set it to zero, and solve for $E$ (I think this derivative is right, but it would be safer to not assume so):
$$ 0 = \gamma B^T ((I - \gamma(A + BE))^{-1})^T r s^T ((I - \gamma(A + BE))^{-1})^T $$
But I started running into problems trying to collect all of the $E$ matrices together when they are embedded in the inverse.
Any ideas?