Editing Polar decomposition (section)

==Numerical determination of the matrix polar decomposition==
To compute an approximation of the polar decomposition ''A'' = ''UP'', usually the unitary factor ''U'' is approximated.<ref name="higham1986" /><ref name="byers2008" /> The iteration is based on [[Heron's method]] for the square root of ''1'' and computes, starting from <math>U_0 = A</math>, the sequence
<math display="block">U_{k+1} = \frac{1}{2}\left(U_k + \left(U_k^*\right )^{-1}\right),\qquad k = 0, 1, 2, \ldots</math>

The combination of inversion and Hermite conjugation is chosen so that in the singular value decomposition, the unitary factors remain the same and the iteration reduces to Heron's method on the singular values.

This basic iteration may be refined to speed up the process:
{{bulleted list
| Every step or in regular intervals, the range of the singular values of <math>U_k</math> is estimated and then the matrix is rescaled to <math>\gamma_k U_k</math> to center the singular values around ''1''. The scaling factor <math>\gamma_k</math> is computed using matrix norms of the matrix and its inverse. Examples of such scale estimates are:
<math display="block">\gamma_k = \sqrt[4]{\frac{\left\|U_k^{-1}\right\|_1 \left\|U_k^{-1}\right\|_\infty}{\left\|U_k\right\|_1 \left\|U_k\right\|_\infty} }</math>
using the row-sum and column-sum [[induced norm|matrix norms]] or
<math display="block">\gamma_k = \sqrt{\frac{\left\|U_k^{-1}\right\|_F}{\left\|U_k\right\|_F}}</math>
using the [[Frobenius norm]]. Including the scale factor, the iteration is now
<math display="block">U_{k+1} = \frac{1}{2}\left(\gamma_k U_k + \frac{1}{\gamma_k} \left(U_k^*\right)^{-1}\right), \qquad k = 0, 1, 2, \ldots</math>
| The [[QR decomposition]] can be used in a preparation step to reduce a singular matrix ''A'' to a smaller regular matrix, and inside every step to speed up the computation of the inverse.
| Heron's method for computing roots of <math>x^2 - 1 = 0</math> can be replaced by higher order methods, for instance based on [[Halley's method]] of third order, resulting in
<math display="block">U_{k+1} = U_k\left(I + 3U_k^* U_k\right)^{-1}\left(3I + U_k^* U_k\right),\qquad k = 0, 1, 2,  \ldots</math>

This iteration can again be combined with rescaling. This particular formula has the benefit that it is also applicable to singular or rectangular matrices ''A''.
}}