pythonnumpysvddecomposition

Finding a polar decomposition of a matrix


I am trying to find the polar decomposition of a matrix. I tried studying it from a video lecture : https://www.youtube.com/watch?v=VjBEFOwHJoo I tried implementing the way they have done, trusting their calculations of Eigen values and Eigen vectors as shown below and it gives correct value of U and P matrix.

A = np.array([[4,-14,1],[22,23,18],[15,10,20]])
A_T_multiply_A = np.array(A.T @ A)
print("A transpose into A")
print(A_T_multiply_A)


x = np.array([[-1,1,1],[0,-2,1],[1,1,1]])           # Eigen vectors
y = np.array([[25,0,0],[0,225,0],[0,0,2025]])      # Eigen values
z = np.linalg.inv(np.array([[-1,1,1],[0,-2,1],[1,1,1]]))    # Inverse Eigen vectors
A_transpose_A = x @ y @ z
print("A transpose A by multiplying Eigen vectors, Eigen values and Inv of Eigen vectors")
print(A_transpose_A)  
P  = x @ np.sqrt(y) @ z
print("Matrix P")
print(P)
U = A @ np.linalg.inv(P)
print("Matrix U")
print(U)

But when I try to implement the Eigen vectors and Eigen values using the API from numpy, the thing do not match:

import numpy as np

# Define the matrix
A = np.array([[4,-14,1],
              [22,23,18],
              [15,10,20]])
A_T_multiply_A = np.array(A.T @ A)

# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A_T_multiply_A)

# Print the results
print("Eigenvalues:")
print(eigenvalues)

print("\nEigenvectors:")
print(eigenvectors)

Not sure why I am getting this difference? It would be great if someone can explain me this or provide me some good link.

The values of Eigen values and Eigen vectors I got using numpy :

   Eigenvalues:
[2025.   25.  225.]

Eigenvectors:
[[-5.77350269e-01 -7.07106781e-01  4.08248290e-01]
 [-5.77350269e-01 -3.91901695e-16 -8.16496581e-01]
 [-5.77350269e-01  7.07106781e-01  4.08248290e-01]]

Also, providing the python code for calculating polar decomposition using a numpy API.

    import numpy as np
    from scipy.linalg import polar
    A = np.array([[4,-14,1],[22,23,18],[15,10,20]])
    U, P = polar(A)
    print("Matrix U=")
    print(U)
    print("Matrix P=")
    print(P)

The result is:

Matrix U=
[[ 6.00000000e-01 -8.00000000e-01  2.40023768e-16]
 [ 8.00000000e-01  6.00000000e-01  3.21346710e-16]
 [ 2.32191141e-16  6.61562711e-17  1.00000000e+00]]
Matrix P=
[[20. 10. 15.]
 [10. 25. 10.]
 [15. 10. 20.]]

Solution

  • Here is a solution using svd:

    import numpy as np
    
    vecU, vals, vecV = np.linalg.svd(A)
    
    P = vecV.T @ np.diag(vals) @ vecV
    U = vecU @ vecV
    P
    array([[20., 10., 15.],
           [10., 25., 10.],
           [15., 10., 20.]])
    U
    array([[ 6.00000000e-01, -8.00000000e-01,  5.73090676e-16],
           [ 8.00000000e-01,  6.00000000e-01,  3.76857861e-16],
           [ 2.32191141e-16, -5.92653245e-17,  1.00000000e+00]])
    

    Compare the results to what you have and you will notice they are the same


    Solution using Eigen decomposition:

    Evals, Evec = np.linalg.eig(A.T @ A)
    P = Evec @ np.diag(np.sqrt(Evals)) @ Evec.T
    U = A @ Evec @ np.diag(1/np.sqrt(Evals)) @ Evec.T
    
    P
    array([[20., 10., 15.],
           [10., 25., 10.],
           [15., 10., 20.]])
    U
    array([[ 6.00000000e-01, -8.00000000e-01,  3.33066907e-16],
           [ 8.00000000e-01,  6.00000000e-01, -2.63677968e-16],
           [-9.71445147e-16, -1.66533454e-16,  1.00000000e+00]])