In [1]:
#imports 
import numpy as np
from numpy import linalg,random

In this tutorial, we will make use of numpy's linalg module.

In [2]:
# first, initialize new numpy array 
a = np.array([[1],[2],[3]])
a
Out[2]:
array([[1],
       [2],
       [3]])

We can compute various norms using the linalg.norm function:

In [3]:
linalg.norm(a)  # L2 / Euclidean norm
Out[3]:
3.7416573867739413
In [4]:
linalg.norm(a, ord=1)  # L1 norm
Out[4]:
6.0
In [5]:
linalg.norm(a, ord=np.inf)  # inf-norm / max norm
Out[5]:
3.0

Computing the Determinant:

In [6]:
# initialize matrix A
A = np.array([[1,0,3],[4,-1,6],[7,8,3]])
A
Out[6]:
array([[ 1,  0,  3],
       [ 4, -1,  6],
       [ 7,  8,  3]])
In [7]:
# compute determinant
linalg.det(A)
Out[7]:
65.99999999999997
In [8]:
# initialize matrix B
B = np.array([[1,2,3],[2,4,6],[3,6,9]])
B
Out[8]:
array([[1, 2, 3],
       [2, 4, 6],
       [3, 6, 9]])
In [9]:
# compute determinant
linalg.det(B)
Out[9]:
0.0
In [10]:
# compute inverse
try:
    linalg.inv(B)  # throws exception because matrix is singular 
except Exception as e:
    print(e)
Singular matrix

Computing eigenvalues and eigenvectors:

In [11]:
eigvals, eigvecs = linalg.eig(A)  
In [12]:
print(eigvals) # vector of eigenvalues
[10.36660027 -1.         -6.36660027]
In [13]:
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.

In [14]:
# initialize symmetric matrix
S = np.matmul(A, A.transpose())
S
Out[14]:
array([[ 10,  22,  16],
       [ 22,  53,  38],
       [ 16,  38, 122]])
In [15]:
eigvals, eigvecs = linalg.eigh(S)   
In [16]:
print(eigvals)
[  0.73882306  41.21711541 143.04406153]
In [17]:
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:

In [18]:
A
Out[18]:
array([[ 1,  0,  3],
       [ 4, -1,  6],
       [ 7,  8,  3]])
In [19]:
L, V = linalg.eig(A)  # eigenvalues and eigenvectors
In [20]:
A_recomposed = np.matmul(V, np.matmul(np.diag(L), linalg.inv(V)))  # recompose matrix 
A_recomposed
Out[20]:
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]])
In [21]:
linalg.norm(A-A_recomposed)
Out[21]:
6.954641202201828e-15

Let's do the same thing with SVD:

In [22]:
A = random.rand(15).reshape(3,5)
A
Out[22]:
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]])
In [23]:
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]]
In [24]:
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]]
In [25]:
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:

In [26]:
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]]
In [27]:
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.

In [28]:
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]]
In [29]:
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):

In [30]:
A = random.rand(25).reshape(5,5)
A = np.matmul(A, A.T)
A
Out[30]:
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]])
In [31]:
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.

In [32]:
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]]
In [33]:
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.