I'm new to Mathematica 10 (half-a-month user). I'm looking for a function of polar decomposition but for now I cannot find it. Can anyone tell me about related functions. Thanks in advance.

  • 6
    $\begingroup$ from tutorial/LinearAlgebra/MatrixManipulation.html, it says "Compute the polar decomposition of a matrix:" {u, w, v} = SingularValueDecomposition[{{1, 1}, {1, I}}]; {u.ConjugateTranspose[v], v.w.ConjugateTranspose[v]} // RootReduce $\endgroup$
    – Nasser
    Commented Jan 10, 2015 at 11:58

2 Answers 2


As Nasser pointed out, the polar decomposition of a matrix is not pre-built into Mathematica, but it can be easily computed from the singular value decomposition. For reference, here is a method to do it:

polarDecomposition[m_] := {#.#3\[ConjugateTranspose], #3.#2.#3\[ConjugateTranspose]} & @@ SingularValueDecomposition[m];

This tests that it actually works:

{U, P} = polarDecomposition[RandomComplex[{-1 - I, 1 + I}, {10, 10}]];


  • $\begingroup$ +1 but I suggest you rename $PolarDecomposition. That naming convention is used for system parameters like $Path or $Post, and IMHO it should be reserved for such. $\endgroup$
    – Mr.Wizard
    Commented Jan 10, 2015 at 15:25
  • $\begingroup$ @Mr.Wizard: Yeah, I know that J.M. tried to promote the use of $ as a prefix for user-defined functions, but it does seem to clash with system stuff. $\endgroup$ Commented Jan 10, 2015 at 15:38
  • $\begingroup$ I tried the code and for my random values, HermitianMatrixQ[P] return False. But P - Conjugate[Transpose[P]] // Chop returned only zeros. $\endgroup$
    – anderstood
    Commented Jun 25, 2015 at 15:51
  • $\begingroup$ Really? I don't remember ever doing that convention except for constants, like $defaultSetting or some such. $\endgroup$ Commented Nov 6, 2015 at 5:09

In the polar decomposition of the rectangular complex matrix $A_{m\times n}=UH$, ($m\geq n$), the unitary polar factor $U$ is unique if $A$ is non-singular. On the other hand, it is enough to compute $U$ since $H=U^*A$. Therefore, one approach is to apply iterative methods specially when working in high precision computing.

In what follows, I provide the Newton's iteration as simply as possible to calculate the unitary polar factor $$U_{k+1}=\frac{1}{2}\left(U_k+U_k^{-*}\right),$$ wherein $U_0=A$, $U_k^{-*}=U_k^{{-}^{*}}$, $U_k^{-}$ is the pseudo-inverse, see e.g. Chapter 8 of [N.J. Higham, Functions of Matrices: Theory and Computation, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2008.]

This can be done as

m = 1010; n = 1000;
A = RandomComplex[{-10. - I, 10. + I}, {m, n}];
Id = SparseArray[{i_, i_} -> 1., {n, n}];
max = 25;

k = 0; tolerance = 10^-6; R[0] = 1; U = A;
While[k < max && R[k] >= tolerance,
   {u = U; U1 = PseudoInverse[U]; U2 = ConjugateTranspose[U1];
    U = 1/2 (U + U2); 
    R[k + 1] = Norm[U - u, Infinity]/Norm[U, Infinity]; 
    k++;}]; // AbsoluteTiming

In the above code, I considered solving a complex random $1010\times 1000$ example, with the tolerance $10^{-6}$ while the maximum number of iterates set to 25.

To check the results, it is easy to observe that

H = Chop[ConjugateTranspose[U].A];

set = ScientificForm@Abs@N[Table[R[k], {k, 1, k}], 3]
Grid[{{"Number of iterations", (k)}, {"Stop termination", 
   ScientificForm@N[R[k], 3]},
   N[Log[Abs[R[k]]/Abs[R[k - 1]]]/Log[Abs[R[k - 1]]/Abs[R[k - 2]]], 
    3]}}, Alignment -> Left, Frame -> All]

enter image description here

Here, COC is an approximate calculation of the rate of convergence which should be around 2.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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