#imports
import numpy as np
from numpy import linalg,random
In this tutorial, we will make use of numpy's linalg
module.
# first, initialize new numpy array
a = np.array([[1],[2],[3]])
a
array([[1], [2], [3]])
We can compute various norms using the linalg.norm
function:
linalg.norm(a) # L2 / Euclidean norm
3.7416573867739413
linalg.norm(a, ord=1) # L1 norm
6.0
linalg.norm(a, ord=np.inf) # inf-norm / max norm
3.0
Computing the Determinant:
# initialize matrix A
A = np.array([[1,0,3],[4,-1,6],[7,8,3]])
A
array([[ 1, 0, 3], [ 4, -1, 6], [ 7, 8, 3]])
# compute determinant
linalg.det(A)
65.99999999999997
# initialize matrix B
B = np.array([[1,2,3],[2,4,6],[3,6,9]])
B
array([[1, 2, 3], [2, 4, 6], [3, 6, 9]])
# compute determinant
linalg.det(B)
0.0
# compute inverse
try:
linalg.inv(B) # throws exception because matrix is singular
except Exception as e:
print(e)
Singular matrix
Computing eigenvalues and eigenvectors:
eigvals, eigvecs = linalg.eig(A)
print(eigvals) # vector of eigenvalues
[10.36660027 -1. -6.36660027]
print(eigvecs) # eigenvectors are columns of the matrix
[[-0.2603943 -0.75856745 -0.30110589] [-0.5207886 0.4108907 -0.60221177] [-0.8130031 0.50571163 0.73937557]]
Note that eig
does not guarantee that eigenvalues are returned in sorted order.
For symmetric matrices, use eigh
instead. This is more efficient and numerically reliable because it takes advantage of the fact that the matrix is symmetric. Moreover, the eigenvalues are returned in ascending order.
# initialize symmetric matrix
S = np.matmul(A, A.transpose())
S
array([[ 10, 22, 16], [ 22, 53, 38], [ 16, 38, 122]])
eigvals, eigvecs = linalg.eigh(S)
print(eigvals)
[ 0.73882306 41.21711541 143.04406153]
print(eigvecs)
[[ 9.21673813e-01 3.45544439e-01 1.76398478e-01] [-3.87965696e-01 8.20934884e-01 4.18985124e-01] [-3.36839069e-05 -4.54604175e-01 8.90693574e-01]]
If you only want eigenvalues, you can use the linalg.eigvals
and linalg.eigvalsh
functions.
Let's try to recompose the matrix from eigenvalues and eigenvectors:
A
array([[ 1, 0, 3], [ 4, -1, 6], [ 7, 8, 3]])
L, V = linalg.eig(A) # eigenvalues and eigenvectors
A_recomposed = np.matmul(V, np.matmul(np.diag(L), linalg.inv(V))) # recompose matrix
A_recomposed
array([[ 1.00000000e+00, 2.88657986e-15, 3.00000000e+00], [ 4.00000000e+00, -1.00000000e+00, 6.00000000e+00], [ 7.00000000e+00, 8.00000000e+00, 3.00000000e+00]])
linalg.norm(A-A_recomposed)
6.954641202201828e-15
Let's do the same thing with SVD:
A = random.rand(15).reshape(3,5)
A
array([[0.35457622, 0.40306362, 0.16841218, 0.93304962, 0.45510363], [0.6547846 , 0.99911555, 0.15688209, 0.67689878, 0.24453189], [0.91685889, 0.4879655 , 0.3128614 , 0.76569169, 0.15619979]])
U, S, V_T = linalg.svd(A) #Notice S is a one dimensional vector
print(f"U =\n{U}\n{15*'-'}\nS =\n{S}\n{15*'-'}\nV_T = \n{V_T}")
U = [[-0.50644226 0.83263508 -0.2241318 ] [-0.62425083 -0.53335249 -0.57082924] [-0.5948337 -0.14917759 0.78988538]] --------------- S = [2.17310785 0.50113546 0.42141805] --------------- V_T = [[-0.52169559 -0.5145099 -0.1699526 -0.62148301 -0.21906223] [-0.3806819 -0.53891491 0.01971627 0.6019134 0.44940353] [ 0.64299912 -0.65309715 0.28433818 0.02204218 -0.28050344] [-0.41163103 0.12980763 0.63247988 0.27075682 -0.58341222] [-0.00519902 0.03826712 0.6998917 -0.42150209 0.57532269]]
A_recomposed = np.matmul(U, np.matmul(np.diag([*S, 0, 0])[0:3], V_T))
print(f"U * S* V_T = \n{A_recomposed}")
U * S* V_T = [[0.35457622 0.40306362 0.16841218 0.93304962 0.45510363] [0.6547846 0.99911555 0.15688209 0.67689878 0.24453189] [0.91685889 0.4879655 0.3128614 0.76569169 0.15619979]]
print(f"Norm of the diff is: {linalg.norm(A-A_recomposed)}")
Norm of the diff is: 1.3423487087819982e-15
Let's see what happens if the number of rows were more than columns:
A = A.T
print(A)
[[0.35457622 0.6547846 0.91685889] [0.40306362 0.99911555 0.4879655 ] [0.16841218 0.15688209 0.3128614 ] [0.93304962 0.67689878 0.76569169] [0.45510363 0.24453189 0.15619979]]
U, S, V_T = linalg.svd(A) #Notice S is a one dimensional vector
print(f"U =\n{U}\n{15*'-'}\nS =\n{S}\n{15*'-'}\nV_T = \n{V_T}")
U = [[-0.52169559 0.3806819 -0.64299912 -0.41163103 -0.00519902] [-0.5145099 0.53891491 0.65309715 0.12980763 0.03826712] [-0.1699526 -0.01971627 -0.28433818 0.63247988 0.6998917 ] [-0.62148301 -0.6019134 -0.02204218 0.27075682 -0.42150209] [-0.21906223 -0.44940353 0.28050344 -0.58341222 0.57532269]] --------------- S = [2.17310785 0.50113546 0.42141805] --------------- V_T = [[-0.50644226 -0.62425083 -0.5948337 ] [-0.83263508 0.53335249 0.14917759] [ 0.2241318 0.57082924 -0.78988538]]
Notice the relation with the previous example U and V changed places.
A_recomposed = np.matmul(U, np.matmul(np.diag([*S, 0, 0])[:,0:3], V_T))
print(f"U * S* V_T = \n{A_recomposed}")
U * S* V_T = [[0.35457622 0.6547846 0.91685889] [0.40306362 0.99911555 0.4879655 ] [0.16841218 0.15688209 0.3128614 ] [0.93304962 0.67689878 0.76569169] [0.45510363 0.24453189 0.15619979]]
print(f"Norm of the diff is: {linalg.norm(A-A_recomposed)}")
Norm of the diff is: 1.64976489155517e-15
Finally if a symmetric matrix has distinct eigenvalues, Then SVD coincides with eigendecomposition(Up to sign change):
A = random.rand(25).reshape(5,5)
A = np.matmul(A, A.T)
A
array([[2.53244456, 1.61325776, 1.29089283, 1.74799725, 1.16011574], [1.61325776, 1.76330046, 0.77240397, 1.47852544, 0.84832181], [1.29089283, 0.77240397, 0.81080074, 0.95987403, 0.92398355], [1.74799725, 1.47852544, 0.95987403, 1.47211859, 1.05536618], [1.16011574, 0.84832181, 0.92398355, 1.05536618, 1.53650662]])
U, S, V_T = linalg.svd(A)
eigvals, eigvecs = linalg.eigh(A)
print(list(S), list(reversed(eigvals)), sep='\n')
[6.611344461630094, 0.9090431034896154, 0.5411320677863699, 0.04292794173897154, 0.010723392314945436] [6.611344461630095, 0.909043103489615, 0.5411320677863696, 0.042927941738971416, 0.01072339231494516]
Here, numbers are not exactly equal. That is because floating point arithmatic is not exact.
print(eigvecs[:,-1::-1])
[[-0.58434444 -0.20128852 -0.65267107 -0.4148717 -0.1411611 ] [-0.45171902 -0.46700378 0.62606541 -0.24322496 0.35600719] [-0.32515894 0.25564588 -0.25885168 0.49702312 0.71755186] [-0.46436722 -0.11149579 0.15748727 0.65306448 -0.56624758] [-0.36486471 0.8146191 0.3004266 -0.30854157 -0.13347441]]
print(f"U =\n{U}\n{15*'-'}\n V_T =\n{V_T}")
U = [[-0.58434444 0.20128852 -0.65267107 0.4148717 0.1411611 ] [-0.45171902 0.46700378 0.62606541 0.24322496 -0.35600719] [-0.32515894 -0.25564588 -0.25885168 -0.49702312 -0.71755186] [-0.46436722 0.11149579 0.15748727 -0.65306448 0.56624758] [-0.36486471 -0.8146191 0.3004266 0.30854157 0.13347441]] --------------- V_T = [[-0.58434444 -0.45171902 -0.32515894 -0.46436722 -0.36486471] [ 0.20128852 0.46700378 -0.25564588 0.11149579 -0.8146191 ] [-0.65267107 0.62606541 -0.25885168 0.15748727 0.3004266 ] [ 0.4148717 0.24322496 -0.49702312 -0.65306448 0.30854157] [ 0.1411611 -0.35600719 -0.71755186 0.56624758 0.13347441]]
U and eigenvectors are same(up to sign). Also, U and V are equal.
When $A$ is symmetric we have $A = A^T$. By SVD we have $A = U \times diag(S) \times V^T$ so $A^T = V \times diag(S) \times U^T$. For reasons out of scope of this class SVD is unique (up to change of sign) so:
$A = A^T \implies V = U$
in this example numbers are not exactly the same but $V ≈ E$. Again becuase of floating point arithmatic.
By defenition of SVD, V and U are orthogonal and we just saw $U = V$:
$V^T \times V = I \implies V^T = V^{-1} \implies V^T = U^{-1}$
Since $A = U \times diag(S) \times V^T$ We have:
$A = U \times diag(S) \times U^{-1}$
If eigenvalues are different this is the only eigendecomposition.