from numpy import diag
from numpy import zerosFactorization
Matrix Decompositions
Many complex matrix operations cannot be solved efficiently or with stability using the limited precision of computers. Matrix decompositions are methods that reduce a matrix into constituent parts that make it easier to calculate more complex matrix operations. Matrix decomposition methods, also called matrix factorization methods, are a foundation of linear algebra in computers, even for basic operations such as solving systems of linear equations, calculating the inverse, and calculating the determinant of a matrix. In this tutorial, you will discover matrix decompositions and how to calculate them in Python.
What is a Matrix Decomposition
A matrix decomposition is a way of reducing a matrix into its constituent parts. It is an approach that can simplify more complex matrix operations that can be performed on the decomposed matrix rather than on the original matrix itself. A common analogy for matrix decomposition is the factoring of numbers, such as the factoring of 10 into 2 × 5. For this reason, matrix decomposition is also called matrix factorization. Like factoring real values, there are many ways to decompose a matrix, hence there are a range of different matrix decomposition techniques. Two simple and widely used matrix decomposition methods are the LU matrix decomposition and the QR matrix decomposition. Next, we will take a closer look at each of these methods.
LU Decomposition
from numpy import array
from scipy.linalg import luA = array([
[1,2,3],
[4,5,6],
[7,8,9]])
print(A)[[1 2 3]
[4 5 6]
[7 8 9]]P, L, U = lu(A)print(P)[[0. 1. 0.]
[0. 0. 1.]
[1. 0. 0.]]print(L)[[1. 0. 0. ]
[0.14285714 1. 0. ]
[0.57142857 0.5 1. ]]print(U)[[7. 8. 9. ]
[0. 0.85714286 1.71428571]
[0. 0. 0. ]]B = P.dot(L).dot(U)
print(B)[[1. 2. 3.]
[4. 5. 6.]
[7. 8. 9.]]QR Decomposition
from numpy import array
from scipy.linalg import qrA = array([
[1,2],
[3,4],
[5,6]])
print(A)[[1 2]
[3 4]
[5 6]]Q, R = qr(A)print(Q)[[-0.16903085 0.89708523 0.40824829]
[-0.50709255 0.27602622 -0.81649658]
[-0.84515425 -0.34503278 0.40824829]]print(R)[[-5.91607978 -7.43735744]
[ 0. 0.82807867]
[ 0. 0. ]]B = Q.dot(R)
print(B)[[1. 2.]
[3. 4.]
[5. 6.]]Cholesky Decomposition
from numpy import array
from numpy.linalg import choleskyA = array([
[2,1,1],
[1,2,1],
[1,1,2]])
print(A)[[2 1 1]
[1 2 1]
[1 1 2]]L = cholesky(A)
print(L)[[1.41421356 0. 0. ]
[0.70710678 1.22474487 0. ]
[0.70710678 0.40824829 1.15470054]]B = L.dot(L.T)
print(B)[[2. 1. 1.]
[1. 2. 1.]
[1. 1. 2.]]Eigendecomposition
Perhaps the most used type of matrix decomposition is the eigendecomposition that decomposes a matrix into eigenvectors and eigenvalues. This decomposition also plays a role in methods used in machine learning, such as in the Principal Component Analysis method or PCA.
Eigendecomposition of a Matrix
Eigenvectors and Eigenvalues
Calculation of Eigendecomposition
from numpy import array
from numpy.linalg import eigA = array([
[1,2,3],
[4,5,6],
[7,8,9]])
print(A)[[1 2 3]
[4 5 6]
[7 8 9]]values, vectors = eig(A)print(values)[ 1.61168440e+01 -1.11684397e+00 -1.30367773e-15]print(vectors)[[-0.23197069 -0.78583024 0.40824829]
[-0.52532209 -0.08675134 -0.81649658]
[-0.8186735 0.61232756 0.40824829]]Confirm an Eigenvector and Eigenvalue
from numpy import array
from numpy.linalg import eigA = array([
[1,2,3],
[4,5,6],
[7,8,9]])
print(A)[[1 2 3]
[4 5 6]
[7 8 9]]values, vectors = eig(A)B = A.dot(vectors[:,0])
print(B)[ -3.73863537 -8.46653421 -13.19443305]C = vectors[:,0] * values[0]
print(C)[ -3.73863537 -8.46653421 -13.19443305]Reconstruct Matrix
from numpy import array, diag
from numpy.linalg import eig, invA = array([
[1,2,3],
[4,5,6],
[7,8,9]])
print(A)[[1 2 3]
[4 5 6]
[7 8 9]]values, vectors = eig(A)Q = vectorsR = inv(Q)L = diag(values)B = Q.dot(L).dot(R)
print(B)[[1. 2. 3.]
[4. 5. 6.]
[7. 8. 9.]]Singular Value Decomposition
Calculate Singular-Value Decomposition
from numpy import array
from numpy.linalg import svdA = array([
[1,2],
[3,4],
[5,6]])
print(A)[[1 2]
[3 4]
[5 6]]U, s, VT = svd(A)print(U)[[-0.2298477 0.88346102 0.40824829]
[-0.52474482 0.24078249 -0.81649658]
[-0.81964194 -0.40189603 0.40824829]]print(s)[9.52551809 0.51430058]print(VT)[[-0.61962948 -0.78489445]
[-0.78489445 0.61962948]]Reconstruct Matrix
A = array([
[1,2],
[3,4],
[5,6]])
print(A)[[1 2]
[3 4]
[5 6]]U, s, VT = svd(A)Sigma = zeros((A.shape[0], A.shape[1]))Sigma[:A.shape[1], :A.shape[1]] = diag(s)B = U.dot(Sigma.dot(VT))
print(B)[[1. 2.]
[3. 4.]
[5. 6.]]A = array([
[1,2,3],
[4,5,6],
[7,8,9]])
print(A)[[1 2 3]
[4 5 6]
[7 8 9]]U, s , VT = svd(A)Sigma = diag(s)B = U.dot(Sigma.dot(VT))
print(B)[[1. 2. 3.]
[4. 5. 6.]
[7. 8. 9.]]Pseudoinverse
from numpy import array
from numpy.linalg import pinvA = array([
[0.1, 0.2],
[0.3, 0.4],
[0.5, 0.6],
[0.7, 0.8]])
print(A)[[0.1 0.2]
[0.3 0.4]
[0.5 0.6]
[0.7 0.8]]B = pinv(A)
print(B)[[-1.00000000e+01 -5.00000000e+00 1.42385628e-14 5.00000000e+00]
[ 8.50000000e+00 4.50000000e+00 5.00000000e-01 -3.50000000e+00]]U, s, VT = svd(A)d = 1.0 / sD = zeros(A.shape)D[:A.shape[1], :A.shape[1]] = diag(d)B = VT.T.dot(D.T).dot(U.T)
print(B)[[-1.00000000e+01 -5.00000000e+00 1.42578328e-14 5.00000000e+00]
[ 8.50000000e+00 4.50000000e+00 5.00000000e-01 -3.50000000e+00]]Dimensionality Reduction
A = array([
[1,2,3,4,5,6,7,8,9,10],
[11,12,13,14,15,16,17,18,19,20],
[21,22,23,24,25,26,27,28,29,30]])
print(A)[[ 1 2 3 4 5 6 7 8 9 10]
[11 12 13 14 15 16 17 18 19 20]
[21 22 23 24 25 26 27 28 29 30]]U, s, VT = svd(A)Sigma = zeros((A.shape[0], A.shape[1]))Sigma[:A.shape[0], :A.shape[0]] = diag(s)n_elemnts = 2Sigma = Sigma[:, :n_elemnts]VT = VT[:n_elemnts, :]B = U.dot(Sigma.dot(VT))
print(B)[[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.]
[11. 12. 13. 14. 15. 16. 17. 18. 19. 20.]
[21. 22. 23. 24. 25. 26. 27. 28. 29. 30.]]T = U.dot(Sigma)
print(T)[[-18.52157747 6.47697214]
[-49.81310011 1.91182038]
[-81.10462276 -2.65333138]]T = A.dot(VT.T)
print(T)[[-18.52157747 6.47697214]
[-49.81310011 1.91182038]
[-81.10462276 -2.65333138]]