Scipy

Scipy provides many more computational algorithms and is built on Numpy.

%config InlineBackend.figure_format = 'svg'
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as la

Use LU decomposition to solve \(Ax=b\)

Create a random \(5 \times 5\) matrix

n = 5
A = np.random.rand(n,n)
print(A)

[[0.35402452 0.49935002 0.29029255 0.28445016 0.17101075]
 [0.12359041 0.77317746 0.33987877 0.69495421 0.55888056]
 [0.05678832 0.32513424 0.94160814 0.17099188 0.15532328]
 [0.32628048 0.67175687 0.67524423 0.05514395 0.69653413]
 [0.78898    0.54252607 0.38044898 0.2413887  0.65688315]]

Compute its LU factorization with pivoting using scipy.linalg.lu_factor

lu,piv = la.lu_factor(A)

Another option is to use scipy.linalg.lu which returns different arguments.

Create a right hand side vector

b = np.random.rand(n)

Solve \(A x = b\)

x = la.lu_solve((lu,piv),b)
print(x)

[ 0.20579505  3.46572465 -0.36764783 -2.20790792 -1.75929343]

Check that \(x\) solves the problem by computing \(Ax-b\)

print(A@x-b)

[-2.22044605e-16  2.77555756e-16  0.00000000e+00  1.11022302e-16
 -1.94289029e-16]

If we do not want the LU decomposition, we can directly solve, using scipy.linalg.solve

y = la.solve(A,b)
print(y)

[ 0.20579505  3.46572465 -0.36764783 -2.20790792 -1.75929343]

Note: There are similar methods in Numpy, see numpy.linalg.solve.

Sparse matrix

Scipy provides methods that can work on sparse matrices, see scipy.sparse.linalg

Sparse matrix formats are provided by

• csc_matrix

• csr_matrix

• coo_matrix

csc_matrix stores the entries column-wise.

from scipy.sparse import csc_matrix, csr_matrix
A = csc_matrix([[3, 2, 0], 
                [1, -1, 0], 
                [0, 5, 1]], dtype=float)
print(A)
print(A.todense())

<Compressed Sparse Column sparse matrix of dtype 'float64'
	with 6 stored elements and shape (3, 3)>
  Coords	Values
  (0, 0)	3.0
  (1, 0)	1.0
  (0, 1)	2.0
  (1, 1)	-1.0
  (2, 1)	5.0
  (2, 2)	1.0
[[ 3.  2.  0.]
 [ 1. -1.  0.]
 [ 0.  5.  1.]]

csr_matrix stores them row-wise.

A = csr_matrix([[3, 2, 0], 
                [1, -1, 0], 
                [0, 5, 1]], dtype=float)
print(A)
print(A.todense())

<Compressed Sparse Row sparse matrix of dtype 'float64'
	with 6 stored elements and shape (3, 3)>
  Coords	Values
  (0, 0)	3.0
  (0, 1)	2.0
  (1, 0)	1.0
  (1, 1)	-1.0
  (2, 1)	5.0
  (2, 2)	1.0
[[ 3.  2.  0.]
 [ 1. -1.  0.]
 [ 0.  5.  1.]]

We can construct sparse matrix by giving indices and values of non-zero entries. Indices not specified are assumed to be zero and need not be stored.

row = np.array([0, 0, 1, 1, 2, 2])
col = np.array([0, 1, 0, 1, 1, 2])
data = np.array([3.0, 2.0, 1.0, -1.0, 5.0, 1.0])
A = csc_matrix((data, (row, col)), shape=(3, 3))
print(A)
print(A.todense())

<Compressed Sparse Column sparse matrix of dtype 'float64'
	with 6 stored elements and shape (3, 3)>
  Coords	Values
  (0, 0)	3.0
  (1, 0)	1.0
  (0, 1)	2.0
  (1, 1)	-1.0
  (2, 1)	5.0
  (2, 2)	1.0
[[ 3.  2.  0.]
 [ 1. -1.  0.]
 [ 0.  5.  1.]]

Curve fitting

Suppose we are given some data set

\[ (x_i,y_i), \qquad i=0,1,...,n-1 \]

and we suspect there is a linear relation ship between them

\[ y = a + b x \]

But possibly the data also has some noise. Let us generate such a data set.

# Sample points
n = 10
x = np.linspace(0,1,n)

# Exact data
a = 1.0
b = 2.0
ye = a + b*x

# Add some noise
y = ye + 0.1*(2*np.random.rand(n)-1)
plt.plot(x,y,'o',label='Noisy Data')
plt.plot(x,ye,label='True function')
plt.legend();

Define the function we want to fit. The first argument is the independent variable, and remaining are the parameters that we want to find.

def f(x,a,b):
    return a + b*x

We can use curve_fit function from scipy to do the fitting.

from scipy.optimize import curve_fit
popt, pcov = curve_fit(f, x, y)
print('Fitted parameters = ',*popt)
plt.plot(x,y,'o',label='Noisy Data')
plt.plot(x,f(x,*popt),label='Fitted function')
plt.plot(x,ye,label='True function')
plt.legend();

Fitted parameters =  0.9874456177923452 2.044456278310629

We can also use the polyfit function or fit function from numpy to do this.

c = np.polyfit(x, y, deg=1)
print('Polynomial coefficients = ',c)

Polynomial coefficients =  [2.04445628 0.98744562]

The polynomials coefficients are returned in \(c\) in the following order

\[ p(x) = c[0] * x^{deg} + c[1] * x^{deg-1} + \ldots + c[-1] \]

p = np.poly1d(c)
plt.plot(x,y,'o',label='Noisy Data')
plt.plot(x,p(x),label='Fitted function')
plt.plot(x,ye,label='True function')
plt.legend();