Solve matrix equation $X = I + A + AXA$ Notations
$X$ : p times p positive semi definite, unknown.
$I$ : p times p identity matrix, known.
$A$ : p times p positive semi definite, known.

I want to find X, but it is difficult to me.
Thus, instead of exact solution, I used following iterative method to find the approximate solution.
Step 1) : initialize $X(0)$
Step 2) : $X(k) = I + A + A X(k-1) A$
Step 3) : until $|X(k) - X(k-1)|_F < tol$
But, this method highly depends on the initial value X(0) and is vulnerable to diverge.
Please modify my method or propose another method.
Regard, Han
 A: From induction we have $$X=I+A+\cdots +A^{2n-1}+A^nXA^n$$by tending $n\to \infty$ and defining $B=\sum_{n=0}^\infty A^n$, we have $B=(I-A)^{-1}$ only if $I-A$ is invertible, unless $X$ does not exist. Also all the eigenvalues of $A$ must have $|\cdot |<1$. Therefore $A^n\to 0$ and $$X=(I-A)^{-1}$$
A: I suggest you to follow this for a strict justification that $X=\left(I-A\right)^{-1}$.
First it is easy to see that each eigenvalue $\lambda$ of $A$ satisfies $0\le\lambda <1$ for the given equation to have a solution. To illustrate, let $Av = \lambda v$ for some $v\in\mathbb{R}^p\setminus\left\{0\right\}$. Then the given equation yields
\begin{align*}
v^t Xv &= v^t v + v^t Av + v^t AXAv
\\&= v^t v + \lambda v^t v + \lambda^2 v^t Xv,
\end{align*}
or,
\begin{align*}
\left(1-\lambda^2 \right)v^t Xv = \left(1+\lambda\right)v^t v.
\end{align*}
If $\lambda\ge1$, as $v^t v>0$, the right hand side must be positive. While the positive semidefiniteness of $X$ gives $v^t Xv \ge0$ and hence the left hand side is $\le 0$, yielding a contradiction. This together with the positive semidefiniteness of $A$ proves the claim.
Now multiplying $\left(I-A\right)$ to the given equation by left, it becomes
$$
\left(I-A\right)X = I-A^2 + A \left(I-A\right)XA
$$
or we can write
$$
I-\left(I-A\right)X = A\left(I-\left(I-A\right)X\right)A.
$$
Replacing $I-\left(I-A\right)X$ with $Y$, we arrive at a simpler form
$$
Y=AYA.
$$
Let $v$ and $w$ be eigenvectors of $A$ corresponding to the eigenvalues $\lambda$ and $\mu$. Then
\begin{align*}
v^t Y w = v^t AYA w = \lambda\mu\left( v^t Y w\right)
\end{align*}
holds, which together with the fact that $0\le\lambda, \mu<1$ force $v^t Y w =0$. Note, however, that $A$ is symmetric and thus is diagonalizable, that is, its eigenvectors span any vectors in $\mathbb{R}^p$. This implies $v^t Y w =0$ for any $v,w\in\mathbb{R}^p$ which holds only if $Y=O$. From this we have $\left(I-A\right)X = I$ and therefore $X=\left(I-A\right)^{-1}$.
A: More generally, let $K$ be an algebraically closed field and $A\in M_n(K)$ with $spectrum(A)=(\lambda_i)_i$.
We consider the equation in $X\in M_n(K)$: $(*)$ $X=I+A+AXA$.
$\textbf{Proposition}$. If  for every $(i,j),$ $\lambda_i\lambda_j\not= 1$ then $(*)$ has the sole solution 
$X=(I-A)^{-1}$. 
$\textbf{Proof}$. We use the Ramanasa's trick. If $Y=I-(I-A)X$, then $f(Y)=AYA-Y=0$. $f$ is linear and $spectrum(f)=(\lambda_i\lambda_j-1)_{i,j}$.
(cf. https://en.wikipedia.org/wiki/Kronecker_product).
Thus $0\notin spectrum(f)$, $f$ is one to one and $Y=0$. Since $1\notin spectrum(A)$, we deduce that $X=(I-A)^{-1}$. $\square$
