I’m stuck in a problem that concerns a nonlinear iterative matrix algorithm. Although the problem is quite complicated to explain I’ll try to describe it in a simple way, neglecting unnecessary details.

The matrix iteration is the following one: $$ X_{k+1} = \int_{-\pi}^{\pi} X_k^{1/2}G(e^{j\omega})\frac{\Psi(e^{j\omega})}{G^\top(e^{-j\omega})X_k G(e^{j\omega})}G^\top(e^{-j\omega})X_k^{1/2}\frac{\mathrm{d}\omega}{2\pi}\quad\quad\quad (\star) $$ where $\{X_k\}$ is a sequence of $n\times n$ matrices, and $G(e^{j\omega})$, $\Psi(e^{j\omega})=\Psi^\top(e^{-j\omega})$ are $n\times 1$ and $1\times 1$, respectively, matrix-valued functions analytic on the complex unit circle which satisfies certain conditions so that there exists at least one fixed point of the iterative algorithm (actually a set of fixed points) and $\int_{-\pi}^{\pi} \Psi(e^{j\omega})\frac{\mathrm{d}\omega}{2\pi}=1$. Furthermore it can be proved that the above iteration maps positive-definite matrices into positive-definite matrices and if we initialize the algorithm using a trace-one matrix the algorithm is trace preserving.

Experimentally, I noticed that the iteration converges to a positive semi-definite matrix for every positive-definite initialization of the algorithm. Moreover, from the simulations, it seems that the iteration is a contraction in the Hilbert (and Thompson) metric. Hence, my hope is that there exists a way to formally prove this contractive property. I recall the the Hilbert metric on the cone of positive definite matrices is defined as $$ d_H(X,Y) = \log \frac{\lambda_{\mathrm{max}}(Y^{-1/2}XY^{-1/2})}{\lambda_{\mathrm{min}}(Y^{-1/2}XY^{-1/2})}, $$ where $X$ and $Y$ are positive definite matrices and $\lambda_{\mathrm{max}}(\cdot)$, $\lambda_{\mathrm{min}}(\cdot)$ denote the maximum and minimum eigenvalue, respectively. By plugging the latter expression into the initial iteration, we get $$ d_H(X_{k+1},X_k) = \log \frac{\lambda_{\mathrm{max}}\left(\int_{-\pi}^{\pi} G(e^{j\omega})\frac{\Psi(e^{j\omega})}{G^\top(e^{-j\omega})X_k G(e^{j\omega})}G^\top(e^{-j\omega})\frac{\mathrm{d}\omega}{2\pi}\right)}{\lambda_{\mathrm{min}}\left(\int_{-\pi}^{\pi} G(e^{j\omega})\frac{\Psi(e^{j\omega})}{G^\top(e^{-j\omega})X_k G(e^{j\omega})}G^\top(e^{-j\omega})\frac{\mathrm{d}\omega}{2\pi}\right)}. $$ Starting from this equation the only thing that I was able to prove, using trivial matrix facts, is that $$ d_H(X_{k+1},X_k)\leq 2d_H(X_{k},X_{k-1}) $$ but clearly this bound is not useful at all.

I’m sorry for the fuzziness of my question, but I have no ideas about how to tackle this problem. Every suggestion or comment is highly appreciated. Thank you.

EDIT. Some further considerations. Assuming that $\Psi(e^{j\omega})=1$ for all $\omega\in[-\pi,\pi)$ and $G(e^{j\omega})=(e^{j\omega}I_n-A)^{-1}B$ where $A$ is $n\times n$ Schur stable (i.e. spectrum lies inside the unit circle) and $B$ is $n\times 1$ and such that the pair $(A,B)$ is controllable, the above iteration can be written as a Lyapunov-like matrix equation $$ X_{k+1} = X_k^{1/2}Z_k X_k^{-1/2}X_{k+1}X_k^{-1/2}Z_k X_k^{1/2}+X_k^{1/2}\frac{BB^*}{B^*P_kB}X_k^{1/2}, $$ where $P_k$ is the unique stabilizing solution of the following discrete-time algebraic Riccati equation $$ \Pi = A^*\Pi A-A^*\Pi B(B^*\Pi B)^{-1}B^*\Pi A+X_k $$ and $Z_k:=A-B(B^*P_k B)^{-1}B^*P_k A$ is Schur stable.

Now it remains to prove that the above matrix equation converges to a positive semi-definite matrix for every positive definite initialization.

EDIT #2. Additional remarks. Another way to prove convergence is to find a suitable Lyapunov function and prove that this function is decreasing along the trajectories generated by the iteration. Consider the candidate Lyapunov function $$ L(X)=\mathrm{trace}(X)-\int_{-\pi}^{\pi} \Psi(e^{j\omega}) \log G^\top(e^{-j\omega})X G(e^{j\omega}) \frac{\mathrm{d}\omega}{2\pi}. $$ Extensive simulations suggest that this is a valid Lyapunov function. A good deal more is true. If we rewrite $(\star)$ as $$ X_{k+1}=X_k+\Delta X_k $$ where $\Delta X_k = X_{k+1}-X_k$, it can be proved that $$ \langle \nabla L(X_k), \Delta X_K\rangle =\mathrm{trace}(\nabla L(X_k) \Delta X_K)\leq 0, $$ i.e. the iteration $(\star)$ can be seen as a gradient descent algorithm with fixed step size equal to 1. Obviously, this is not sufficient to prove convergence, but can it be useful somehow?