Let's assume that $Q_0$ is $3\times 3$ matrix with $\det Q_0\neq0 $ and

$$ Q_{i+1} = \frac{1}{2}\left[ Q_i+(Q_i^{-1})^T \right] $$

I need to find next limit: $\lim _{i \to +\infty}$$Q_i$. In other words I need to find such $Q_i$ that $\left|Q_{i+1} - Q_i\right| < \epsilon$. Does anybody know, is it possible to find such kind of limits with Mathematica?

  • $\begingroup$ I don't know the answer, but do you need more assumptions on the determinant at each iteration ? $\endgroup$ Oct 30 '13 at 19:24
  • $\begingroup$ @b.gatessucks that just says that $Q_0$ is invertible. The trick is then showing that $\det (Q_1 + (Q^{-1}_1)^T) \neq 0$. $\endgroup$
    – rcollyer
    Oct 30 '13 at 19:41
  • 1
    $\begingroup$ To see the limit appearing for random matrices one can try n = 3; Q0= RandomReal[{-100, 100}, {n, n}]; res = NestWhileList[(# + Transpose@Inverse@#)/2 &, Q0, Norm[#1 - #2] > 10^-8 &, 5, 10^4]; limit = Last@res ListAnimate[MatrixPlot[#] & /@ res] $\endgroup$ Oct 30 '13 at 20:30
  • $\begingroup$ q0={{-1, -1, 1/5 (-2 - 4 I Sqrt[6])}, {1, -1, -1}, {0, 1, 1}} is invertible but q1 = 1/2 (q0 + Transpose[Inverse[q0]]) is not. $\endgroup$ Oct 31 '13 at 7:11

I don't know how to prove that starting from an invertible matrix each iteration will give an invertible matrix. Also, I'm going to consider real matrices only (because of what I have in the comments).

We can start by the singular value decompostion of q0 = Dot[u, w0, Transpose[v]] where w is diagonal. Then you can easily convince yourself that the same decomposition holds at each iteration, i.e. q1 = Dot[u, w1, Transpose[v]] where w1 = 1/2 (w0 + Transpose[Inverse[w0]]) and so on. All matrices in the last equation are diagonal, so the equality holds equivalently for each diagonal element.

The solution to that equation is given by :

RSolve[{a[i + 1] == 1/2 (a[i] + 1/a[i]), a[0] == a0}, a[i], i]
(* {{a[i] -> Coth[2^i ArcCoth[a0]]}} *)

and you can then take the limit from here.

Check :

q0 = RandomReal[UniformDistribution[], {3, 3}];
svd = SingularValueDecomposition[q0];

Nest[1/2 (# + Transpose[Inverse[#]]) &, q0, 10] == 
      DiagonalMatrix[Coth[2^(10) ArcCoth[#]] & /@ Diagonal[svd[[2]]]],            
  ] // Chop)
(* True *)
  • $\begingroup$ In general starting from an invertible need not always provide another. Take {{0,1},{-1,0}} for example, as its inverse is also its negation. In general I believe that when the sequence stays in the invertible family, the convergence is to a diagonal matrix with +-1 in each position on the main diagonal. $\endgroup$ Nov 2 '13 at 16:49
  • $\begingroup$ Thanks, this code seems to work well, however as I understand in order to obtain $Q$ matrix I always need to perform SVD. SVD works more slowly than this iteration formula which actually perform Polar Decomposition. $\endgroup$ Nov 2 '13 at 19:37
  • $\begingroup$ You need to do it only once for each matrix Q0 you need to consider. $\endgroup$ Nov 2 '13 at 19:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.