#
Advanced Matrix Algebra and Applications - Python Module
¶

## 1. An overview of Python 101¶

## Fundamental Mathematical Operations¶

In [193]:

    """Importing the required modules"""

    import numpy as np #module in which most of the operations of our interest are defined
    import matplotlib.pyplot as plt #for plotting
    from math import * #module for exp,log,power, trig functions etc...
    from numpy.linalg import * #module with linear algebraic functions
    import scipy.linalg as sp
    import matplotlib.image as mpimg
    from sklearn.decomposition import *

In [18]:

    a = 2.5
    b = 3e01
    c = sqrt(2)+1.55555j

In [19]:

    d = a + 2*b - np.conj(c) #linear combination of complex numbers
    print(d) #plain print statement
    print("Result : {}".format(np.round(d,3))) #formatted print statement

    (61.08578643762691+1.55555j)
    Result : (61.086+1.556j)

$\ $

# Lists¶

A list is one of the prominent data types in Python. It's a finite sequence of python objects (number/strings/matrices/other lists etc...) which is very used to define matrices.

In [59]:

    a = [1,2,3]
    b = [1,'vector',3]
    c = [[1,0],[0,1]]
    print(a+b+c) #prints the concatenated list
    print(type(a+b+c)) #outputs the data type
    print("Number of elements in a : {}".format(len(a)))

    [1, 2, 3, 1, 'vector', 3, [1, 0], [0, 1]]
    <class 'list'>
    Number of elements in a : 3

In [62]:

    #a[i] returns the i+1^th element of the list a
    print(a[1])
    print(b[1][0])
    print(c[1][0])

    2
    v
    0

In [25]:

    a + 1 #cannot perform element-wise addition in lists directly

    ---------------------------------------------------------------------------
    TypeError                                 Traceback (most recent call last)
    <ipython-input-25-ca42ed42e993> in <module>()
    ----> 1 a + 1

    TypeError: can only concatenate list (not "int") to list

In [27]:

    np.add(a,1) #This serves the purpose. However, it returns a different data type namely a numpy array which will be described later

Out[27]:

    array([2, 3, 4])

In [60]:

    np.add(a,list(range(3))) #range()

Out[60]:

    array([1, 3, 5])

In [61]:

    np.multiply(a,a[1]) #element-wise multiplication

Out[61]:

    array([2, 4, 6])

In [58]:

    np.multiply(a,[2,3,1/4])

Out[58]:

    array([2.  , 6.  , 0.75])

In [63]:

    """List methods"""

    a = list(range(3,10))
    print(np.sum(a))
    print(np.mean(a))
    print(sorted(a,reverse=True))
    print(a.index(5)) #returns the index of the specified element in the list
    a.append(100)  #appends the list
    print(a)
    print(a.pop()) #removes the last element
    print(a)
    print(a.pop(a[2])) #removes the element of the specified index
    a

    42
    6.0
    [9, 8, 7, 6, 5, 4, 3]
    2
    [3, 4, 5, 6, 7, 8, 9, 100]
    100
    [3, 4, 5, 6, 7, 8, 9]
    8

Out[63]:

    [3, 4, 5, 6, 7, 9]

## Arrays(Matrices)¶

Arrays are the most useful data type for one who works with matrices. A matrix is defined as a 2-dimensional array in python.

In [64]:

    A = np.array([[1,0],[0,1],[0,0]])
    print(A)

    print('\n',A.shape) #dimension of the array

    [[1 0]
     [0 1]
     [0 0]]

     (3, 2)

In [65]:

    A.reshape(2,3) #reshapes the array provided the new dim is compatible; does not alter the original array

Out[65]:

    array([[1, 0, 0],
           [1, 0, 0]])

In [66]:

    A.reshape(2,1,3)

Out[66]:

    array([[[1, 0, 0]],

           [[1, 0, 0]]])

In [67]:

    A.tolist()

Out[67]:

    [[1, 0], [0, 1], [0, 0]]

In [68]:

    B = 2*A #scalar multiplication
    B

Out[68]:

    array([[2, 0],
           [0, 2],
           [0, 0]])

In [69]:

    C = A + B #matrix addition
    print(C)

    [[3 0]
     [0 3]
     [0 0]]

In [70]:

    C[1][0] = 2
    C**2 #element-wise operation

Out[70]:

    array([[9, 0],
           [4, 9],
           [0, 0]], dtype=int32)

In [71]:

    C@C.T #Matrix multiplication

Out[71]:

    array([[ 9,  6,  0],
           [ 6, 13,  0],
           [ 0,  0,  0]])

In [72]:

    D = np.array([1,2]).reshape(2,1)
    C@D

Out[72]:

    array([[3],
           [8],
           [0]])

In [73]:

    print('Outer product : {}'.format(D@D.T))
    print('Inner product : {}'.format(D.T@D))

    Outer product : [[1 2]
     [2 4]]
    Inner product : [[5]]

In [76]:

    """Numpy methods for linear algebraic computation"""

    print('| I_3 | = {}'.format(det(np.identity(3))))

    | I_3 | = 1.0

In [77]:

    print('rank(A) = {}'.format(matrix_rank(A)))

    rank(A) = 2

In [78]:

    print('eigenvalues of CC^T: {}'.format(eig(C@C.T)[0]))

    eigenvalues of CC^T: [ 4.67544468 17.32455532  0.        ]

In [79]:

    P = eig(C@C.T)[1]
    print('eigenvectors of CC^T corresponding to lambda = {}: {}'.format(P[0][0],P[:,0]))
    print('eigenvectors of CC^T corresponding to lambda = {}: {}'.format(P[0][1],P[:,1]))
    print('eigenvectors of CC^T corresponding to lambda = {}: {}'.format(P[0][2],P[:,2]))

    eigenvectors of CC^T corresponding to lambda = -0.8112421851755608: [-0.81124219  0.58471028  0.        ]
    eigenvectors of CC^T corresponding to lambda = -0.584710284663765: [-0.58471028 -0.81124219  0.        ]
    eigenvectors of CC^T corresponding to lambda = 0.0: [0. 0. 1.]

## Loops, Logical Operators and Conditionals¶

In [107]:

    """Logical operations"""

    a,b = 2,3
    print('\n',a > b and a == b,'\n', a > b or a != b ,'\n',not(a==b))

     False
     True
     True

In [110]:

    """if - else statements"""

    if a < b:
        print ('a < b')

    elif a == b:
        print('a = b')

    else:
        print('a > b')

    #Such a conditional executes the 'first' set of statements for which the condition returns 'True'.
    #Note that the indentation is necessary for the compiler to understand the hierarchy.

    if a <= b:
        print ('a <= b')

    elif a < b:
        print('a < b')

    else:
        print('a > b') #The first two conditions return 'True', but only the first statement is executed.

    a < b
    a <= b

In [114]:

    """for loop""" #the statements are executed until the loop exhausts the list of elements in the iterable object

    for j in 'Python':
        print (j)

    print('\n')

    for [k,l] in [[1,2],[2,3]]:
        print (k+l)

    P
    y
    t
    h
    o
    n

    3
    5

In [115]:

    """while loop"""
    n = [1,2,3]
    while(len(n)):
        print(n.pop())

    3
    2
    1

## User defined functions¶

In [118]:

    def IP(a,b): #computes the euclidean inner product of two 1D arrays

        if len(a)==len(b):
            return np.sum([np.conj(a[k])*b[k] for k in range(len(a))])
        else:
            return("Error: The vectors should be of the same dimension")

    IP(np.array([0+1j,2,3]),np.array([1,1+2j,1]))

Out[118]:

    (5+3j)

## 2D Plots¶

In [95]:

    fig, axs = plt.subplots(1, 2, figsize=(16, 6))

    ax = axs[0]

    x = np.random.normal(size=5000)
    ax.hist(x,density=True,bins=50)

    ax = axs[1]

    n = np.arange(100)
    f = 1/n

    ax.plot(n,f,'r-',label='sequence')
    ax.set_xlabel('n')
    ax.set_ylabel('f(n)')
    ax.hlines(y=0,xmin=0,xmax=100,colors='g',label='limit')
    ax.set_title('$ f(n) = 1/n $')
    ax.legend()

    plt.show()

    C:\Users\Jaiprakash R\Anaconda3\lib\site-packages\ipykernel_launcher.py:12: RuntimeWarning: divide by zero encountered in true_divide
      if sys.path[0] == '':

## 2. Applications of Advanced Linear Algebra¶

## Polynomial Regression (Least Squares Problem)¶

In [89]:

    x = np.random.normal(5,2,size = 1000)
    y = np.array([0.5 + 3*x[i] - 0.5*x[i]**2 + np.random.normal(0,1,size=1) for i in range(len(x))]) #data with noise
    #true parameters b0= 0.5, b1 = 3 + b2 = -0.5, noise variance = 1

    plt.plot(x,y,'ro')

Out[89]:

    [<matplotlib.lines.Line2D at 0x1335d2b0>]

In [108]:

    x = x.reshape(-1,1)
    ones = np.ones(len(y)).reshape(-1,1)
    X = np.concatenate((ones,x,x**2),axis=1)

    B = []

    for i in range(1000):
        y = np.array([0.5 + 3*x[i] - 0.5*x[i]**2 + np.random.normal(0,1,size=1) for i in range(len(x))])

        b = inv((X.T@X))@X.T@y
        B.append(b)

In [109]:

    B = np.array(B)

In [110]:

    B.shape

Out[110]:

    (1000, 3, 1)

In [115]:

    fig, axs = plt.subplots(1, 3, figsize=(16, 6))

    axs[0].hist(B[:,0],density = True,bins=100)
    axs[0].set_title('$b_0$')

    axs[1].hist(B[:,1],density = True,bins=100)
    axs[1].set_title('$b_1$')

    axs[2].hist(B[:,2],density = True,bins=100)
    axs[2].set_title('$b_2$')

    #Note the large variance in data however

Out[115]:

    Text(0.5,1,'$b_2$')

## SVD to recover from Multicollinearity¶

In [184]:

    B = []
    U,s,Vh=sp.svd(X)

    #r=2
    Uhat = U[:,0:2]
    s_hat = np.diag(s[0:2])
    Vhat = Vh[0:2].T
    Vhat.shape

    for i in range(1000):

        y = np.array([0.5 + 3*x[i] - 0.5*x[i]**2 + np.random.normal(0,1,size=1) for i in range(len(x))])

        b = Vhat@inv(s_hat)@Uhat.T@y
        B.append(b)

In [186]:

    B = np.array(B)

    fig, axs = plt.subplots(1, 3, figsize=(16, 6))

    axs[0].hist(B[:,0],density = True,bins=100)
    axs[0].set_title('$b_0$')

    axs[1].hist(B[:,1],density = True,bins=100)
    axs[1].set_title('$b_1$')

    axs[2].hist(B[:,2],density = True,bins=100)
    axs[2].set_title('$b_2$')

    #Notice the reduction in variance. Although, the estimators are not centred at the true values anymore.

Out[186]:

    Text(0.5,1,'$b_2$')

## Image Compression - SVD¶

In [190]:

    fig = plt.figure(figsize=(16,6))
    img = mpimg.imread('12437-large.jpeg')
    print (img.shape)
    plt.axis('off')
    plt.imshow(img)

    (720, 1280, 3)

Out[190]:

    <matplotlib.image.AxesImage at 0x16de5128>

In [191]:

    img_r = np.reshape(img, (720, 1280*3))
    print (img_r.shape )

    (720, 3840)

In [194]:

    pca = PCA(n_components=256)
    ipca = pca.fit(img_r)
    img_c = pca.transform(img_r)
    print (img_c.shape)
    print (np.sum(ipca.explained_variance_ratio_))

    (720, 256)
    0.9995461532614928

In [199]:

    fig = plt.figure(figsize=(16,6))
    temp = pca.inverse_transform(img_c)
    print (temp.shape)
    temp = np.reshape(temp, (720,1280,3))
    print (temp.shape)
    plt.axis('off')
    plt.imshow(temp/255)
    plt.show()

    Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

    (720, 3840)
    (720, 1280, 3)

In [201]:

    pca = PCA(n_components=64)
    ipca = pca.fit(img_r)
    img_c = pca.transform(img_r)
    print (img_c.shape)
    print (np.sum(ipca.explained_variance_ratio_))

    (720, 64)
    0.9807895029617795

In [202]:

    fig = plt.figure(figsize=(16,6))
    temp = pca.inverse_transform(img_c)
    print (temp.shape)
    temp = np.reshape(temp, (720,1280,3))
    print (temp.shape)
    plt.axis('off')
    plt.imshow(temp/255)
    plt.show()

    Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

    (720, 3840)
    (720, 1280, 3)

In [203]:

    pca = PCA(n_components=16)
    ipca = pca.fit(img_r)
    img_c = pca.transform(img_r)
    print (img_c.shape)
    print (np.sum(ipca.explained_variance_ratio_))

    (720, 16)
    0.8891535421517809

In [204]:

    fig = plt.figure(figsize=(16,6))
    temp = pca.inverse_transform(img_c)
    print (temp.shape)
    temp = np.reshape(temp, (720,1280,3))
    print (temp.shape)
    plt.axis('off')
    plt.imshow(temp/255)
    plt.show()

    Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

    (720, 3840)
    (720, 1280, 3)