11
$\begingroup$

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.

$\endgroup$
  • 5
    $\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 Jan 10 '15 at 11:58
13
$\begingroup$

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
$\endgroup$
  • $\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 Jan 10 '15 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$ – DumpsterDoofus Jan 10 '15 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 Jun 25 '15 at 15:51
  • $\begingroup$ Really? I don't remember ever doing that convention except for constants, like $defaultSetting or some such. $\endgroup$ – J. M. will be back soon Nov 6 '15 at 5:09
3
$\begingroup$

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]

enter image description here

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

$\endgroup$

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.