I seek for a symmetric positive definite (PD) solution $\mathbf{X}$ for the following quadratic matrix equation:
\begin{array}{cc} \mathbf{XBX=A}+\lambda \mathbf{X} \tag{1}\label{eq1} \end{array}
where $\mathbf{A}$ and $\mathbf{B}$ are 2 given positive definite symmetric matrices and $\lambda>0$ is a given scalar.
Unfortunately, I couldn't find a direct solution for this equation. However, I know that $\mathbf{X}=\mathbf{B}^{-\frac{1}{2}} (\mathbf{{B}^{\frac{1}{2}}A{B}^{\frac{1}{2}}})^{\frac{1}{2}}\mathbf{B}^{-\frac{1}{2}} $ is a symmetric PD solution for the following equation:
\begin{array}{cc} \mathbf{XBX=A} \tag{2}\label{eq2} \end{array}
It seems that I can solve equation \eqref{eq1} using a fixed point iterative algorithm based on equation \eqref{eq2} as follows:
\begin{array}{cc} \mathbf{X}_{k+1} \mathbf{B} \mathbf{X}_{k+1}=\mathbf{A}+\lambda \mathbf{X}_{k} \tag{3}\label{eq3} \end{array}
\begin{array}{cc} \Rightarrow \mathbf{X_{k+1}}=\mathbf{B}^{-\frac{1}{2}} (\mathbf{{B}^{\frac{1}{2}}G_{k}{B}^{\frac{1}{2}}})^{\frac{1}{2}}\mathbf{B}^{-\frac{1}{2}} \tag{4}\label{eq4} \end{array}
where $\mathbf{G}_{k}=\mathbf{A}+\lambda \mathbf{X_{k}}$. I have tested this procedure many times using different simulation situations and in all cases (when $\mathbf{A}$ and $\mathbf{B}$ were PD and $\lambda>0$), this algorithm converged to a solvent of equation \eqref{eq1}.
However, I couldn't prove the convergence of this iterative algorithm. I know that for convergence analysis, I need to seek for an upper bound as follows: \begin{array}{cc} ||\mathbf{X_{k+1}}-\mathbf{X_{k}}||<=\alpha ||\mathbf{X_{k}}-\mathbf{X_{k-1}}|| \end{array} where $\alpha<1$.
Can anyone help me in proving this convergence?
It's not the good way.
Your equation is equivalent to $(BX)^2=BA+\lambda BX$.
Put $Y=BX$; then $Y^2=BA+\lambda Y$, that is,
$(Y-\dfrac{\lambda}{2}I)^2=BA+\dfrac{\lambda^2}{4}I$.
Note that $(\mu_i)=spectrum(BA)\subset (0,+\infty)$ and is diagonalizable. Then we can calculate
$BA=Pdiag(\mu_i)P^{-1}$, $BA+\dfrac{\lambda^2}{4}I=Pdiag(\mu_i+\lambda^2/4)P^{-1}$.
Therefore $Y-\dfrac{\lambda}{2}I=Pdiag(\pm\sqrt{\mu_i+\lambda^2/4})P^{-1}$ ($2^n$ values).
Note that, necessarily, $spectrum(Y)\subset (0,+\infty)$ and $\sqrt{\mu_i+\lambda^2/4}>\lambda/2$. Then
$Y=\dfrac{\lambda}{2}I+Pdiag(\sqrt{\mu_i+\lambda^2/4})P^{-1}$ (only $+$ signum) and
$X=B^{-1}Y$ is the unique $>0$ solution of the considered equation.