# Polar Decomposition

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.

• 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 – Nasser Jan 10 '15 at 11:58

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}]];
UnitaryMatrixQ[U]
HermitianMatrixQ[P]
PositiveSemidefiniteMatrixQ[P]


producing

True
True
True

• +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. – Mr.Wizard Jan 10 '15 at 15:25 • @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. – DumpsterDoofus Jan 10 '15 at 15:38
• I tried the code and for my random values, HermitianMatrixQ[P] return False. But P - Conjugate[Transpose[P]] // Chop returned only zeros. – anderstood Jun 25 '15 at 15:51
• Really? I don't remember ever doing that convention except for constants, like $defaultSetting or some such. – J. M. will be back soon Nov 6 '15 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 ClearAll["Global*"] SeedRandom[12345]; 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];
MatrixPlot[Chop[ConjugateTranspose[U].U]]
UnitaryMatrixQ[U]
PositiveSemidefiniteMatrixQ[H]

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

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