Revision history [back]

Here is a way to do this for some matrix with exact rational coefficients. The computations are done in the field $\overline{\Bbb Q}$. Here is the sample matrix:

import random
random.seed(r'https://ask.sagemath.org/question/76804/polar-decomposition-of-a-matrix')

A = matrix(QQbar, 3, 3, [random.choice([1..100]) for j in range(9)])

So it is:

sage: A
[68 30 10]
[99 28 56]
[15 76  7]
sage: A.det()
-200630

We build the positive matrix $B = A^*A$:

sage: B = A.conjugate_transpose() * A
sage: B.is_positive_definite()
True
sage: B.eigenvalues()
[449.9840259214669?, 4341.671055858722?, 20603.34491821981?]

We need a square root of $B$, call it $P$. Then $U:=AP^{-1}$ is unitary, and we are done:

A.conjugate_transpose() * A
D, S = B.base_extend(QQbar).jordan_form(transformation=True)
T = S^-1    # so S*D*T == B
P = S * diagonal_matrix([sqrt(D[j, j]) for j in range(D.nrows())]) * T
U = A * P^-1

And with the above code we have:

sage: U.is_unitary()
True
sage: A == U*P
True
sage: U
[ 0.788723890595904?  0.163845291028464?  -0.59251105053920?]
[ 0.596245514676019?   0.03079691068744?   0.80221121690041?]
[-0.149686040204124?  0.986005208353337?  0.073401760660776?]
sage: P
[110.41623991038545?  28.98045207329225?  40.22918544638722?]
[ 28.98045207329225?   80.7140680649559?  10.26511636725468?]
[ 40.22918544638722?  10.26511636725468?  39.51252996565655?]

Note: I was working over QQbar in order to be able to test unitarity (orthogonality) exactly. Just work over some inexact field (like RR) if this is not an issue.