I have the following quadratic matrix equation:

$$ XAX+X = B $$

where both $A$ and $B$ are given positive definite matrices, and $X$ is a covariance matrix and, hence, positive definite.

When there is no constraint, the equation can be solved via Bernoulli iteration in the following form:

$$X_{k+1} = -A^{-1}(I-BX_k^{-1})$$

However, this does not seems to preserve positive semidefinite.

Any guidance would be appreciated. Thank you.