Author avatar

Vivek Kumar

Working with SciPy Linear Algebra Module

Vivek Kumar

  • May 9, 2019
  • 11 Min read
  • 7 Views
  • May 9, 2019
  • 11 Min read
  • 7 Views
Data
SciPy

Introduction

In the previous guide in this series, you were introduced to the SciPy module, got familiarized with its organization, and learned a few basic methods. In this guide, you will learn about how SciPy extends the functionality of NumPy to provide rich linear algebra methods for developers and mathematicians.

By the end of this guide you'll have a strong understanding of the following topics:

  1. The differences between NumPy and SciPy Linear Algebra Modules
  2. Linear Algebra basic routines
  3. Decompositions
  4. Matrix Functions
  5. Special Matrices

The Baseline

Throughout this guide, we will be using the following libraries:

1
2
3
import numpy as np
import scipy as sp
from scipy import linalg
python

The Differences Between the NumPy and the SciPy Linear Algebra Modules

The table below highlights the difference between both the modules:

SciPy's linalg moduleNumPy's linalg module
It has all the methods which are available in NumPy's linalg module, including some extra methods.All the methods here are also available in the SciPy's linalg module.
It is always compiled using the BLAS/LAPACK support.It is optional for this module to compile using the BLAS/LAPACK support.

Linear Algebra Basic Routines

All the inputs and outputs to a SciPy linalg module are in a two-dimensional matrix. Therefore, before we understand some basic routines of the linear algebra module, let us first clear the differentiation between np.matrix and matrix formed using np.ndarray.

The np.matrix is a truly convenient method to constuct a two-dimensional matrix. On top of that, it also supports attributes like I and T for inverse and transpose respectively.

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
# Construction of a 2D matrix using np.matrix
mat = np.matrix([[1, 2], [3, 4]])
mat

# Output:
# matrix([[1, 2],
#         [3, 4]])

# Matrix inverse using I attribute
mat.I

# Output:
# matrix([[-2. ,  1. ],
#         [ 1.5, -0.5]])

# Matrix transpose using T attribute
mat.T

# Output:
# matrix([[1, 3],
#         [2, 4]])

# Dot product
mat * mat

# Output:
# array([[ 7, 10],
#        [15, 22]])
python

Even being quite convenient, the use of np.matrix method is still discouraged because whatever we can achieve using this method, we can also do it using the np.array method with some dependencies as shown:

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
# Construction of a 2D matrix using np.array
mat = np.array([[1, 2], [3, 4]]).reshape(2, 2)
mat

# Output:
# matrix([[1, 2],
#         [3, 4]])

# Matrix inverse using inv method
np.linalg.inv(mat)

# linalg.inv(mat) # Same answer using SciPy linalg

# Output:
# array([[-2. ,  1. ],
#        [ 1.5, -0.5]])

# Matrix transpose using the T attribute
mat.T

# Output:
# array([[1, 3],
#        [2, 4]])

# Dot product
np.dot(mat, mat)

# Output:
# array([[ 7, 10],
#        [15, 22]])
python

Now that we have an understanding of the pros and cons of each method, we can advance to the basic routines of linear algebra:

  1. Finding the inverse
  2. Solving equations
  3. Finding determinants
  4. Computing norms

We assume that you have some knowledge about these terms and will be skipping the introduction to these basic terms. Therefore, let us take an example for each and implement them using the linalg module of the SciPy.

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# Construction of a 3x3 matrix using np.array
mat = np.array([5, 4, 1, 2, 8, 6, 3, 7, 1]).reshape(3, 3)
mat

# array([[5, 4, 1],
#        [2, 8, 6],
#        [3, 7, 1]])

# ====================================
# 1. Finding the inverse of the matrix
# ====================================

linalg.inv(mat)

# Output:
# array([[ 0.29310345, -0.02586207, -0.13793103],
#        [-0.13793103, -0.01724138,  0.24137931],
#        [ 0.0862069 ,  0.19827586, -0.27586207]])

# ====================================
# 2. Finding determinant of the matrix
# ====================================

linalg.det(mat)

# Output:
# -116.0

# ================================
# 3. Computing norm of the matrix
# ================================

linalg.norm(mat) # Frobenius norm
# 14.317821063276353

linalg.norm(mat, 1) # L1 norm
# 19.0

linalg.norm(mat, -1) # Min column sum
# 8.0

linalg.norm(mat, np.inf) # L inf norm
# 16.0


# ===================================
# 4. Solving system of linear equations
# ===================================

A = np.array([[1, 3, 5], [2, 5, 1], [2, 3, 8]])
b = np.array([[10], [8], [3]])

linalg.solve(A, b)

# Solution:
# array([[-9.28],
#        [ 5.16],
#        [ 0.76]])
python

Decompositions

Decompositions are the most common procedures to be found in every other mathematical problem. Finding eigenvalue and eigenvector is one of the main components in the decomposition methods. We assume you are fairly well-versed with the various decomposition techniques as well as with the understanding of eigenvalue and eigenvectors. However, just to refresh your memories, you can refer to the given resource:

Here, we are going to take you through the implementation of SciPy in the computation of eigenvalues and eigenvectors, along with various decomposition methods.

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
# ============================
# Constructing some input data
# ============================
mat = np.array([5, 10, 15, 20, 25, 30, 35, 40, 45]).reshape(3, 3)
mat

# ===================================
# Finding eigenvalue and eigenvectors
# ===================================
eig = linalg.eig(mat)

eig[0] # Eigenvalues
# array([ 8.05842198e+01+0.j, -5.58421985e+00+0.j,  1.03401344e-15+0.j])

eig[1] # Eigenvectors
# array([[-0.23197069, -0.78583024,  0.40824829],
#        [-0.52532209, -0.08675134, -0.81649658],
#        [-0.8186735 ,  0.61232756,  0.40824829]])


# =============================
# Singular Value Decompositions
# =============================
U, s, Vh = linalg.svd(mat)
Sig = linalg.diagsvd(s, mat.shape[0], mat.shape[1])

U
# array([[-0.21483724,  0.88723069,  0.40824829],
#        [-0.52058739,  0.24964395, -0.81649658],
#        [-0.82633754, -0.38794278,  0.40824829]])

Sig
# array([[8.42405168e+01, 0.00000000e+00, 0.00000000e+00],
#        [0.00000000e+00, 5.34184757e+00, 0.00000000e+00],
#        [0.00000000e+00, 0.00000000e+00, 1.96374433e-15]])

Vh
# array([[-0.47967118, -0.57236779, -0.66506441],
#        [-0.77669099, -0.07568647,  0.62531805],
#        [-0.40824829,  0.81649658, -0.40824829]])

U.dot(Sig.dot(Vh)) # Verification
# array([[ 5., 10., 15.],
#        [20., 25., 30.],
#        [35., 40., 45.]])


# ================
# LU Decomposition
# ================
p, l, u = linalg.lu(mat)

p
# array([[0., 1., 0.],
#        [0., 0., 1.],
#        [1., 0., 0.]])

l
# array([[1.        , 0.        , 0.        ],
#        [0.14285714, 1.        , 0.        ],
#        [0.57142857, 0.5       , 1.        ]])

u
# array([[ 3.50000000e+01,  4.00000000e+01,  4.50000000e+01],
#        [ 0.00000000e+00,  4.28571429e+00,  8.57142857e+00],
#        [ 0.00000000e+00,  0.00000000e+00, -3.55271368e-15]])

p @ l @ u # Verification (Dot product)
# array([[ 5., 10., 15.],
#        [20., 25., 30.],
#        [35., 40., 45.]])

# ================
# QR Decomposition
# ================

q, r = linalg.qr(mat)

q
# array([[-0.12309149,  0.90453403,  0.40824829],
#        [-0.49236596,  0.30151134, -0.81649658],
#        [-0.86164044, -0.30151134,  0.40824829]])

r
# array([[-4.06201920e+01, -4.80056815e+01, -5.53911709e+01],
#        [ 0.00000000e+00,  4.52267017e+00,  9.04534034e+00],
#        [ 0.00000000e+00,  0.00000000e+00, -8.88178420e-15]])

q @ r # Verification (dot product)
# array([[ 5., 10., 15.],
#        [20., 25., 30.],
#        [35., 40., 45.]])

# ======================
# Cholesky Decomposition
# ======================

L = linalg.cholesky(np.array([[1,-2j],[2j,5]]), lower=True)

L
# array([[1.+0.j, 0.+0.j],
#        [0.+2.j, 1.+0.j]])

L @ L.T.conj() # Verification (dot prod)
# array([[1.+0.j, 0.-2.j],
#        [0.+2.j, 5.+0.j]])
python

Matrix Functions

While dealing with the matrices, there's often a requirement to implement trignometric, exponential, or related functions. Such functions can be implemented using the following methods:

ClassMathematical FunctionsSciPy Methods
Exponentialexpexpm
Logarithmloglogm
Trignometricsin, cos, tansinm, cosm, tanm
Hyperbolic trigonometricsinh, cosh, tanhsinhm, coshm, tanhm

Special Matrices

Both SciPy and NumPy provide separate methods to create special matrices which are frequently used by the developers. A few of them are listed below:

To get the complete list of other usable special matrices, you can refer to the given link:

Conclusion

In this guide, you have learned about the difference between SciPy and NumPy linalg modules, along with the implementation of various SciPy linalg methods.

To learn more about SciPy, you can refer the following guides:

0