The linear algebra module contains a lot of matrix related functions, including linear equation solving, eigenvalue solvers, matrix functions (for example matrix-exponentiation), a number of different decomposition (SVD, LU, cholesky), etc. In this post we will look at how to use some of these functions:
1. Linear equation systems
Linear equation systems on the matrix form Ax = b where A is a matrix and x; b are vectors can be solved as shown in the program below:
from scipy.linalg import *
from numpy import array
A = array([[1,2,3], [4,5,6], [7,8,9]])
b = array([1,2,3])
x = solve(A, b)
print('\n')
print(x)
print('\n')
print(N.dot(A, x) - b)
The output of this program is shown below:
[-0.33333333 0.66666667 0. ]
[ 0.00000000e+00 -2.22044605e-16 0.00000000e+00]
------------------
(program exited with code: 0)
Press any key to continue . . .
Linear equation systems on the matrix form AX = B where A;B;X are matrices, can be solved as shown in the program below:
from scipy.linalg import *
import numpy as N
from numpy import array
A = N.random.rand(3,3)
B = N.random.rand(3,3)
X = solve(A, B)
print('\n')
print(X)
The output of this program is shown below:
[[ 0.07324201 2.28354346 0.78679028]
[ 0.310878 -0.99860178 -0.56368695]
[ 0.24402283 -0.72520934 0.32563823]]
------------------
(program exited with code: 0)
Press any key to continue . . .
2. Eigenvalues and eigenvectors
To calculate eigenvalues of a matrix, use the eigvals and for calculating both eigenvalues and eigenvectors, use the function eig as shown in the program below:
from scipy.linalg import *
import numpy as N
from numpy import array
A = N.random.rand(3,3)
evals = eigvals(A)
print('\n')
print(evals)
evals, evecs = eig(A)
print('\n')
print(evals)
print('\n')
print(evecs)
The output of this program is shown below:
[ 1.45980689+0.j -0.12519707+0.j 0.51421045+0.j]
[ 1.45980689+0.j -0.12519707+0.j 0.51421045+0.j]
[[ 0.52090604 0.49412854 -0.27850149]
[ 0.45795772 -0.77200849 -0.19138167]
[ 0.72036909 0.3997998 0.94117479]]
------------------
(program exited with code: 0)
Press any key to continue . . .
The eigenvectors corresponding to the nth eigenvalue (stored in evals[n]) is the nth column in evecs,
i.e., evecs[:,n]. We can verify this by multiplying eigenvectors with the matrix and compare to the product of the eigenvector and the eigenvalue:
n = 1
print('\n')
print(norm(N.dot(A, evecs[:,n]) - evals[n] * evecs[:,n]))
Add this code to the above program and the output will be as shown below:
5.453105385130469e-16
------------------
(program exited with code: 0)
Press any key to continue . . .
3. Matrix operations
The following program shows some of the matrix operations performed using scipy:
from scipy.linalg import *
import numpy as N
import math
from numpy import array
A = N.random.rand(3,3)
# the matrix inverse
inv(A)
print('\n')
print(A)
# determinant
det(A)
print('\n')
print(A)
print('\n')
# norms of various orders
print(norm(A, ord=2), norm(A, ord=math.inf))
The output of this program is shown below:
[[0.32326371 0.7398846 0.15224097]
[0.27898164 0.65034488 0.77975491]
[0.49160609 0.71130144 0.5823993 ]]
[[0.32326371 0.7398846 0.15224097]
[0.27898164 0.65034488 0.77975491]
[0.49160609 0.71130144 0.5823993 ]]
1.6385878124754003 1.7853068270024153
------------------
(program exited with code: 0)
Press any key to continue . . .
4. Sparse matrices
Scipy has a good support for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc). There are many possible strategies for storing sparse matrices in an efficient way. Some of the most common are the so-called coordinate form (COO), list of list (LIL) form, and compressed-sparse column CSC (and row, CSR). Each format has some advantages and disadvantages. Most computational algorithms (equation solving, matrix-matrix multiplication, etc) can be efficiently implemented using CSR or CSC formats, but they are not so intuitive and not so easy to initialize. So often a sparse matrix is initially created in COO or LIL format (where we can e efficiently add elements to the sparse matrix data), and then converted to CSC or CSR before used in real calculations.
When we create a sparse matrix we have to choose which format it should be stored in. See the following example:
from scipy.linalg import *
from scipy.sparse import *
import numpy as N
import math
from numpy import array
# dense matrix
M = array([[1,0,0,0], [0,3,0,0], [0,1,1,0], [1,0,0,1]])
print('\n')
print(M)
# convert from dense to sparse
A = csr_matrix(M)
print('\n')
print(A)
print('\n')
# convert from sparse to dense
print(A.todense())
The output of this program is shown below:
[[1 0 0 0]
[0 3 0 0]
[0 1 1 0]
[1 0 0 1]]
(0, 0) 1
(1, 1) 3
(2, 1) 1
(2, 2) 1
(3, 0) 1
(3, 3) 1
[[1 0 0 0]
[0 3 0 0]
[0 1 1 0]
[1 0 0 1]]
------------------
(program exited with code: 0)
Press any key to continue . . .
There is a more e efficient way to create sparse matrices: first create an empty matrix and populate it with using matrix indexing (avoids creating a potentially large dense matrix) . See the following program:
from scipy.linalg import *
from scipy.sparse import *
import numpy as N
import math
from numpy import array
A = lil_matrix((4,4)) # empty 4x4 sparse matrix
A[0,0] = 1
A[1,1] = 3
A[2,2] = A[2,1] = 1
A[3,3] = A[3,0] = 1
print('\n')
print(A)
print('\n')
# convert from sparse to dense
print(A.todense())
print('\n')
print(A)
print('\n')
print(csr_matrix(A))
print('\n')
print(csc_matrix(A))
print('\n')
The output of this program is shown below:
(0, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 0) 1.0
(3, 3) 1.0
[[1. 0. 0. 0.]
[0. 3. 0. 0.]
[0. 1. 1. 0.]
[1. 0. 0. 1.]]
(0, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 0) 1.0
(3, 3) 1.0
(0, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 0) 1.0
(3, 3) 1.0
(0, 0) 1.0
(3, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 3) 1.0
------------------
(program exited with code: 0)
Press any key to continue . . .
It is possible to compute with sparse matrices like we compute with dense matrices. See the following program:
from scipy.linalg import *
from scipy.sparse import *
import numpy as N
import math
from numpy import array
A = lil_matrix((4,4)) # empty 4x4 sparse matrix
print('\n')
print(A.todense())
print('\n')
print((A * A).todense())
print('\n')
print(A.todense())
print('\n')
print(A.dot(A).todense())
print('\n')
v = array([1,2,3,4])[:,N.newaxis]
print(v)
print('\n')
# sparse matrix - dense vector multiplication
print(A*v)
print('\n')
# same result with dense matrix - dense vector multiplcation
print(A.todense() * v)
The output of this program is shown below:
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[1]
[2]
[3]
[4]]
[[0.]
[0.]
[0.]
[0.]]
[[0.]
[0.]
[0.]
[0.]]
------------------
(program exited with code: 0)
Press any key to continue . . .
Here I am ending this post regarding SciPy linear algebra module. Next post will discuss about optimization with SciPy. Till we meet next keep learning and practicing Python as Python is easy to learn!
1. Linear equation systems
Linear equation systems on the matrix form Ax = b where A is a matrix and x; b are vectors can be solved as shown in the program below:
from scipy.linalg import *
from numpy import array
A = array([[1,2,3], [4,5,6], [7,8,9]])
b = array([1,2,3])
x = solve(A, b)
print('\n')
print(x)
print('\n')
print(N.dot(A, x) - b)
The output of this program is shown below:
[-0.33333333 0.66666667 0. ]
[ 0.00000000e+00 -2.22044605e-16 0.00000000e+00]
------------------
(program exited with code: 0)
Press any key to continue . . .
Linear equation systems on the matrix form AX = B where A;B;X are matrices, can be solved as shown in the program below:
from scipy.linalg import *
import numpy as N
from numpy import array
A = N.random.rand(3,3)
B = N.random.rand(3,3)
X = solve(A, B)
print('\n')
print(X)
The output of this program is shown below:
[[ 0.07324201 2.28354346 0.78679028]
[ 0.310878 -0.99860178 -0.56368695]
[ 0.24402283 -0.72520934 0.32563823]]
------------------
(program exited with code: 0)
Press any key to continue . . .
2. Eigenvalues and eigenvectors
To calculate eigenvalues of a matrix, use the eigvals and for calculating both eigenvalues and eigenvectors, use the function eig as shown in the program below:
from scipy.linalg import *
import numpy as N
from numpy import array
A = N.random.rand(3,3)
evals = eigvals(A)
print('\n')
print(evals)
evals, evecs = eig(A)
print('\n')
print(evals)
print('\n')
print(evecs)
The output of this program is shown below:
[ 1.45980689+0.j -0.12519707+0.j 0.51421045+0.j]
[ 1.45980689+0.j -0.12519707+0.j 0.51421045+0.j]
[[ 0.52090604 0.49412854 -0.27850149]
[ 0.45795772 -0.77200849 -0.19138167]
[ 0.72036909 0.3997998 0.94117479]]
------------------
(program exited with code: 0)
Press any key to continue . . .
The eigenvectors corresponding to the nth eigenvalue (stored in evals[n]) is the nth column in evecs,
i.e., evecs[:,n]. We can verify this by multiplying eigenvectors with the matrix and compare to the product of the eigenvector and the eigenvalue:
n = 1
print('\n')
print(norm(N.dot(A, evecs[:,n]) - evals[n] * evecs[:,n]))
Add this code to the above program and the output will be as shown below:
5.453105385130469e-16
------------------
(program exited with code: 0)
Press any key to continue . . .
3. Matrix operations
The following program shows some of the matrix operations performed using scipy:
from scipy.linalg import *
import numpy as N
import math
from numpy import array
A = N.random.rand(3,3)
# the matrix inverse
inv(A)
print('\n')
print(A)
# determinant
det(A)
print('\n')
print(A)
print('\n')
# norms of various orders
print(norm(A, ord=2), norm(A, ord=math.inf))
The output of this program is shown below:
[[0.32326371 0.7398846 0.15224097]
[0.27898164 0.65034488 0.77975491]
[0.49160609 0.71130144 0.5823993 ]]
[[0.32326371 0.7398846 0.15224097]
[0.27898164 0.65034488 0.77975491]
[0.49160609 0.71130144 0.5823993 ]]
1.6385878124754003 1.7853068270024153
------------------
(program exited with code: 0)
Press any key to continue . . .
4. Sparse matrices
Scipy has a good support for sparse matrices, with basic linear algebra operations (such as equation solving, eigenvalue calculations, etc). There are many possible strategies for storing sparse matrices in an efficient way. Some of the most common are the so-called coordinate form (COO), list of list (LIL) form, and compressed-sparse column CSC (and row, CSR). Each format has some advantages and disadvantages. Most computational algorithms (equation solving, matrix-matrix multiplication, etc) can be efficiently implemented using CSR or CSC formats, but they are not so intuitive and not so easy to initialize. So often a sparse matrix is initially created in COO or LIL format (where we can e efficiently add elements to the sparse matrix data), and then converted to CSC or CSR before used in real calculations.
When we create a sparse matrix we have to choose which format it should be stored in. See the following example:
from scipy.linalg import *
from scipy.sparse import *
import numpy as N
import math
from numpy import array
# dense matrix
M = array([[1,0,0,0], [0,3,0,0], [0,1,1,0], [1,0,0,1]])
print('\n')
print(M)
# convert from dense to sparse
A = csr_matrix(M)
print('\n')
print(A)
print('\n')
# convert from sparse to dense
print(A.todense())
The output of this program is shown below:
[[1 0 0 0]
[0 3 0 0]
[0 1 1 0]
[1 0 0 1]]
(0, 0) 1
(1, 1) 3
(2, 1) 1
(2, 2) 1
(3, 0) 1
(3, 3) 1
[[1 0 0 0]
[0 3 0 0]
[0 1 1 0]
[1 0 0 1]]
------------------
(program exited with code: 0)
Press any key to continue . . .
There is a more e efficient way to create sparse matrices: first create an empty matrix and populate it with using matrix indexing (avoids creating a potentially large dense matrix) . See the following program:
from scipy.linalg import *
from scipy.sparse import *
import numpy as N
import math
from numpy import array
A = lil_matrix((4,4)) # empty 4x4 sparse matrix
A[0,0] = 1
A[1,1] = 3
A[2,2] = A[2,1] = 1
A[3,3] = A[3,0] = 1
print('\n')
print(A)
print('\n')
# convert from sparse to dense
print(A.todense())
print('\n')
print(A)
print('\n')
print(csr_matrix(A))
print('\n')
print(csc_matrix(A))
print('\n')
The output of this program is shown below:
(0, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 0) 1.0
(3, 3) 1.0
[[1. 0. 0. 0.]
[0. 3. 0. 0.]
[0. 1. 1. 0.]
[1. 0. 0. 1.]]
(0, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 0) 1.0
(3, 3) 1.0
(0, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 0) 1.0
(3, 3) 1.0
(0, 0) 1.0
(3, 0) 1.0
(1, 1) 3.0
(2, 1) 1.0
(2, 2) 1.0
(3, 3) 1.0
------------------
(program exited with code: 0)
Press any key to continue . . .
It is possible to compute with sparse matrices like we compute with dense matrices. See the following program:
from scipy.linalg import *
from scipy.sparse import *
import numpy as N
import math
from numpy import array
A = lil_matrix((4,4)) # empty 4x4 sparse matrix
print('\n')
print(A.todense())
print('\n')
print((A * A).todense())
print('\n')
print(A.todense())
print('\n')
print(A.dot(A).todense())
print('\n')
v = array([1,2,3,4])[:,N.newaxis]
print(v)
print('\n')
# sparse matrix - dense vector multiplication
print(A*v)
print('\n')
# same result with dense matrix - dense vector multiplcation
print(A.todense() * v)
The output of this program is shown below:
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
[[1]
[2]
[3]
[4]]
[[0.]
[0.]
[0.]
[0.]]
[[0.]
[0.]
[0.]
[0.]]
------------------
(program exited with code: 0)
Press any key to continue . . .
Here I am ending this post regarding SciPy linear algebra module. Next post will discuss about optimization with SciPy. Till we meet next keep learning and practicing Python as Python is easy to learn!
0 comments:
Post a Comment