Home

Linear Algebra and Polynomials

Rajesh Majumder

March 23, 2023

Numpy and Scipy packages are used for dealing with most of the mathematical operations ilke matrix operations polynomial, numarical integration, fourier series and many more.

A. Linear Algebra

1. Creating matrices and checking dimenssion

import numpy as np

A= np.array([[1,2,3],
            [4,5,6]])

B= np.array([[1,2],
             [3,4],
             [5,6]])

C= np.array([[7,8,9],
             [10,11,12]])

print("Matrix A:")
print(A)
print("Dimension:")
print(A.shape)

print("\nMatrix B:")
print(B)
print("Dimension:")
print(B.shape)

print("\nMatrix C:")
print(C)
print("Dimension:")
print(C.shape)

Matrix A:
[[1 2 3]
 [4 5 6]]
Dimension:
(2, 3)

Matrix B:
[[1 2]
 [3 4]
 [5 6]]
Dimension:
(3, 2)

Matrix C:
[[ 7  8  9]
 [10 11 12]]
Dimension:
(2, 3)

2. Transpose of a Matrix

A.transpose()

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

3. Matrix Addition and Subtraction

D= A+C

print(D)

[[ 8 10 12]
 [14 16 18]]

A+B

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_9560\3267467007.py in <module>
----> 1 A+B

ValueError: operands could not be broadcast together with shapes (2,3) (3,2) 

D-A

array([[ 7,  8,  9],
       [10, 11, 12]])

4. Adding each element of Matrix

$$A= \left[\begin{array}{ccc} a_{11} & a_{12} & . & . & . & a_{1p}\\ a_{21} & a_{22} & . & . & . & a_{2p}\\ . & . & . & . & . & .\\ . & . & . & . & . & .\\ . & . & . & . & . & .\\ a_{n1} & a_{n2} & . & . & . & a_{np}\\ \end{array}\right]$$

then, $$\sum_{i=0}^{n-1}{\sum_{j=0}^{p-1}{a_{(i,j)}}}$$

print("Element wise Sum of A : ",sum(sum(A)) )

Element wise Sum of A :  21

5. Element wise Matrix Multiplication (Hadmad Product)

A*C

array([[ 7, 16, 27],
       [40, 55, 72]])

6. Matrix Multiplication (Dotproduct)

$(AB)_{ij}= \sum_{k=0}^{p-1}{A_{ik}B_{kj}}$

print("A*B =", A.dot(B))

# Working Alternative 1
print("\nOutput by Alternative 1")
print("A*B =", np.dot(A,B))

# Working Alternative 2
print("\nOutput by Alternative 2")
print("A*B =", np.mat(A)*np.mat(B))

# Working Alternative 2
print("\nOutput by Alternative 3")
print("A*B =",  np.matmul(A,B))

A*B = [[22 28]
 [49 64]]

Output by Alternative 1
A*B = [[22 28]
 [49 64]]

Output by Alternative 2
A*B = [[22 28]
 [49 64]]

Output by Alternative 3
A*B = [[22 28]
 [49 64]]

• Note: np.mat() interpret the input as a matrix.

7. Cross Product

np.cross(A,C)

array([[-6, 12, -6],
       [-6, 12, -6]])

For calculating $X^TY$:

X= np.array([[2,1],
             [3,2],
             [-2,2]])

Y= np.array([[1,1],
             [4,2],
             [-2,1]])

np.dot(np.transpose(X),Y)

array([[18,  6],
       [ 5,  7]])

8. Diagonal Matrix

a= [1,2,3,4]

Diag1= np.diag(a)
print(Diag1)

np.diag()

[[1 0 0 0]
 [0 2 0 0]
 [0 0 3 0]
 [0 0 0 4]]

x = np.arange(16).reshape((4,4))

print(x)

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]

np.diag(x)

array([ 0,  5, 10, 15])

9. Upper and Lower Triangular Matrix

np.tril() for lower triangle

np.triu() for upper triangle

## Creating lower triangular matrix including the diagonal
np.tril([1,2,3,4],0)

array([[1, 0, 0, 0],
       [1, 2, 0, 0],
       [1, 2, 3, 0],
       [1, 2, 3, 4]])

x

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

np.tril(x,k=-1) ## Lower triangle of x

array([[ 0,  0,  0,  0],
       [ 4,  0,  0,  0],
       [ 8,  9,  0,  0],
       [12, 13, 14,  0]])

np.tril(x,k=1)

array([[ 0,  1,  0,  0],
       [ 4,  5,  6,  0],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

np.triu(x,k=0) #upper triangular matrix including the diagonal

array([[ 0,  1,  2,  3],
       [ 0,  5,  6,  7],
       [ 0,  0, 10, 11],
       [ 0,  0,  0, 15]])

np.triu(x,k=1) ## Upper triangle of x

array([[ 0,  1,  2,  3],
       [ 0,  0,  6,  7],
       [ 0,  0,  0, 11],
       [ 0,  0,  0,  0]])

np.triu(x,k=-1)

array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 0,  9, 10, 11],
       [ 0,  0, 14, 15]])

10. Determinant of a Matrix

To doing more complecated linear algebra operations NumPy has a submodule called linalg.

We can import linalg globally or directly we can use it.

import numpy as np
import numpy.linalg as la

xdet= la.det(x)
print(xdet)

-2.9582283945787796e-30

M= np.array([[-1,3,0],
             [2,1,-5],
             [1,4,-2]])

print("Determinant of M: ", np.linalg.det(M))

Determinant of M:  -21.00000000000001

11. Inverse of Matrix

print(la.inv(M))

[[-0.85714286 -0.28571429  0.71428571]
 [ 0.04761905 -0.0952381   0.23809524]
 [-0.33333333 -0.33333333  0.33333333]]

12. Trace of a Matrix

A= np.array(np.random.randint(low=10,high=50,size=16).reshape(4,4))
print("A:")
print(A)

print("\nTrace Of A is:")
print(np.trace(A))

A:
[[10 11 19 40]
 [18 19 28 37]
 [39 27 20 30]
 [17 39 23 33]]

Trace Of A is:
82

13. Eigen values, Eigen vectors

la.eig(x)

(array([ 3.24642492e+01, -2.46424920e+00,  4.52933966e-17, -1.52275132e-15]),
 array([[-0.11417645, -0.7327781 ,  0.01936986, -0.48846795],
        [-0.3300046 , -0.28974835,  0.38216644,  0.83598112],
        [-0.54583275,  0.15328139, -0.82244247, -0.20655839],
        [-0.76166089,  0.59631113,  0.42090617, -0.14095478]]))

la.eigvals(x)

array([ 3.24642492e+01, -2.46424920e+00,  4.52933966e-17, -1.52275132e-15])

14. Rank of a Matrix

la.matrix_rank(x)

2

15. Row Sum and Column Sum of Matrix

print(A)

## Row sum of A:
print("\nRow Sum of A:")
print(np.sum(A,1))

## Coluumn sum of A:
print("\nColumn Sum of A:")
print(np.sum(A,0))

[[1 2]
 [3 4]]

Row Sum of A:
[3 7]

Column Sum of A:
[4 6]

16. Add New Column or Row

np.append() is used for adding a new row.

np.insert() is used for adding a new column.

m= np.array([["Mon",18,20,22,17],
          ["Tue",11,13,15,12],
          ["Wed",15,21,20,29],
          ["Thu",11,10,22,21],
          ["Fir",12,14,20,18]])

print("Original Matrix:")
print(m)

## Add a New Row
newrow= ["Sat",12,22,20,28] # New row

print("\nNew Row Added")
print(np.append(m,[newrow],0)) # Adding the new row

## Add a New Column
newcol= [13,16,22,15,27] # New column

print("\nNew column Added but not in a proper way")
print(np.insert(m,[5],[newcol],1)) # Adding the new column

print("\nNew column Added with the proper way")
print(np.insert(m,[5],[[13],[16],[22],[15],[27]],1)) # Adding the new column

Original Matrix:
[['Mon' '18' '20' '22' '17']
 ['Tue' '11' '13' '15' '12']
 ['Wed' '15' '21' '20' '29']
 ['Thu' '11' '10' '22' '21']
 ['Fir' '12' '14' '20' '18']]

New Row Added
[['Mon' '18' '20' '22' '17']
 ['Tue' '11' '13' '15' '12']
 ['Wed' '15' '21' '20' '29']
 ['Thu' '11' '10' '22' '21']
 ['Fir' '12' '14' '20' '18']
 ['Sat' '12' '22' '20' '28']]

New column Added but not in a proper way
[['Mon' '18' '20' '22' '17' '13' '16' '22' '15' '27']
 ['Tue' '11' '13' '15' '12' '13' '16' '22' '15' '27']
 ['Wed' '15' '21' '20' '29' '13' '16' '22' '15' '27']
 ['Thu' '11' '10' '22' '21' '13' '16' '22' '15' '27']
 ['Fir' '12' '14' '20' '18' '13' '16' '22' '15' '27']]

New column Added with the proper way
[['Mon' '18' '20' '22' '17' '13']
 ['Tue' '11' '13' '15' '12' '16']
 ['Wed' '15' '21' '20' '29' '22']
 ['Thu' '11' '10' '22' '21' '15']
 ['Fir' '12' '14' '20' '18' '27']]

17. Linear Equations

Given the following equations:

$a_{11}x_1+a_{12}x_2+a_{13}x_3+....+a_{1m}x_m=b_1$ $a_{21}x_1+a_{22}x_2+a_{23}x_3+....+a_{2m}x_m=b_2$

$\:\:\:\:\;\;\;\;\;\;\;\;\;\;\;\:\;\:....$

$a_{n1}x_1+a_{n2}x_2+a_{n3}x_3+....+a_{nm}x_m=b_n$

In matrix notation:

$$Ax=b$$

Solution should be:

$$x= A^{-1}b$$

• Note: The $A$ matrix must be square and of fullrank, i.e. the inverse matrix needs to exists.

## Method 1:

A = np.array([[1, 2],
[3, 4]])
b = np.array([[5],
[6]])
Ainv = la.inv(A)
x = Ainv.dot(b)
print(x)

[[-4. ]
 [ 4.5]]

## Method 2:

x = np.linalg.solve(A, b)
print(x)

[[-4. ]
 [ 4.5]]

B. Polynomials

Form NumPy, poly1d sunmodule is supports to work with polynomials.

polu1d --> Its a one-dimensional polynomial class.

• Note: This forms part of the old polynomial API. Since version 1.4, the new polynomial API defined in numpy.polynomial is preferred.

1. Creating Polynomials

Creating 2nd order polynomial: $y= 4x^2 + 2x + 6$

## Creating 2nd order polynomial y= 4x^2 + 2x + 6
p1= np.poly1d([4,2,6])
print(p1)

   2
4 x + 2 x + 6

Creating 3rd order polynomial: $y= 7x^3 + 4x^2 + 5x + 6$

## Creating 3rd order polynomial y= 7x^3 + 4x^2 + 5x + 6 
p2=np.poly1d([7,4,5,6])
print(p2)

   3     2
7 x + 4 x + 5 x + 6

2. Solution of the Polynomials

p1(2) ## solution of the 2nd order polynomial when x=2 

26

p2(-1) ## solution of the 3rd order polynomial when x=-1

-2

3. Addition and Subtraction Polynomials

print(p1+p2) ## Addition

   3     2
7 x + 8 x + 7 x + 12

print(p2-p1) ## Subtraction

   3
7 x + 3 x

4. Multiplication and Division of Polynomials

print(p1*p2) ## Multiplication

    5      4      3      2
28 x + 30 x + 70 x + 58 x + 42 x + 36

print(p2/p1) ## Division

(poly1d([1.75 , 0.125]), poly1d([-5.75,  5.25]))

5. Coefficients & Order of the polynimial

## get the Coefficients & order of the polynimial

print("Coefficients of 2nd order polynomial:",p1.c)
print("And the order is:",p1.order)

print("\nCoefficients of 2nd order polynomial:",p2.c)
print("And the order is:",p2.order)

Coefficients of 2nd order polynomial: [4 2 6]
And the order is: 2

Coefficients of 2nd order polynomial: [7 4 5 6]
And the order is: 3

6. Derivative of Polynomial

print("1st Order Derivative of 2nd order Polunomial:\n",p1.deriv(m=1))

print("\n1st Order Derivative of 3nd order Polunomial:\n",p2.deriv(m=1))

1st Order Derivative of 2nd order Polunomial:
  
8 x + 2

1st Order Derivative of 3nd order Polunomial:
     2
21 x + 8 x + 5

7. Integration of Polynomial

print(p1.integ(m=1))
print(p1.integ(m=2))

print(p2.integ(m=1))
print(p2.integ(m=2))

       3     2
1.333 x + 1 x + 6 x
        4          3     2
0.3333 x + 0.3333 x + 3 x
      4         3       2
1.75 x + 1.333 x + 2.5 x + 6 x
      5          4          3     2
0.35 x + 0.3333 x + 0.8333 x + 3 x

8. To plot Polynomials

import numpy as np               # for creating range
import matplotlib.pyplot as plt  # for ploting

x= np.linspace(start=-10,stop=10,num=100)

## 2nd order Polynolinal
poly2= np.poly1d([4,5,6])

## Output
y= poly2(x)

## Plot of y
plt.plot(x,y,lw=3)    # Plot y versus x as lines

plt.show()

Example: Curve fitting by Polynomial

To fit the followibg data by a polynomial...

Var	Value
x	0, 10, 20, 30, 40, 50, 60, 70, 80, 90
y	76, 92, 106, 123, 132, 151, 179, 203, 227, 249

##import packages...
import numpy as np
import numpy.polynomial.polynomial as poly
import matplotlib.pyplot as plt

x= np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90])
y= np.array([76, 92, 106, 123, 132, 151, 179, 203, 227, 249])

## coefficients...
coeffs= poly.polyfit(x,y,2)
coeffs

array([7.81909091e+01, 1.10204545e+00, 9.12878788e-03])

## fitted values...
yfit= poly.polyval(x,coeffs)
yfit

array([ 78.19090909,  90.12424242, 103.88333333, 119.46818182,
       136.87878788, 156.11515152, 177.17727273, 200.06515152,
       224.77878788, 251.31818182])

## Plot...
plt.plot(x,y,x,yfit)
plt.show()

plt.plot(x,y,x,yfit,'ko')
plt.title('Fitting by polyfit',size=20)
plt.show()

Example: Curve fitting by user defined function

x= np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90])
y= np.array([76, 92, 106, 123, 132, 151, 179, 203, 227, 249])

## Creating a function
def fun(a,b,c,x):
    return a*x**2 + b*x + c

## Optimize the parameters
from scipy.optimize import curve_fit
par, var=curve_fit(fun,x,y)
a,b,c =par

## Plot
plt.plot(x,y,x,fun(a,b,c,x))
plt.show()