Neural Networks



In this last chapter, we will take a look at the mathematical tools used to model and understand biological neural networks. We're going to see how to use Python to do operations on matrices and then we will take a look at some of the fundamental concepts in Graph Theory which will be the starting point for network modeling.
What we're going to do in this Notebook, is that we're going to see how to use SymPy and NumPy to implement some matrix properties and operations, and after that we're going to use another Python Package "NetworkX" which is used to study and implement graphs and networks.

Table of contents

The necessary libraries and functions are :

In [5]:
import numpy as np
import sympy as sp
import networkx as nx
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from itertools import combinations
from collections import Counter

Matrices

Matrix addition and substraction

  • If we have the following matrices :
$$ A = \begin{bmatrix} 2 & 5 & 7 \\ 4 & 8 & 3 \\ 9 & 1 & 4 \\\end{bmatrix} \ et \ B = \begin{bmatrix} 1 & 5 & 3 \\ 3 & 6 & 1 \\ 7 & 5 & 9 \\\end{bmatrix}$$

And we want to find $A+B$ or $A-B$

  • We can proceed with NumPy with the following manner :
In [6]:
# Creating the arrays 
A = np.array([2,5,7,4,8,3,9,1,4]) 
B = np.array([1,5,3,3,6,1,7,5,9])
# Reshaping them to a 3x3 array
for mat in (A,B) : 
    mat.shape = (3,3) #

print("A + B \n" , A + B, "\nA - B \n", A - B) # Matrix addition and substraction

A + B 
 [[ 3 10 10]
 [ 7 14  4]
 [16  6 13]] 
A - B 
 [[ 1  0  4]
 [ 1  2  2]
 [ 2 -4 -5]]
  • Or we can choose SymPy :
In [ ]:
A = sp.Matrix(3, 3,[2, 5, 7, 4, 8, 3, 9, 1, 4])  # First argument is number of row and second is number of columns
B = sp.Matrix(3, 3, [1, 5, 3, 3, 6, 1, 7, 5, 9])

And the addition/substration are as simple as :

In [14]:
A + B
Out[14]:
$\displaystyle \left[\begin{matrix}3 & 10 & 10\\7 & 14 & 4\\16 & 6 & 13\end{matrix}\right]$
In [15]:
A - B
Out[15]:
$\displaystyle \left[\begin{matrix}1 & 0 & 4\\1 & 2 & 2\\2 & -4 & -5\end{matrix}\right]$

Matrix multiplication or dot product

  • If we take the following vectors and we want to calculate their dot product :
$$ A = \begin{bmatrix} 2 & 3 & 4 \end{bmatrix} \ et \ B = \begin{bmatrix} 4 \\ 1 \\ 5 \end{bmatrix}$$

With NumPy, this can be done using the function np.dot() :

In [11]:
A = np.array([[2,3,4]])
B = np.array([4,1,5])
B.shape = (3,1)
print('A.B \n', np.dot(A,B))

A.B 
 [[31]]

And with SymPy, the multiplication symbol will do the trick :

In [16]:
A = sp.Matrix(1,3, [2,3,4]) # 1x3
B = sp.Matrix(3,1, [4,1,5]) # 3x1
A * B
Out[16]:
$\displaystyle \left[\begin{matrix}31\end{matrix}\right]$
  • Let's take another example
$$ A = \begin{bmatrix} 2 & 8 \\ 5 & -6 \end{bmatrix} \ et \ B = \begin{bmatrix} 7 & -4 \\ 3 & 12 \end{bmatrix}$$
  • This time let's calculate $A.B$ and $B.A$

In SymPy :

In [104]:
A = sp.Matrix(2,2,[2,8,5,-6])
B = sp.Matrix(2,2,[7,-4,3,12])
A*B
Out[104]:
$\displaystyle \left[\begin{matrix}38 & 88\\17 & -92\end{matrix}\right]$
In [105]:
B*A
Out[105]:
$\displaystyle \left[\begin{matrix}-6 & 80\\66 & -48\end{matrix}\right]$

And with NumPy :

In [110]:
A = np.array([2,8,5,-6])
B = np.array([7,-4,3,12])
for mat in (A,B) : 
    mat.shape = (2,2)
print('AB \n', np.dot(A,B), '\n BA \n', np.dot(B,A))

AB 
 [[ 38  88]
 [ 17 -92]] 
 BA 
 [[ -6  80]
 [ 66 -48]]

Solving equations with matrices

Let's imagine that we have the following equation :

$$ 2 \begin{bmatrix} x+2 & y+3 \\ 3 & 0 \end{bmatrix} = \begin{bmatrix} 3 & 6 \\ y & z \end{bmatrix}^\top$$

Solving it at hand would look like this : $$\Rightarrow \begin{bmatrix} 2x+4 & 2y+6 \\ 6 & 0 \end{bmatrix} = \begin{bmatrix} 3 & y \\ 6 & z \end{bmatrix}$$ $$\Rightarrow \begin{cases} z = 0 \\ y = -6 \\ x = \frac{-1}{2} \end{cases}$$

  • But we can use SymPy to solve it in the following manner :
In [3]:
x, y, z = sp.symbols('x,y,z') # Initialize our variable symbols
lhs = 2 * sp.Matrix([[x+2, y+3],[3,0]]) # Create the matrix on the left hand side
lhs
Out[3]:
$\displaystyle \left[\begin{matrix}2 x + 4 & 2 y + 6\\6 & 0\end{matrix}\right]$
In [4]:
rhs = sp.transpose(sp.Matrix([[3,6],[y,z]])) # Create the matrix on the right hand side
rhs
Out[4]:
$\displaystyle \left[\begin{matrix}3 & y\\6 & z\end{matrix}\right]$
In [5]:
sp.solve(sp.Eq(lhs,rhs)) # And solve the equality between the two
Out[5]:
{x: -1/2, y: -6, z: 0}

Identity matrix

$$ A = \begin{bmatrix} 2 & 5 & 7 \\ 4 & 8 & 3 \\ 9 & 1 & 4 \\\end{bmatrix} \ et \ I = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\\end{bmatrix}$$

Let's verify that $AI = A$

With NumPy :

In [18]:
A = np.array([2,5,7,4,8,3,9,1,4])
A.shape = (3,3)
I = np.eye(3) # Identity matrix 3x3
print('A.I \n', np.dot(A,I))

A.I 
 [[2. 5. 7.]
 [4. 8. 3.]
 [9. 1. 4.]]

And with SymPy :

In [19]:
A = sp.Matrix(3,3,[2,5,7,4,8,3,9,1,4])
I = sp.eye(3)
A * I
Out[19]:
$\displaystyle \left[\begin{matrix}2 & 5 & 7\\4 & 8 & 3\\9 & 1 & 4\end{matrix}\right]$
In [20]:
sp.Eq(A*I, A) # Verifying that A.I is A
Out[20]:
$\displaystyle \text{True}$

Matrix transposition

With the following matrices, let's verify an important property of matrices and their transpose : $(AB)^\top = B^\top A^\top$ $$ A = \begin{bmatrix} 3 & 2 & 3 \\ 2 & 2 & 1 \\ 2 & 3 & 1 \\\end{bmatrix} \ et \ B = \begin{bmatrix} 2 & -2 \\ -3 & 4 \\ 1 & 1\\\end{bmatrix}$$

SymPy :

In [21]:
A = sp.Matrix(3,3,[3,2,3,2,2,1,2,3,1])
B = sp.Matrix(3,2,[2,-2,-3,4,1,1])
(A*B).T # We can use sp.transpose(A*B) too
Out[21]:
$\displaystyle \left[\begin{matrix}3 & -1 & -4\\5 & 5 & 9\end{matrix}\right]$
In [22]:
B.T * A.T
Out[22]:
$\displaystyle \left[\begin{matrix}3 & -1 & -4\\5 & 5 & 9\end{matrix}\right]$
In [23]:
(A*B).T == B.T * A.T # Let's verify the property
Out[23]:
True

NumPy :

In [24]:
A = np.array([3,2,3,2,2,1,2,3,1])
A.shape = (3,3)
B = np.array([2,-2,-3,4,1,1])
B.shape = (3,2)
print('(AB)T \n', np.dot(A,B).T, '\n BT.AT \n', np.dot(B.T,A.T)) # .T method will give the array transpose 

(AB)T 
 [[ 3 -1 -4]
 [ 5  5  9]] 
 BT.AT 
 [[ 3 -1 -4]
 [ 5  5  9]]
In [26]:
np.array_equal(np.dot(A,B).T, np.dot(B.T,A.T)) # Verifying the equality
Out[26]:
True

Matrix Determinant

  • If we want to calculate the determinant for the follwing matrix : $$ A = \begin{bmatrix} 2 & 3 \\ 4 & -2 \end{bmatrix}$$

With SymPy :

In [13]:
A = sp.Matrix(2,2,[2,3,4,-2])
A.det()
Out[13]:
$\displaystyle -16$

And with NumPy :

In [16]:
A = np.array([[2,3],[4,-2]])
np.linalg.det(A)
Out[16]:
-15.999999999999998
  • Another example : $$ I = \begin{bmatrix} 2 & -2 & 4 \\ 1 & -3 & 2 \\ -2 & 5 & 1 \end{bmatrix}$$

The determinant is :

$$det = +2 \begin{pmatrix} -3 & 2 \\ 5 & 1 \end{pmatrix} - (-2)\begin{pmatrix} 1 & 2 \\ -2 & 1 \end{pmatrix} +4\begin{pmatrix} 1 & -3 \\ -2 & 5 \end{pmatrix}$$$$\Rightarrow det = +2(-13) +2(5) +4(-1) = -20$$

Easily with SymPy :

In [27]:
I = sp.Matrix(3,3,[2,-2,4,1,-3,2,-2,5,1])
I.det()
Out[27]:
$\displaystyle -20$

And with NumPy :

In [28]:
I = np.array([2,-2,4,1,-3,2,-2,5,1])
I.shape = (3,3)
np.linalg.det(I)
Out[28]:
-19.999999999999996

Matrix inverse

  • If we have the following matrix :
$$ A = \begin{bmatrix} 1 & 1 & 1 \\ 1 & -1 & -1 \\ 1 & 1 & -1 \end{bmatrix}$$

the matrix inverse has a property which is that $AA^{−1} = A^{−1}A = I$.

We can verify it with SymPy :

In [29]:
A = sp.Matrix(3,3,[1,1,1,1,-1,-1,1,1,-1])
A_1 = A.inv() # Calculating the inverse
A_1
Out[29]:
$\displaystyle \left[\begin{matrix}\frac{1}{2} & \frac{1}{2} & 0\\0 & - \frac{1}{2} & \frac{1}{2}\\\frac{1}{2} & 0 & - \frac{1}{2}\end{matrix}\right]$
In [30]:
A*A_1 == A_1*A == sp.eye(3) # Verifying the property
Out[30]:
True

And with NumPy, it goes as follows :

In [31]:
A = np.array([1,1,1,1,-1,-1,1,1,-1])
A.shape = (3,3)
A_1 = np.linalg.inv(A)
print('A^-1 \n', A_1)

A^-1 
 [[ 0.5  0.5  0. ]
 [-0.  -0.5  0.5]
 [ 0.5 -0.  -0.5]]
In [32]:
if np.array_equal(np.dot(A,A_1), np.dot(A_1,A)) and np.array_equal(np.dot(A,A_1), np.eye(3)) : 
    print("The property has been verified")

The property has been verified

Solving linear systems of equations

Let's take the following system of equations :

$$x + y + z = 8$$$$ x−y−z = −4$$$$ x + y−z = 6$$

We can solve them in the following manner with NumPy :

In [35]:
a = np.array([[1,1,1],[1,-1,-1],[1,1,-1]])
b = np.array([[8],[-4],[6]])
sol = np.linalg.solve(a,b)
print('The solution is the vector column [x,y,z] \n',sol)

The solution is the vector column [x,y,z] 
 [[2.]
 [5.]
 [1.]]
In [37]:
np.array_equal(np.dot(a,sol),b) # Verifying the solution is true
Out[37]:
True

Eigenvalues and Eigenvectors

  • If we want to calculate the eigenvalues and the eignvectors for the following matrix :
$$ B = \begin{bmatrix} 2 & 16 & 8 \\ 4 & 14 & 8 \\ -8 & -32 & -18 \end{bmatrix}$$
In [39]:
B = np.array([2,16,8,4,14,8,-8,-32,-18])
B.shape = (3,3) # A 3x3 matrix 
eig_vals, eig_vects = np.linalg.eig(B)
print("The Eigenvalues are \n", eig_vals, "\nAnd the Eigenvectors are \n", eig_vects)

The Eigenvalues are 
 [ 2. -2. -2.] 
And the Eigenvectors are 
 [[ 0.40824829 -0.84689344  0.29640505]
 [ 0.40824829 -0.05284196  0.36680613]
 [-0.81649658  0.52913064 -0.88181478]]

Introduction to Graph Theory

Graph creation and structure description

  • Let's write a function that creates a NetworkX graph and gives the degree distribution as well as the adjacency matrix
In [6]:
def graph(nodes, edges, type="simple"):
    """This function creates a graph and sends back the degree distribution as well
    as the sparse matrix. graph can be simple or oriented by using the parameter type.
    nodes : a list of graph vertices.
    edges : a list of tuples, each tuple contains two vertices to represent an edge.
    """

    if type == "simple":
        g = nx.Graph()
    elif type == "oriented":
        g = nx.DiGraph()

    g.add_nodes_from(nodes)  # Add nodes from the list of nodes
    g.add_edges_from(edges)  # Add edges

    a = nx.adjacency_matrix(g)  # The corresponding adjacency matrix

    degree, freq = zip( \
        *Counter(sorted([deg for node, deg in g.degree()])).items()) # Sorting the degrees and counting them

    fig, ax = plt.subplots(1, 3, figsize=(15, 5), dpi=150)
    plt.subplot(131)  # Drawing the graph
    nx.draw(g, with_labels=True)
    plt.title("Graph")
    plt.subplot(132)  # The degree distribution
    plt.bar(degree, freq, width=0.2)
    plt.xlabel("Degrees")
    plt.ylabel("Frequency")
    plt.xticks(nodes)
    plt.title("Degree distribution of the graph")
    plt.subplot(133)  # The sparse matrix
    plt.spy(a, marker='x')
    plt.xticks([i for i in range(len(nodes))], nodes)
    plt.yticks([i for i in range(len(nodes))], nodes)
    plt.title("Sparse Matrix")

    return g, a  # Graph and adjacency matrix
  • Now Let's create a graph that has 4 vertices (1,2,3,4) where each vertex is connnected to the other 3.
In [8]:
vertices = [n for n in range(1, 5)]
g1, a1 = graph(vertices, list(combinations(
    vertices, 2)))  # combinations function creates the list of edges needed

From the degree distribution diagram, we can see that our graph has 4 nodes of degree 3.

  • We can also create an oriented graph :
In [10]:
g2, a2 = graph(vertices, [(1,2),(2,1),(2,3),(2,1),(4,1),(4,2),(4,3)], type = "oriented")

Centrality measures

  • Centrality measures in graph theory are very important because they can indicate to some very important aspects of the vertices/nodes in a graph.
  • Some of the indicators of centrality are : Degree centrality, Eigenvector centrality, Closeness centrality and betweenness centrality.

Now let's write a function that can represent those types of indicators for a certain graph.

In [19]:
def centrality(graph, centr_type, title):
    """This function creates a represntation for a given centrality type and for a given graph.
    graph : networkx graph object
    centr_type : type of centrality on a given graph, for example : nx.degree_centrality(graph)
    title : string for the title of the representation.
    """

    plt.figure(dpi=150)
    pos = nx.spring_layout(graph)
    nodes = nx.draw_networkx_nodes(
        graph,
        pos,
        node_size=250,
        cmap=plt.cm.plasma,
        node_color=list(centr_type.values()
                        ),  # color nodes according to their centrality values
        nodelist=list(centr_type.keys()))
    nodes.set_norm(mcolors.SymLogNorm(linthresh=0.01, linscale=1, base=10))

    labels = nx.draw_networkx_labels(graph, pos)
    edges = nx.draw_networkx_edges(graph, pos)

    plt.title(title)
    plt.colorbar(nodes)
    plt.axis('off')

Now let's initialize a random graph so we can try our function on.

In [13]:
ver = [n for n in range(1, 15)] # nodes from 1 to 14
edges = [(l1,l2) for l1,l2 in map(lambda n : (n,8), [n for n in range(11,15)])] + \
    [(l1,l2) for l1,l2 in map(lambda n : (n,6), [n for n in range(1,6)])] + \
    [(7,6), (8,7),(9,7), (10,7)] # list of edges
g3, a3 = graph(ver, edges)

Degree centrality

Let's check the degree centrality for our 14 vertices.

In [22]:
for node in g3.nodes():
    print(node, nx.degree_centrality(g3)[node])

1 0.07692307692307693
2 0.07692307692307693
3 0.07692307692307693
4 0.07692307692307693
5 0.07692307692307693
6 0.46153846153846156
7 0.3076923076923077
8 0.38461538461538464
9 0.07692307692307693
10 0.07692307692307693
11 0.07692307692307693
12 0.07692307692307693
13 0.07692307692307693
14 0.07692307692307693
In [23]:
centrality(g3, nx.degree_centrality(g3), 'Degree Centrality')

Eigenvector centrality

In [24]:
for node in g3.nodes():
    print(node, nx.eigenvector_centrality(g3)[node])

1 0.2005897407903768
2 0.2005897407903768
3 0.2005897407903768
4 0.2005897407903768
5 0.2005897407903768
6 0.543357022804172
7 0.4688940252630566
8 0.3805718175986731
9 0.17310146940310608
10 0.17310146940310608
11 0.1404965947281777
12 0.1404965947281777
13 0.1404965947281777
14 0.1404965947281777
In [25]:
centrality(g3, nx.eigenvector_centrality(g3), 'Eigenvector centrality')

Closeness centrality

In [26]:
for node in g3.nodes():
    print(node, nx.closeness_centrality(g3)[node])

1 0.3611111111111111
2 0.3611111111111111
3 0.3611111111111111
4 0.3611111111111111
5 0.3611111111111111
6 0.5416666666666666
7 0.5909090909090909
8 0.5
9 0.38235294117647056
10 0.38235294117647056
11 0.34210526315789475
12 0.34210526315789475
13 0.34210526315789475
14 0.34210526315789475
In [28]:
centrality(g3, nx.closeness_centrality(g3), 'Closeness centrality')

Betweennes centrality

In [29]:
for node in g3.nodes():
    print(node, nx.betweenness_centrality(g3)[node])

1 0.0
2 0.0
3 0.0
4 0.0
5 0.0
6 0.641025641025641
7 0.6794871794871795
8 0.5384615384615384
9 0.0
10 0.0
11 0.0
12 0.0
13 0.0
14 0.0
In [30]:
centrality(g3, nx.betweenness_centrality(g3),'Betweenness centrality')