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}
$$
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 to 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.