I have the following quadratic matrix equation in symmetric positive semidefinite covariance matrix $\mathbf{C}$
$$\mathbf{C}\mathbf{S}\mathbf{C}+\mathbf{C}+\mathbf{A}=\mathbf{0}$$
where $\mathbf{S}$ is a singular matrix and $\mathbf{A} = \mathbf{a} \mathbf{a}^{T}$. I need to iteratively find $\mathbf{C}$.
1) If I try to set up the problem using Bernoulli iteration, in order to find the dominant solvent I would have to go for
$$ \mathbf{S}+\{\mathbf{I}+\mathbf{C}_{k}^{-1}\mathbf{A}\}*\mathbf{C}_{k+1}^{-1} = \mathbf{0} $$
$$ \Longrightarrow \mathbf{C}_{k+1}^{-1} = \{\mathbf{I}+\mathbf{C}_{k}^{-1}\mathbf{A}\}^{-1}*(-\mathbf{S}) $$
Which won't make sense since $\mathbf{S}$ is singular
2) If I try to to find the minimal solvent, I set it up as
$$ (\mathbf{C}_{k}\mathbf{S}+\mathbf{I})*\mathbf{C}_{k+1}+\mathbf{A}=\mathbf{0} $$
$$ \Longrightarrow \mathbf{C}_{k+1}=(\mathbf{C}_{k}\mathbf{S}+\mathbf{I})^{-1}*(-\mathbf{A}) $$
But this doesn't seem to be working either.
3) A last, non-Bernoulli iteration I set up was
$$ \mathbf{C}_{k}\mathbf{S}\mathbf{C}_{k}+\mathbf{C}_{k+1}+\mathbf{A}=\mathbf{0} $$
$$ \Longrightarrow \mathbf{C}_{k+1}=-\mathbf{A}-\mathbf{C}_{k}\mathbf{S}\mathbf{C}_{k} $$
Which seems to be working but I'm not sure which root it converges to. Also the final estimate $\mathbf{C}_{k}$ converges to, somehow appears to be twice the actual value (upon performing a ground truth test).
My queries are
Am I on the right track on dealing with this problem? If so, where am I making mistakes?
Is there any procedure to find the dominant and minimal solvent when the coefficient of the square term is singular?
Is the last, non-Bernoulli iteration valid? If so, why is the result twice the actual one? Which root does it converge to?
Any guidance will be highly appreciated. Thank you.