Difficult matrix equation $A(x-x_0) + \alpha\left(\frac{\Sigma x}{\mu^\top x}\right) = 0$

I can provide some hints for the case where $A,\Sigma$ are strictly positive definite. This assumption guarantees us that both matrices are invertible; let us set henceforth $\lambda=\frac{\alpha}{\mu^T x}$. By multiplying on the right by $\Sigma^{-1}$ the problem transforms into the following one:

$$(B-\lambda I)x=B x_0~~,~~B=\Sigma^{-1}A$$

Of course $B x_0\neq 0$, since $B$ is itself positive definite. Now it is finally obvious that $B-\lambda I$ can only be singular for $\lambda $ in the set of eigenvalues of the matrix $B$ itself, in which case a solution does not exist. For all other values of the parameter a unique solution is guaranteed; if one uses the eigenvalue decomposition of $B=S\Lambda S^{-1}$ the solution is given by the relation

$$x=S(\Lambda-\lambda I)^{-1}S^{-1}Bx_0$$

This solution is naturally, a function of $\lambda$ itself, which can be written in the concise form $$x=\sum_{k}\frac{v_k}{\lambda_k-\lambda}$$ where the sum is performed over the distinct eigenvalues of $B$ with the vectors $(v_k)_i=S_{ik}\sum_{lm}S^{-1}_{kl}B_{lm}(x_0)_{m}$ whenever all the eigenvalues have multiplicity 1 (but with trivial modifications if it's greater than 1-just take the limit of coinciding eigenvalues).

We are not done- $\lambda$ depends on $x$ itself and that imposes a constraint on the possible solutions. Dotting from the left with the vector $\mu$, we obtain the rational equation $$\frac{\alpha}{\lambda}=\sum_{k}\frac{\mu^Tv_k}{\lambda_k-\lambda}$$ which generally has $n$ complex roots. However $\lambda$ is a real parameter. Is a solution guaranteed? It is only guaranteed for odd $n$, since the polynomial equation then is forced to have a real solution. A counterexample to this is given by the function $$f(\lambda)=-\frac{1}{1-\lambda}+\frac{1}{2-\lambda}-\frac{3}{\lambda}$$ which has no real roots. However, it is true that there may be found a large enough value of $|\alpha|$ for which a real solution exists, regardless of the values of $n$ and the projections $\mu^T v_k$.

EDIT: I believe that this is as close algorithmically as one can get to a "closed form". The problem has certainly been reduced to a simpler one, but the impasse is obvious- we cannot solve arbitrary degree polynomial equations in closed form.


By multiplying with $(\mu^T x)$ we see that the equation is of the form of a system of second degree polynomials. These systems are generally not solvable "in closed form", since (1) higher degree polynomial systems can be reduced to second degree systems in more variables (2) polynomials of degree $\ge 5$ cannot be solved "in closed form". (cf. discussion on mathoverflow).

So you shouldn't expect a "closed form" solution. Depending on your use case, you can try two things:

  1. Try solving particular instantiations of the systems using computer algebra packages such as sympy

  2. Use Numerical methods such as Newton's method or fixed point iteration

At the end of the day, the question you should always ask yourself is whether do you really need a "closed form" solution.