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?

  • I just noticed this question of yours. Is it still open? or any further simplifications? (Also in L(X), the log is applies to $G'XG$? right?) – Suvrit Sep 14 '16 at 19:32
  • @Suvrit: Very recently I managed to prove the convergence of the algorithm. As soon as I double-check the proof I will write the (a trace of) solution here. (Yes, the log applies to $G^\top X G$). – Ludwig Sep 14 '16 at 19:49
  • Ah good; I'd be interested to know your solution; I had an idea for a solution but since you've already proved it, I'll just wait for your answer :-) – Suvrit Sep 14 '16 at 19:53

Your Answer


By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.