Mathematics Archives

Instructions on Transformer for people outside NLP field, but with examples of NLP

December 30, 2020/in Artificial Intelligence, Data Mining, Data Science, Deep Learning, Machine Learning, Main Category, Mathematics, Natural Language Processing, Predictive Analytics/by Yasuto Tamura

I found it quite difficult to explain mathematical details of long short-term memory (LSTM) in my previous article series. But when I was studying LSTM, a new promising algorithm was already attracting attentions. The algorithm is named Transformer. Its algorithm was a first announced in a paper named “Attention Is All You Need,” and it outperformed conventional translation algorithms with lower computational costs.

In this article series, I am going to provide explanations on minimum prerequisites for understanding deep learning in NLP (natural language process) tasks, but NLP is not the main focus of this article series, and actually I do not study in NLP field. I think Transformer is going to be a new major model of deep learning as well as CNN or RNN, and the model is now being applied in various fields.

Even though Transformer is going to be a very general deep learning model, I still believe it would be an effective way to understand Transformer with some NLP because language is a good topic we have in common. Unlike my previous article series, in which I tried to explain theoretical side of RNN as precisely as possible, in this article I am going to focus on practical stuff with my toy implementations of NLP tasks, largely based on Tensorflow official tutorial. But still I will do my best to make it as straightforward as possible to understand the architecture of Transformer with various original figures.

This series is going to be composed of the articles below.

On the difficulty of language: prerequisites for NLP with Transformer
Seq2seq model and attention mechanism: a backbone of NLP with deep learning
Multi-head attention: the key component of Transformer
Positional encoding, residual connections, padding masks: covering the rest of Transformer components
How to make a toy English-German translator with multi-head attention heat maps: the overall architecture of Transformer
Transformer in image processing (Coming soon)

If you are in the field and can read the codes in the official tutorial with no questions, this article series is not for you, but if you want to see how a Transformer works but do not want to go too much into details of NLP, this article would be for you.

Rethinking linear algebra: visualizing linear transformations and eigen vectors

October 27, 2020/in Data Mining, Data Science, Machine Learning, Main Category, Mathematics/by Yasuto Tamura

In terms of calculation processes of Principal Component Analysis (PCA) or Linear Discriminant Analysis (LDA), which are the dimension reduction techniques I am going to explain in the following articles, diagonalization is what they are all about. Throughout this article, I would like you to have richer insight into diagonalization in order to prepare for understanding those basic dimension reduction techniques.

When our professor started a lecture on the last chapter of our textbook on linear algebra, he said “It is no exaggeration to say that everything we have studied is for this ‘diagonalization.'” Until then we had to write tons of numerical matrices and vectors all over our notebooks, calculating those products, adding their rows or columns to other rows or columns, sometimes transposing the matrices, calculating their determinants.

It was like the scene in “The Karate Kid,” where the protagonist finally understood the profound meaning behind the prolonged and boring “wax on, wax off” training given by Miyagi (or “jacket on, jacket off” training given by Jackie Chan). We had finally understood why we had been doing those seemingly endless calculations.

Source: http://thinkbedoleadership.com/secret-success-wax-wax-off/

But usually you can do those calculations easily with functions in the Numpy library. Unlike Japanese college freshmen, I bet you are too busy to reopen textbooks on linear algebra to refresh your mathematics. Thus I am going to provide less mathematical and more intuitive explanation of diagonalization in this article.

1, The mainstream ways of explaining diagonalization.

*The statements below are very rough for mathematical topics, but I am going to give priority to offering more visual understanding on linear algebra in this article. For further understanding, please refer to textbooks on linear algebra. If you would like to have minimum understandings on linear algebra needed for machine learning, I recommend the Appendix C of Pattern Recognition and Machine Learning by C. M. Bishop.

In most textbooks on linear algebra, the explanations on dioagonalization is like this (if you are not sure what diagonalization is or if you are allergic to mathematics, you do not have to read this seriously):

Let $V (dimV = D)$ be a vector space and let $T_A : V \rightarrow V$ be a mapping of $V$ into itself, defined as $T_A(v) = A \cdot \boldsymbol{v}$ , where $A$ is a $D\times D$ matrix and $\boldsymbol{v}$ is $D$ dimensional vector. An element $\boldsymbol{v} \in V$ is called an eigen vector if there exists a number $\lambda$ such that $A \cdot \boldsymbol{v}= \lambda \cdot \boldsymbol{v}$ and $\boldsymbol{v} \neq \boldsymbol{0}$ . In this case $\lambda$ is uniquely determined and is called an eigen value of $A$ belonging to the eigen vector $\boldsymbol{v}$ .

Any matrix $A$ has $D$ eigen values $\lambda_{i}$ , belonging to $\boldsymbol{v}_{i} (i=1, 2, …., D)$ . If $\boldsymbol{v}_{i}$ is basis of the vector space $V$ , then $A$ is diagonalizable.

When $A$ is diagonalizable, with $D \times D$ matrices $P = (\boldsymbol{v}_{1}, \dots, \boldsymbol{v}_{D})$ , whose column vectors are eigen vectors $\boldsymbol{v}_{i} (i=1, 2, …., D)$ , the following equation holds: $P^{-1}AP = \Lambda$ , where $\Lambda = diag(\lambda_{1}, \dots, \lambda_{D})= \begin{pmatrix} \lambda_{1} & 0& \ldots &0\\ 0 & \lambda_{2} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_{D} \end{pmatrix}$ .

And when $A$ is diagonalizable, you can diagonalize $A$ as below.

Most textbooks keep explaining these type of stuff, but I have to say they lack efforts to make it understandable to readers with low mathematical literacy like me. Especially if you have to apply the idea to data science field, I believe you need more visual understanding of diagonalization. Therefore instead of just explaining the definitions and theorems, I would like to take a different approach. But in order to understand them in more intuitive ways, we first have to rethink waht linear transformation $T_A$ means in more visible ways.

2, Linear transformations

Even though I did my best to make this article understandable to as little prerequisite knowledge, you at least have to understand linear transformation of numerical vectors and with matrices. Linear transformation is nothing difficult, and in this article I am going to use only 2 or 3 dimensional numerical vectors or square matrices. You can calculate linear transformation of $\boldsymbol{v}$ by A as equations in the figure. In other words, $\boldsymbol{u}$ is a vector transformed by $A$ .

*I am not going to use the term “linear transformation” in a precise way in the context of linear algebra. In this article or in the context of data science or machine learning, “linear transformation” for the most part means products of matrices or vectors.

*Forward/back propagation of deep learning is mainly composed of this linear transformation. You keep linearly transforming input vectors, frequently transforming them with activation functions, which are for the most part not linear transformation.

As you can see in the equations above, linear transformation with A transforms a vector to another vector. Assume that you have an original vector $\boldsymbol{v}$ in grey and that the vector $\boldsymbol{u}$ in pink is the transformed $\boldsymbol{v}$ by A is. If you subtract $\boldsymbol{v}$ from $\boldsymbol{u}$ , you can get a displacement vector, which I displayed in purple. A displacement vector means the transition from a vector to another vector.

Let’s calculate the displacement vector with more vectors $\boldsymbol{v}$ . Assume that $A =\begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}$ , and I prepared several grid vectors $\boldsymbol{v}$ in grey as you can see in the figure below. If you transform those grey grid points with $A$ , they are mapped into the vectors $\boldsymbol{u}$ in pink. With those vectors in grey or pink, you can calculate the their displacement vectors $\boldsymbol{u} = \boldsymbol{v}$ in purple.

I think you noticed that the displacement vectors in the figure above have some tendencies. In order to see that more clearly, let’s calculate displacement vectors with several matrices $A$ and more grid points. Assume that you have three $2 \times 2$ square matrices $A_1 =\begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}, A_2 =\begin{pmatrix} 3 & 1 \\ -1 & 1 \end{pmatrix}, A_3 =\begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix}$ , and I plotted displace vectors made by the matrices respectively in the figure below.

I think you noticed some characteristics of the displacement vectors made by those linear transformations: the vectors are swirling and many of them seem to be oriented in certain directions. To be exact, some displacement vectors have extend in the same directions as some of original vectors in grey. That means linear transformation by $A$ did not change the direction of the original vector $\boldsymbol{v}$ , and the unchanged vectors are called eigen vectors. Real eigen vectors of each A are displayed as arrows in yellow in the figure above. But when it comes to $A_3$ , the matrix does not have any real eigan values.

In linear algebra, depending on the type matrices $A$ , you have consider various cases such as whether the matrices have real or imaginary eigen values, whether the matrices are diagonalizable, whether the eigen vectors are orthogonal, or whether they are unit vectors. But those topics are out of the scope of this article series, so please refer to textbooks on linear algebra if you are interested.

Luckily, however, in terms of PCA or LDA, you only have to consider a type of matrices named positive semidefinite matrices, which $A_1$ is classified to, and I am going to explain positive semidefinite matrices in the fourth section.

3, Eigen vectors as coordinate system

Source: Ian Stewart, “Professor Stewart’s Cabinet of Mathematical Curiosities,” (2008), Basic Books

Let me take Fibonacci numbers as an example to briefly see why diagonalization is useful. Fibonacci is sequence is quite simple and it is often explained using an example of pairs of rabbits increasing generation by generation. Let $a_n (n=0, 1, 2, …)$ be the number of pairs of grown up rabbits in the $n^{th}$ generation. One pair of grown up rabbits produce one pair of young rabbit The concrete values of $a_n$ are $a_0 = 0$ , $a_1 = 1$ , $a_2=1$ , $a_3=2$ , $a_4=3$ , $a_5=5$ , $a_6=8$ , $a_7=13, \dots$ . Assume that $A =\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$ and that $\begin{pmatrix} a_1 \\ a_0 \end{pmatrix} =\begin{pmatrix} 1 \\ 0 \end{pmatrix}$ , then you can calculate the number of the pairs of grown up rabbits in the next generation with the following recurrence relation. $\begin{pmatrix} a_{n+1} \\ a_{n} \end{pmatrix}=\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \cdot \begin{pmatrix} a_{n+1} \\ a_{n} \end{pmatrix}$ .Let $\boldsymbol{a}_n$ be $\begin{pmatrix} a_{n+1} \\ a_{n} \end{pmatrix}$ , then the recurrence relation can be written as $\boldsymbol{a}_{n+1} = A \boldsymbol{a}_n$ , and the transition of $\boldsymbol{a}_n$ are like purple arrows in the figure below. It seems that the changes of the purple arrows are irregular if you look at the plots in normal coordinate.

Assume that $\lambda _1, \lambda_2 (\lambda _1< \lambda_2)$ are eigen values of $A$ , and $\boldsymbol{v}_1, \boldsymbol{v}_2$ are eigen vectors belonging to them respectively. Also let $\alpha, \beta$ scalars such that $\begin{pmatrix} a_{1} \\ a_{0} \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} = \alpha \boldsymbol{v}_1 + \beta \boldsymbol{v}_2$ . According to the definition of eigen values and eigen vectors belonging to them, the following two equations hold: $A\boldsymbol{v}_1 = \lambda_1 \boldsymbol{v}_1, A\boldsymbol{v}_2 = \lambda_2 \boldsymbol{v}_2$ . If you calculate $\boldsymbol{a}_1$ is, using eigen vectors of $A$ , $\boldsymbol{a}_1 = A\boldsymbol{a}_0 = A (\alpha \boldsymbol{v}_1 + \beta \boldsymbol{v}_2) = \alpha\lambda _1 \boldsymbol{v}_1 + \beta \lambda_2 \boldsymbol{v}_2$ . In the same way, $\boldsymbol{a}_2 = A\boldsymbol{a}_1 = A (\alpha\lambda _1 \boldsymbol{v}_1 + \beta \lambda_2 \boldsymbol{v}_2) = \alpha\lambda _{1}^{2} \boldsymbol{v}_1 + \beta \lambda_{2}^{2} \boldsymbol{v}_2$ , and $\boldsymbol{a}_3 = A\boldsymbol{a}_2 = A (\alpha\lambda _{1}^{2} \boldsymbol{v}_1 + \beta \lambda_{2}^{2} \boldsymbol{v}_2) = \alpha\lambda _{1}^{3} \boldsymbol{v}_1 + \beta \lambda_{2}^{3} \boldsymbol{v}_2$ . These equations show that in coordinate system made by eigen vectors of $A$ , linear transformation by $A$ is easily done by just multiplying eigen values with each eigen vector. Compared to the graph of Fibonacci numbers above, in the figure below you can see that in coordinate system made by eigen vectors the plots changes more systematically generation by generation.

In coordinate system made by eigen vectors of square matrices, the linear transformations by the matrices can be much more straightforward, and this is one powerful strength of eigen vectors.

*I do not major in mathematics, so I am not 100% sure, but vectors in linear algebra have more abstract meanings and various things in mathematics can be vectors, even though in machine learning or data science we mainly use numerical vectors with more concrete elements. We can also say that matrices are a kind of maps. That is just like, at leas in my impression, even though a real town is composed of various components such as houses, smooth or bumpy roads, you can simplify its structure with simple orthogonal lines, like the map of Manhattan. But if you know what the town actually looks like, you do not have to follow the zigzag path on the map.

4, Eigen vectors of positive semidefinite matrices

In the second section of this article I told you that, even though you have to consider various elements when you discuss general diagonalization, in terms of PCA and LDA we mainly use only a type of matrices named positive semidefinite matrices. Let $A$ be a $D \times D$ square matrix. If $\boldsymbol{x}^T A \boldsymbol{x} \geq 0$ for all values of the vector $\boldsymbol{x}$ , the $A$ is said to be a positive semidefinite matrix. And also it is known that $A$ being a semidefinite matrix is equivalent to $\lambda _{i} \geq 0$ for all the eigen values $\lambda_i (i=1, \dots , D)$ .

*I think most people first learn a type of matrices called positive definite matrices. Let $A$ be a $D \times D$ square matrix. If $\boldsymbol{x}^T A \boldsymbol{x} > 0$ for all values of the vector $\boldsymbol{x}$ , the $A$ is said to be a positive definite matrix. You have to keep it in mind that even if all the elements of $A$ are positive, $A$ is not necessarly positive definite/semidefinite.

Just as we did in the second section of this article, let’s visualize displacement vectors made by linear transformation with a $3 \times 3$ square positive semidefinite matrix $A$ .

*In fact $A_1 =\begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}$ , whose linear transformation I visualized the second section, is also positive semidefinite.

Let’s visualize linear transformations by a positive definite matrix $A = \frac{1}{50} \begin{pmatrix} 60.45 & 33.63 & 46.29 \\33.63 & 68.49 & 50.93 \\ 46.29 & 50.93 & 53.61 \end{pmatrix}$ . I visualized the displacement vectors made by the $A$ just as the same way as in the second section of this article. The result is as below, and you can see that, as well as the displacement vectors made by $A_1$ , the three dimensional displacement vectors below are swirling and extending in three directions, in the directions of the three orthogonal eigen vectors $\boldsymbol{v}_1, \boldsymbol{v}_2$ , and $\boldsymbol{v}_3$ .

*It might seem like a weird choice of a matrix, but you are going to see why in the next article.

You might have already noticed $A_1 =\begin{pmatrix} 3 & 1 \\ 1 & 2 \end{pmatrix}$ and $A = \frac{1}{50} \begin{pmatrix} 60.45 & 33.63 & 46.29 \\33.63 & 68.49 & 50.93 \\ 46.29 & 50.93 & 53.61 \end{pmatrix}$ are both symmetric matrices and that their elements are all real values, and that their diagonal elements are all positive values. Super importantly, when all the elements of a $D \times D$ symmetric matrix $A$ are real values and its eigen values are $\lambda_{i} (i=1, \dots , D)$ , there exist orthonormal matrices $U$ such that $U^{-1}AU = \Lambda$ , where $\Lambda = diag(\lambda_{1}, \dots , \lambda_{D})$ .

*The title of this section might be misleading, but please keep it in mind that positive definite/semidefinite matrices are not necessarily real symmetric matrices. And real symmetric vectors are not necessarily positive definite/semidefinite matrices.

5, Orthonormal matrices and rotation of vectors

In this section I am gong to explain orthonormal matrices, as known as rotation matrices. If a $D\times D$ matrix $U$ is an orthonormal matrix, column vectors of $U$ are orthonormal, which means $U = (\boldsymbol{u}_1 \dots \boldsymbol{u}_D)$ , where $\begin{cases} \boldsymbol{u}_{i}^{T}\boldsymbol{u}_{j} = 1 \quad (i = j) \\ \boldsymbol{u}_{i}^{T}\boldsymbol{u}_{j} = 0 \quad (i\neq j) \end{cases}$ . In other words column vectors $\boldsymbol{u}_{i}$ form an orthonormal coordinate system.

Orthonormal matrices $U$ have several important matrices, and one of them is $U^{-1} = U^{T}$ . Combining this fact with what I have told you so far, you we can reach one conclusion that you can orthogonalize a real symmetric matrix $A$ as $U^{T}AU = \Lambda$ . This is known as spectral decomposition or singular value decomposition.

Another important property of $U$ is that $U^{T}$ is also orthonormal. In other words, assume $U$ is orthonormal and that $U = (\boldsymbol{u}_1 \dots \boldsymbol{u}_D) = \begin{pmatrix} -\boldsymbol{v_1}^{T}- \\ \vdots \\ -\boldsymbol{v_D}^{T}- \end{pmatrix}$ , $(\boldsymbol{v}_1 \dots \boldsymbol{v}_D)$ also forms a orthonormal coordinate system.

…It seems things are getting too mathematical and abstract (for me), thus for now I am going to wrap up what I have explained in this article .

We have seen

Numerical matrices linearly transform vectors.
Certain linear transformations do not change the direction of vectors in certain directions, which are called eigen vectors.
Making use of eigen vectors, you can form new coordinate system which can describe the linear transformations in a more straightforward way.
You can diagonalize a real symmetric matrix $A$ with an orthonormal matrix $U$ .

Of our current interest is what kind of linear transformation the real symmetric positive definite matrix enables. I am going to explain why the purple vectors in the figure above is swirling like that in the upcoming articles. Before that, however, we are going to see one application of what we have seen in this article, on dimension reduction. To be concrete the next article is going to be about principal component analysis (PCA), which is very important in many fields.

*In short, the orthonormal matrix $U$ I mentioned above enables rotation of matrix, and the diagonal matrix $diag(\lambda_1, \dots, \lambda_D)$ expands or contracts vectors along each axis. I am going to explain that more precisely in the upcoming articles.

* I make study materials on machine learning, sponsored by DATANOMIQ. I do my best to make my content as straightforward but as precise as possible. I include all of my reference sources. If you notice any mistakes in my materials, including grammatical errors, please let me know (email: yasuto.tamura@datanomiq.de). And if you have any advice for making my materials more understandable to learners, I would appreciate hearing it.

*I attatched the codes I used to make the figures in this article. You can just copy, paste, and run, sometimes installing necessary libraries.

import matplotlib.pyplot as plt 
import numpy as np
import matplotlib.patches as mpatches

T_A = np.array([[1, 1], 
                [1, 0]])

total_step = 5
x = np.zeros((total_step, 2))

x[0]  = np.array([1, 0])

for i in range(total_step - 1):
    x[i + 1] = np.dot(T_A, x[i])

eigen_values, eigen_vectors = np.linalg.eig(T_A)
idx = eigen_values.argsort()[::-1]   
eigen_values = eigen_values[idx]
eigen_vectors = eigen_vectors[:,idx]
for i in range(len(eigen_vectors)):
    if(eigen_vectors.T[i][0] < 0):
        eigen_vectors.T[i] = - eigen_vectors.T[i]   

v_initial = x[0]
v_coefficients = np.zeros((total_step , 2))
v_coefficients[0] = np.dot(v_initial ,  np.linalg.inv(eigen_vectors.T)) 

for i in range(total_step-1):
    v_coefficients[i + 1] = v_coefficients[i] * eigen_values 

for i in range(total_step):
    v_1_list[i+1] = v_coefficients.T[0][i]*eigen_vectors.T[0]
    v_2_list[i+1] = v_coefficients.T[1][i]*eigen_vectors.T[1]

plt.figure(figsize=(20, 15))
fontsize = 20
small_shift = 0.2

plt.plot(x[:, 0], x[:, 1], marker='o', linestyle='none', markersize=10, color='black')

plt.arrow(0, 0, eigen_vectors.T[0][0], eigen_vectors.T[0][1], width=0.05, head_width=0.2, color='orange')
plt.arrow(0, 0, eigen_vectors.T[1][0], eigen_vectors.T[1][1], width=0.05, head_width=0.2, color='orange')

plt.text(eigen_vectors.T[0][0], eigen_vectors.T[0][1]+small_shift, r'', va='center',ha='right', fontsize=fontsize + 10)
plt.text(eigen_vectors.T[1][0] - small_shift, eigen_vectors.T[1][1],r'', va='center',ha='right', fontsize=fontsize + 10)

for i in range(total_step): 
    
    plt.arrow(0, 0, v_1_list[i+1][0], v_1_list[i+1][1], head_width=0.05, color='darkviolet', length_includes_head=True)
    plt.arrow(0, 0, v_2_list[i+1][0], v_2_list[i+1][1], head_width=0.05, color='darkviolet', length_includes_head=True)
    
    plt.text(v_1_list[i+1][0] + 2*small_shift , v_1_list[i+1][1]-2*small_shift,r''.format(1,i+1, 1),va='center',ha='right', fontsize=fontsize)
    plt.text(v_2_list[i+1][0]-0.1, v_2_list[i+1][1],r''.format(2, i+1, 2),va='center',ha='right', fontsize=fontsize)

    plt.arrow(v_1_list[i+1][0], v_1_list[i+1][1], v_2_list[i+1][0], v_2_list[i+1][1], head_width=0, color='black', linestyle='--', length_includes_head=True)
    plt.arrow(v_2_list[i+1][0], v_2_list[i+1][1], v_1_list[i+1][0], v_1_list[i+1][1], head_width=0, color='black', linestyle='--', length_includes_head=True)
    
orange_patch = mpatches.Patch(color='orange', label='Eigen vectors')
purple_patch = mpatches.Patch(color='darkviolet', label='Scalar multiples of the eigen vectors')
plt.legend(handles=[orange_patch, purple_patch], fontsize=25, loc='lower right')

for i in range(total_step):
    plt.text(x[i][0]+0.1, x[i][1]-0.05, r'n={0}'.format(i), fontsize=20)

plt.grid(True)
plt.ylabel("", fontsize=20)
plt.xlabel("", fontsize=20)
plt.title("Fibonacci sequence and its eigen space", fontsize=30)
#plt.savefig("Fibonacci_eigen_space.png")
plt.show()

import matplotlib.pyplot as plt 
import numpy as np 
import matplotlib.patches as mpatches

const_range = 10

X = np.arange(-const_range, const_range + 1, 1)
Y = np.arange(-const_range, const_range + 1, 1)
U_x, U_y = np.meshgrid(X, Y)

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

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

T_A_2 = np.array([[1, -1], 
                [1, 1]])

T_A_list = np.array((T_A_0, T_A_1, T_A_2))

const_range = 5
plt.figure(figsize=(30, 10))
plt.subplots_adjust(wspace=0.1)
labels = ["Grids", "Displacement vectors made by A", "Real eigen vectors of A"]  
title_list = [r" has two different real eigen vectors.", r" has two identical real unit eigen vectors.",  r" has only imaginary eigen vectors."]
for idx in range(len(T_A_list)): 
    
    eigen_values, eigen_vectors = np.linalg.eig(T_A_list[idx])
    sorted_idx = eigen_values.argsort()[::-1]   
    eigen_values = eigen_values[sorted_idx]
    eigen_vectors = eigen_vectors[:,sorted_idx]
    eigen_vectors = eigen_vectors.astype(float)
        
    for i in range(len(eigen_vectors)):
        if(eigen_vectors.T[i][0] < 0):
            eigen_vectors.T[i] = - eigen_vectors.T[i]
        

    X = np.arange(-const_range, const_range + 1, 1)
    Y = np.arange(-const_range, const_range + 1, 1)
    U_x, U_y = np.meshgrid(X, Y)

    V_x = np.zeros((len(U_x), len(U_y)))
    V_y = np.zeros((len(U_x), len(U_y)))

    temp_vec = np.zeros((1, 2))

    W_x = np.zeros((len(U_x), len(U_y)))
    W_y = np.zeros((len(U_x), len(U_y)))

    plt.subplot(1, 3, idx + 1)


    for i in range(len(U_x)):
        for j in range(len(U_y)):
            temp_vec[0][0] = U_x[i][j]
            temp_vec[0][1] = U_y[i][j]
        
            temp_vec[0] = np.dot(T_A_list[idx], temp_vec[0])
        
            V_x[i][j] = temp_vec[0][0]
            V_y[i][j] = temp_vec[0][1]
        
            W_x[i][j] = V_x[i][j] - U_x[i][j]
            W_y[i][j] = V_y[i][j] - U_y[i][j]
            #plt.arrow(0, 0, V_x[i][j], V_y[i][j], head_width=0.1, color='red')
            plt.arrow(0, 0, U_x[i][j], U_y[i][j], head_width=0.3, color='dimgrey', label=labels[0])
            plt.arrow(U_x[i][j], U_y[i][j], W_x[i][j], W_y[i][j], head_width=0.3, color='darkviolet', label=labels[1])
            
            range_const = 20
            plt.xlim([-range_const, range_const])
            plt.ylim([-range_const, range_const])
            plt.title(title_list[idx], fontsize=25)
            
            if idx==2:
                continue
 
            plt.arrow(0, 0, eigen_vectors.T[0][0]*10, eigen_vectors.T[0][1]*10, head_width=1, color='orange', label=labels[2])
            plt.arrow(0, 0, eigen_vectors.T[1][0]*10, eigen_vectors.T[1][1]*10, head_width=1, color='orange', label=labels[2])

grey_patch = mpatches.Patch(color='grey', label='Grids')
purple_patch = mpatches.Patch(color='darkviolet', label='Displacement vectors made by A')
yellow_patch = mpatches.Patch(color='gold', label='Real eigen vectors of A')
plt.legend(handles=[grey_patch, purple_patch, yellow_patch], fontsize=25, loc='lower right', bbox_to_anchor=(-0.1, -.35))
#plt.savefig("linear_transformation.png")
plt.show()

import numpy as np 
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.proj3d import proj_transform
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib.text import Annotation
from matplotlib.patches import FancyArrowPatch
import matplotlib.patches as mpatches

class Annotation3D(Annotation):
    def __init__(self, text, xyz, *args, **kwargs):
        super().__init__(text, xy=(0,0), *args, **kwargs)
        self._xyz = xyz

    def draw(self, renderer):
        x2, y2, z2 = proj_transform(*self._xyz, renderer.M)
        self.xy=(x2,y2)
        super().draw(renderer)

def _annotate3D(ax,text, xyz, *args, **kwargs):
    '''Add anotation `text` to an `Axes3d` instance.'''

    annotation= Annotation3D(text, xyz, *args, **kwargs)
    ax.add_artist(annotation)

setattr(Axes3D,'annotate3D',_annotate3D)

class Arrow3D(FancyArrowPatch):
    def __init__(self, x, y, z, dx, dy, dz, *args, **kwargs):
        super().__init__((0,0), (0,0), *args, **kwargs)
        self._xyz = (x,y,z)
        self._dxdydz = (dx,dy,dz)

    def draw(self, renderer):
        x1,y1,z1 = self._xyz
        dx,dy,dz = self._dxdydz
        x2,y2,z2 = (x1+dx,y1+dy,z1+dz)

        xs, ys, zs = proj_transform((x1,x2),(y1,y2),(z1,z2), renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        super().draw(renderer)

def _arrow3D(ax, x, y, z, dx, dy, dz, *args, **kwargs):
    '''Add an 3d arrow to an `Axes3D` instance.'''

    arrow = Arrow3D(x, y, z, dx, dy, dz, *args, **kwargs)
    ax.add_artist(arrow)

setattr(Axes3D,'arrow3D',_arrow3D)

T_A = np.array([[60.45, 33.63, 46.29], 
                [33.63, 68.49, 50.93], 
                [46.29, 50.93, 53.61]])

T_A = T_A/50
const_range = 2


X = np.arange(-const_range, const_range + 1, 1)
Y = np.arange(-const_range, const_range + 1, 1)
Z = np.arange(-const_range, const_range + 1, 1)

U_x, U_y, U_z = np.meshgrid(X, Y, Z)

V_x = np.zeros((len(U_x), len(U_y), len(U_z)))
V_y = np.zeros((len(U_x), len(U_y), len(U_z)))
V_z = np.zeros((len(U_x), len(U_y), len(U_z)))

temp_vec = np.zeros((1, 3))

W_x = np.zeros((len(U_x), len(U_y), len(U_z)))
W_y = np.zeros((len(U_x), len(U_y), len(U_z)))
W_z = np.zeros((len(U_x), len(U_y), len(U_z)))

eigen_values, eigen_vectors = np.linalg.eig(T_A)
sorted_idx = eigen_values.argsort()[::-1]   
eigen_values = eigen_values[sorted_idx]
eigen_vectors = eigen_vectors[:,sorted_idx]
eigen_vectors = eigen_vectors.astype(float)

fig = plt.figure(figsize=(15, 15))
ax = fig.add_subplot(111, projection='3d')
grid_range = const_range + 5
ax.set_xlim(-grid_range, grid_range)
ax.set_ylim(-grid_range, grid_range)
ax.set_zlim(-grid_range, grid_range)

eigen_values, eigen_vectors = np.linalg.eig(T_A)
sorted_idx = eigen_values.argsort()[::-1]   
eigen_values = eigen_values[sorted_idx]
eigen_vectors = eigen_vectors[:,sorted_idx]
eigen_vectors = eigen_vectors.astype(float)
     
    
for i in range(len(eigen_vectors)):
    if(eigen_vectors.T[i][0] < 0):
        eigen_vectors.T[i] = - eigen_vectors.T[i]

for i in range(len(U_x)):
    for j in range(len(U_x)):
        for k in range(len(U_x)):
            temp_vec[0][0] = U_x[i][j][k]
            temp_vec[0][1] = U_y[i][j][k]
            temp_vec[0][2] = U_z[i][j][k]
        
            temp_vec[0] = np.dot(T_A, temp_vec[0])
        
            V_x[i][j][k] = temp_vec[0][0]
            V_y[i][j][k] = temp_vec[0][1]
            V_z[i][j][k] = temp_vec[0][2]
        
            W_x[i][j][k] = V_x[i][j][k] - U_x[i][j][k]
            W_y[i][j][k] = V_y[i][j][k] - U_y[i][j][k]
            W_z[i][j][k] = V_z[i][j][k] - U_z[i][j][k]
            ax.arrow3D(0, 0, 0, 
                       U_x[i][j][k], U_y[i][j][k], U_z[i][j][k], 
                       mutation_scale=10, arrowstyle="-|>", fc='dimgrey', ec='dimgrey')
            #ax.arrow3D(0, 0, 0, 
            #          V_x[i][j][k], V_y[i][j][k], V_z[i][j][k], 
            #           mutation_scale=10, arrowstyle="-|>", fc='red', ec='red')
            ax.arrow3D(U_x[i][j][k], U_y[i][j][k], U_z[i][j][k], 
                       W_x[i][j][k], W_y[i][j][k], W_z[i][j][k],
                       mutation_scale=10, arrowstyle="-|>", fc='darkviolet', ec='darkviolet')
            
ax.arrow3D(0, 0, 0, eigen_vectors.T[0][0]*10, eigen_vectors.T[0][1]*10, eigen_vectors.T[0][2]*10,
                       mutation_scale=10,  arrowstyle="-|>", fc='orange', ec='orange')
ax.arrow3D(0, 0, 0, eigen_vectors.T[1][0]*10, eigen_vectors.T[1][1]*10, eigen_vectors.T[1][2]*10,
                       mutation_scale=10, arrowstyle="-|>", fc='orange', ec='orange')
ax.arrow3D(0, 0, 0, eigen_vectors.T[2][0]*10, eigen_vectors.T[2][1]*10, eigen_vectors.T[2][2]*10,
                       mutation_scale=10, arrowstyle="-|>", fc='orange', ec='orange')

ax.text(eigen_vectors.T[0][0]*8 , eigen_vectors.T[0][1]*8, eigen_vectors.T[0][2]*8+1, r'', fontsize=20)
ax.text(eigen_vectors.T[1][0]*8 , eigen_vectors.T[1][1]*8, eigen_vectors.T[1][2]*8, r'', fontsize=20)
ax.text(eigen_vectors.T[2][0]*8 , eigen_vectors.T[2][1]*8, eigen_vectors.T[2][2]*8, r'', fontsize=20)


grey_patch = mpatches.Patch(color='grey', label='Grids')
orange_patch = mpatches.Patch(color='orange', label='Orthogonal eigen vectors of A')
purple_patch = mpatches.Patch(color='darkviolet', label='Displacement vectors made by A')
plt.legend(handles=[grey_patch, orange_patch, purple_patch], fontsize=20, loc='lower right')

ax.set_xlabel(r'', fontsize=25)
ax.set_ylabel(r'', fontsize=25)
ax.set_zlabel(r'', fontsize=25)
#plt.savefig("symmetric_positive_definite_visualizaiton.png")
plt.show()

K Nearest Neighbour For Supervised Learning

October 10, 2020/in Artificial Intelligence, Data Science, Machine Learning, Main Category, Mathematics/by Ram Tavva

K-Nearest Neighbour (KNN) Algorithms is an easy-to-implement & advanced level supervised machine learning algorithm used for both – classification as well as regression problems. However, you can see a wide of its applications in classification problems across various industries.

If you’ve been shopping a lot in e-commerce sites like Amazon, Flipkart, Myntra, or love watching web series over Netflix and Amazon Prime, one common thing you’ve always noticed, and that is recommendations.

Are you wondering how they recommend you following your choice? They use KNN Supervised Learning to find out what you may need the next when you’re buying and recommend you with a few more products.

Imagine you’re looking for an iPhone to purchase. When you scroll down a little, you see some iPhone cases, tempered glasses – saying, “People who purchased an iPhone have also purchased these items. The same applies to Netflix and Amazon Prime. When you finished a show or a series, they give you recommendations of the same genre. And do it all using KNN supervised learning and classify the items for the best user experience.

Advantages Of KNN

Quickest Calculation Time
Simple Algorithms
High Accuracy
Versatile – best use for Regression and Classification.
Doesn’t make any assumptions about data.

Where KNN Are Mostly Used

Simple Recommendation Models
Image Recognition Technology
Decision-Making Models
Calculating Credit Rating

Choosing The Right Value For K

To choose the right value of K, you have to run KNN algorithms several times with different values of K and select the value of K, which reduces the number of errors you’ve come across and come out as the most stable value for K.

Your Step-By-Step Guide For Choosing The Value Of K

As you decrease the value of K to 1 (K = 1), you’ll reach a query point, where you get to see many elements from class A (-) and class B (+) where (-) is the only nearest neighbor. Reasonably, you would think about the query point to be most likely the red one. As K =1, which has a blue color, KNN incorrectly predicts the wrong color blue.
As you increase the value of K to 2 (K=2), you get to see two elements, (-) and (+) are the only nearest neighbor. As you have two values, which are of Class A and Class B, KNN incorrectly predicts the wrong values (Blue and Red).
As you increase the value of K to 3 (K=3), you get to see three elements (-) and (+), (+) are the only nearest neighbor. And this time, you got three values, one from blue and two from red. As your assumption is red, KNN correctly predicts the right value (Blue and Red, Red). Your answer is more stable this time compared to previous ones.

Conclusion

KNN works by finding the nearest distance between a query and all the elements in the database. By choosing the value for K, we get the closest to the query. And then, KNN algorithms look for the most frequent labels in classification and averages of labels in regression.

Illustrative introductions on dimension reduction

October 2, 2020/in Artificial Intelligence, Data Mining, Data Science, Machine Learning, Main Category, Mathematics/by Yasuto Tamura

“What is your image on dimensions?”

….That might be a cheesy question to ask to reader of Data Science Blog, but most people, with no scientific background, would answer “One dimension is a line, and two dimension is a plain, and we live in three-dimensional world.” After that if you ask “How about the fourth dimension?” many people would answer “Time?”

You can find books or writings about dimensions in various field. And you can use the word “dimension” in normal conversations, in many contexts.

*In Japanese, if you say “He likes two dimension.” that means he prefers anime characters to real women, as is often the case with Japanese computer science students.

The meanings of “dimensions” depend on the context, but in data science dimension is usually the number of rows of your Excel data.

When you study data science or machine learning, usually you should start with understanding the algorithms with 2 or 3 dimensional data, and you can apply those ideas to any D dimensional data. But of course you cannot visualize D dimensional data anymore, and you always have to be careful of what happens if you expand degree of dimension.

Conversely it is also important to reduce dimension to understand abstract high dimensional stuff in 2 or 3 dimensional space, which are close to our everyday sense. That means dimension reduction is one powerful way of data visualization.

In this blog series I am going to explain meanings of dimension itself in machine learning context and algorithms for dimension reductions, such as PCA, LDA, and t-SNE, with 2 or 3 dimensional visible data. Along with that, I am going to delve into the meaning of calculations so that you can understand them in more like everyday-life sense.

This article series is going to be roughly divided into the contents below.

Curse of Dimensionality
PCA, LDA (to be published soon)
Rethinking eigen vectors (to be published soon)
KL expansion and subspace method (to be published soon)
Autoencoder as dimension reduction (to be published soon)
t-SNE (to be published soon)

I hope you could see that reducing dimension is one of the fundamental approaches in data science or machine learning.

Simple RNN: the first foothold for understanding LSTM

June 17, 2020/in Artificial Intelligence, Data Science, Deep Learning, Machine Learning, Main Category, Mathematics/by Yasuto Tamura

*In this article “Densely Connected Layers” is written as “DCL,” and “Convolutional Neural Network” as “CNN.”

In the last article, I mentioned “When it comes to the structure of RNN, many study materials try to avoid showing that RNNs are also connections of neurons, as well as DCL or CNN.” Even if you manage to understand DCL and CNN, you can be suddenly left behind once you try to understand RNN because it looks like a different field. In the second section of this article, I am going to provide a some helps for more abstract understandings of DCL/CNN , which you need when you read most other study materials.

My explanation on this simple RNN is based on a chapter in a textbook published by Massachusetts Institute of Technology, which is also recommended in some deep learning courses of Stanford University.

First of all, you should keep it in mind that simple RNN are not useful in many cases, mainly because of vanishing/exploding gradient problem, which I am going to explain in the next article. LSTM is one major type of RNN used for tackling those problems. But without clear understanding forward/back propagation of RNN, I think many people would get stuck when they try to understand how LSTM works, especially during its back propagation stage. If you have tried climbing the mountain of understanding LSTM, but found yourself having to retreat back to the foot, I suggest that you read through this article on simple RNNs. It should help you to gain a solid foothold, and you would be ready for trying to climb the mountain again.

*This article is the second article of “A gentle introduction to the tiresome part of understanding RNN.”

1, A brief review on back propagation of DCL.

Simple RNNs are straightforward applications of DCL, but if you do not even have any ideas on DCL forward/back propagation, you will not be able to understand this article. If you more or less understand how back propagation of DCL works, you can skip this first section.

Deep learning is a part of machine learning. And most importantly, whether it is classical machine learning or deep learning, adjusting parameters is what machine learning is all about. Parameters mean elements of functions except for variants. For example when you get a very simple function $f(x)=a + bx + cx^2 + dx^3$ , then $x$ is a variant, and $a, b, c, d$ are parameters. In case of classical machine learning algorithms, the number of those parameters are very limited because they were originally designed manually. Such functions for classical machine learning is useful for features found by humans, after trial and errors(feature engineering is a field of finding such effective features, manually). You adjust those parameters based on how different the outputs(estimated outcome of classification/regression) are from supervising vectors(the data prepared to show ideal answers).

In the last article I said neural networks are just mappings, whose inputs are vectors, matrices, or sequence data. In case of DCLs, inputs are vectors. Then what’s the number of parameters ? The answer depends on the the number of neurons and layers. In the example of DCL at the right side, the number of the connections of the neurons is the number of parameters(Would you like to try to count them? At least I would say “No.”). Unlike classical machine learning you no longer need to do feature engineering, but instead you need to design networks effective for each task and adjust a lot of parameters.

*I think the hype of AI comes from the fact that neural networks find features automatically. But the reality is difficulty of feature engineering was just replaced by difficulty of designing proper neural networks.

It is easy to imagine that you need an efficient way to adjust those parameters, and the method is called back propagation (or just backprop). As long as it is about DCL backprop, you can find a lot of well-made study materials on that, so I am not going to cover that topic precisely in this article series. Simply putting, during back propagation, in order to adjust parameters of a layer you need errors in the next layer. And in order calculate the errors of the next layer, you need errors in the next next layer.

*You should not think too much about what the “errors” exactly mean. Such “errors” are defined in this context, and you will see why you need them if you actually write down all the mathematical equations behind backprops of DCL.

The red arrows in the figure shows how errors of all the neurons in a layer propagate backward to a neuron in last layer. The figure shows only some sets of such errors propagating backward, but in practice you have to think about all the combinations of such red arrows in the whole back propagation(this link would give you some ideas on how DCLs work).

These points are minimum prerequisites for continuing reading this RNN this article. But if you are planning to understand RNN forward/back propagation at an abstract/mathematical level that you can read academic papers, I highly recommend you to actually write down all the equations of DCL backprop. And if possible you should try to implement backprop of three-layer DCL.

2, Forward propagation of simple RNN

In fact the simple RNN which we are going to look at in this article has only three layers. From now on imagine that inputs of RNN come from the bottom and outputs go up. But RNNs have to keep information of earlier times steps during upcoming several time steps because as I mentioned in the last article RNNs are used for sequence data, the order of whose elements is important. In order to do that, information of the neurons in the middle layer of RNN propagate forward to the middle layer itself. Therefore in one time step of forward propagation of RNN, the input at the time step propagates forward as normal DCL, and the RNN gives out an output at the time step. And information of one neuron in the middle layer propagate forward to the other neurons like yellow arrows in the figure. And the information in the next neuron propagate forward to the other neurons, and this process is repeated. This is called recurrent connections of RNN.

*To be exact we are just looking at a type of recurrent connections. For example Elman RNNs have simpler recurrent connections. And recurrent connections of LSTM are more complicated.

Whether it is a simple one or not, basically RNN repeats this process of getting an input at every time step, giving out an output, and making recurrent connections to the RNN itself. But you need to keep the values of activated neurons at every time step, so virtually you need to consider the same RNNs duplicated for several time steps like the figure below. This is the idea of unfolding RNN. Depending on contexts, the whole unfolded DCLs with recurrent connections is also called an RNN.

In many situations, RNNs are simplified as below. If you have read through this article until this point, I bet you gained some better understanding of RNNs, so you should little by little get used to this more abstract, blackboxed way of showing RNN.

You have seen that you can unfold an RNN, per time step. From now on I am going to show the simple RNN in a simpler way, based on the MIT textbook which I recomment. The figure below shows how RNN propagate forward during two time steps $(t-1), (t)$ .

The input $\boldsymbol{x}^{(t-1)}$ at time step $(t-1)$ propagate forward as a normal DCL, and gives out the output $\hat{\boldsymbol{y}} ^{(t)}$ (The notation on the $\boldsymbol{y} ^{(t)}$ is called “hat,” and it means that the value is an estimated value. Whatever machine learning tasks you work on, the outputs of the functions are just estimations of ideal outcomes. You need to adjust parameters for better estimations. You should always be careful whether it is an actual value or an estimated value in the context of machine learning or statistics). But the most important parts are the middle layers.

*To be exact I should have drawn the middle layers as connections of two layers of neurons like the figure at the right side. But I made my figure closer to the chart in the MIT textbook, and also most other study materials show the combinations of the two neurons before/after activation as one neuron.

$\boldsymbol{a}^{(t)}$ is just linear summations of $\boldsymbol{x}^{(t)}$ (If you do not know what “linear summations” mean, please scroll this page a bit), and $\boldsymbol{h}^{(t)}$ is a combination of activated values of $\boldsymbol{a}^{(t)}$ and linear summations of $\boldsymbol{h}^{(t-1)}$ from the last time step, with recurrent connections. The values of $\boldsymbol{h}^{(t)}$ propagate forward in two ways. One is normal DCL forward propagation to $\hat{\boldsymbol{y}} ^{(t)}$ and $\boldsymbol{o}^{(t)}$ , and the other is recurrent connections to $\boldsymbol{h}^{(t+1)}$ .

These are equations for each step of forward propagation.

$\boldsymbol{a}^{(t)} = \boldsymbol{b} + \boldsymbol{W} \cdot \boldsymbol{h}^{(t-1)} + \boldsymbol{U} \cdot \boldsymbol{x}^{(t)}$
$\boldsymbol{h}^{(t)}= g(\boldsymbol{a}^{(t)})$
$\boldsymbol{o}^{(t)} = \boldsymbol{c} + \boldsymbol{V} \cdot \boldsymbol{h}^{(t)}$
$\hat{\boldsymbol{y}} ^{(t)} = f(\boldsymbol{o}^{(t)})$

*Please forgive me for adding some mathematical equations on this article even though I pledged not to in the first article. You can skip the them, but for some people it is on the contrary more confusing if there are no equations. In case you are allergic to mathematics, I prescribed some treatments below.

*Linear summation is a type of weighted summation of some elements. Concretely, when you have a vector $\boldsymbol{x}=(x_0, x_1, x_2)$ , and weights $\boldsymbol{w}=(w_0,w_1, w_2)$ , then $\boldsymbol{w}^T \cdot \boldsymbol{x} = w_0 \cdot x_0 + w_1 \cdot x_1 +w_2 \cdot x_2$ is a linear summation of $\boldsymbol{x}$ , and its weights are $\boldsymbol{w}$ .

*When you see a product of a matrix and a vector, for example a product of $\boldsymbol{W}$ and $\boldsymbol{v}$ , you should clearly make an image of connections between two layers of a neural network. You can also say each element of $\boldsymbol{u}}$ is a linear summations all the elements of $\boldsymbol{v}}$ , and $\boldsymbol{W}$ gives the weights for the summations.

A very important point is that you share the same parameters, in this case $\boldsymbol{\theta \in \{\boldsymbol{U}, \boldsymbol{W}, \boldsymbol{b}, \boldsymbol{V}, \boldsymbol{c} \}}$ , at every time step.

And you are likely to see this RNN in this blackboxed form.

3, The steps of back propagation of simple RNN

In the last article, I said “I have to say backprop of RNN, especially LSTM (a useful and mainstream type or RNN), is a monster of chain rules.” I did my best to make my PowerPoint on LSTM backprop straightforward. But looking at it again, the LSTM backprop part still looks like an electronic circuit, and it requires some patience from you to understand it. If you want to understand LSTM at a more mathematical level, understanding the flow of simple RNN backprop is indispensable, so I would like you to be patient while understanding this step (and you have to be even more patient while understanding LSTM backprop).

This might be a matter of my literacy, but explanations on RNN backprop are very frustrating for me in the points below.

Most explanations just show how to calculate gradients at each time step.
Most study materials are visually very poor.
Most explanations just emphasize that “errors are back propagating through time,” using tons of arrows, but they lack concrete instructions on how actually you renew parameters with those errors.

If you can relate to the feelings I mentioned above, the instructions from now on could somewhat help you. And I am going to share some study materials on simple RNNs in an external link so that you can gain a clear and mathematical understanding on how simple RNNs work.

Backprop of RNN , as long as you are thinking about simple RNNs, is not so different from that of DCLs. But you have to be careful about the meaning of errors in the context of RNN backprop. Back propagation through time (BPTT) is one of the major methods for RNN backprop, and I am sure most textbooks explain BPTT. But most study materials just emphasize that you need errors from all the time steps, and I think that is very misleading and confusing.

You need all the gradients to adjust parameters, but you do not necessarily need all the errors to calculate those gradients. Gradients in the context of machine learning mean partial derivatives of error functions (in this case $J$ ) with respect to certain parameters, and mathematically a gradient of $J$ with respect to $\boldsymbol{\theta \in \{\boldsymbol{U}, \boldsymbol{W}, \boldsymbol{b}^{(t)}, \boldsymbol{V}, \boldsymbol{c} \}}$ is denoted as $( \frac{\partial J}{\partial \boldsymbol{\theta}} )$ . And another confusing point in many textbooks, including the MIT one, is that they give an impression that parameters depend on time steps. For example some study materials use notations like $\frac{\partial J}{\partial \boldsymbol{\theta}^{(t)}}$ , and I think this gives an impression that this is a gradient with respect to the parameters at time step $(t)$ . In my opinion this gradient rather should be written as $( \frac{\partial J}{\partial \boldsymbol{\theta}} )^{(t)}$ . But many study materials denote gradients of those errors in the former way, so from now on let me use the notations which you can see in the figures in this article.

In order to calculate the gradient $\frac{\partial J}{\partial \boldsymbol{x}^{(t)}}$ you need errors from time steps $s (s \geq t) \quad$ (as you can see in the figure, in order to calculate a gradient in a colored frame, you need all the errors in the same color).

*Another confusing point is that the $\frac{\partial J}{\partial \boldsymbol{\ast ^{(t)}}}, \boldsymbol{\ast} \in \{\boldsymbol{a}^{(t)}, \boldsymbol{h}^{(t)}, \boldsymbol{o}^{(t)}, \dots \}$ are correct notations, because $\boldsymbol{\ast}$ are values of neurons after forward propagation. They depend on time steps, and these are very values which I have been calling “errors.” That is why parameters do not depend on time steps, whereas errors depend on time steps.

As I mentioned before, you share the same parameters at every time step. Again, please do not assume that parameters are different from time step to time step. It is gradients/errors (you need errors to calculate gradients) which depend on time step. And after calculating errors at every time step, you can finally adjust parameters one time, and that’s why this is called “back propagation through time.” (It is easy to imagine that this method can be very inefficient. If the input is the whole text on a Wikipedia link, you need to input all the sentences in the Wikipedia text to renew parameters one time. To solve this problem there is a backprop method named “truncated BPTT,” with which you renew parameters based on a part of a text. )

And after calculating those gradients $\frac{\partial J}{\partial \boldsymbol{\theta}^{(t)}}$ you can take a summation of them: $\frac{\partial J}{\partial \boldsymbol{\theta}}=\sum_{t=0}^{t=\tau}{\frac{\partial J}{\partial \boldsymbol{\theta}^{(t)}}}$ . With this gradient $\frac{\partial J}{\partial \boldsymbol{\theta}}$ , you can finally renew the value of $\boldsymbol{\theta}$ one time.

At the beginning of this article I mentioned that simple RNNs are no longer for practical uses, and that comes from exploding/vanishing problem of RNN. This problem was one of the reasons for the AI winter which lasted for some 20 years. In the next article I am going to write about LSTM, a fancier type of RNN, in the context of a history of neural network history.

Matrix search: Finding the blocks of neighboring fields in a matrix with Python

March 17, 2020/in Mathematics, Python, Tutorial/by Emilia Cheladze

Task

In this article we will look at a solution in python to the following grid search task:

Find the biggest block of adjoining elements of the same kind and into how many blocks the matrix is divided. As adjoining blocks, we will consider field touching by the sides and not the corners.

Input data

For the ease of the explanation, we will be looking at a simple 3×4 matrix with elements of three different kinds, 0, 1 and 2 (see above). To test the code, we will simulate data to achieve different matrix sizes and a varied number of element types. It will also allow testing edge cases like, where all elements are the same or all elements are different.

To simulate some test data for later, we can use the numpy randint() method:

import numpy as np
matrix = [[0,0,1,1], [0,1,2,2], [0,1,1,2]]

matrix_test1 = np.random.randint(3, size = (5,5))
matrix_test2 = np.random.randint(5, size = (10,15))

The code

def find_blocks(matrix):
    visited = []
    block_list = []     
    for x in range(len(matrix)):
        for y in range(len(matrix[0])):
            if (x, y) not in visited:
                field_count, visited = explore_block(x, y, visited)
                block_list.append(field_count)                 
    return print('biggest block: {0}, number of blocks: {1}'
                 .format(max(block_list), len(block_list)))

def explore_block(x, y, visited):
    queue = {(x,y)}
    field_count = 1
    while queue:  
        x,y = queue.pop()
        visited.append((x,y))
        if x+1<len(matrix) and (x+1,y) not in visited and (x+1,y) not in queue:
            if matrix[x+1][y] == matrix[x][y]:
                field_count += 1
                queue.add((x+1,y))         
        if x-1>=0 and (x-1,y) not in visited and (x-1,y) not in queue:
            if matrix[x-1][y] == matrix[x][y]:
                field_count += 1
                queue.add((x-1,y))                    
        if y-1>=0 and (x,y-1) not in visited and (x,y-1) not in queue:
            if matrix[x][y-1] == matrix[x][y]:
                field_count += 1
                queue.add((x,y-1))                    
        if y+1<len(matrix[0]) and (x,y+1) not in visited and (x,y+1) not in queue:
            if matrix[x][y+1] == matrix[x][y]:
                field_count += 1
                queue.add((x,y+1))                                              
    return field_count, visited

How the code works

In summary, the algorithm loops through all fields of the matrix looking for unseen fields that will serve as a starting point for a local exploration of each block of color – the find_blocks() function. The local exploration is done by looking at the neighboring fields and if they are within the same kind, moving to them to explore further fields – the explore_block() function. The fields that have already been seen and counted are stored in the visited list.

find_blocks() function:

Finds a starting point of a new block
Runs a the explore_block() function for local exploration of the block
Appends the size of the explored block
Updates the list of visited points
Returns the result, once all fields of the matrix have been visited.

explore_block() function:

Takes the coordinates of the starting field for a new block and the list of visited points
Creates the queue set with the starting point
Sets the size of the current block (field_count) to 1
Starts a while loop that is executed for as long as the queue is not empty
1. Takes an element of the queue and uses its coordinates as the current location for further exploration
2. Adds the current field to the visited list
3. Explores the neighboring fields and if they belong to the same block, they are added to the queue
4. The fields are taken off the queue for further exploration one by one until the queue is empty
Returns the field_count of the explored block and the updated list of visited fields

Execute the function

find_blocks(matrix)

The returned result is biggest block: 4, number of blocks: 4.

Run the test matrices:

find_blocks(matrix_test1)
find_blocks(matrix_test2)

Visualization

The matrices for the article were visualized with the seaborn heatmap() method.

import seaborn as sns
import matplotlib.pyplot as plt
# use annot=False for a matrix without the number labels
sns.heatmap(matrix, annot=True, center = 0, linewidths=.5, cmap = "viridis",
            xticklabels=False, yticklabels=False, cbar=False, square=True)
plt.show()

Applying Data Science Techniques in Python to Evaluate Ionospheric Perturbations from Earthquakes

February 14, 2018/in Data Science, Data Science Hack, Data Science News, Insights, Main Category, Mathematics, Use Cases, Visualization/by Derry Holding

Multi-GNSS (Galileo, GPS, and GLONASS) Vertical Total Electron Content Estimates: Applying Data Science techniques in Python to Evaluate Ionospheric Perturbations from Earthquakes

1 Introduction

Today, Global Navigation Satellite System (GNSS) observations are routinely used to study the physical processes that occur within the Earth’s upper atmosphere. Due to the experienced satellite signal propagation effects the total electron content (TEC) in the ionosphere can be estimated and the derived Global Ionosphere Maps (GIMs) provide an important contribution to monitoring space weather. While large TEC variations are mainly associated with solar activity, small ionospheric perturbations can also be induced by physical processes such as acoustic, gravity and Rayleigh waves, often generated by large earthquakes.

In this study Ionospheric perturbations caused by four earthquake events have been observed and are subsequently used as case studies in order to validate an in-house software developed using the Python programming language. The Python libraries primarily utlised are Pandas, Scikit-Learn, Matplotlib, SciPy, NumPy, Basemap, and ObsPy. A combination of Machine Learning and Data Analysis techniques have been applied. This in-house software can parse both receiver independent exchange format (RINEX) versions 2 and 3 raw data, with particular emphasis on multi-GNSS observables from GPS, GLONASS and Galileo. BDS (BeiDou) compatibility is to be added in the near future.

Several case studies focus on four recent earthquakes measuring above a moment magnitude (MW) of 7.0 and include: the 11 March 2011 MW 9.1 Tohoku, Japan, earthquake that also generated a tsunami; the 17 November 2013 MW 7.8 South Scotia Ridge Transform (SSRT), Scotia Sea earthquake; the 19 August 2016 MW 7.4 North Scotia Ridge Transform (NSRT) earthquake; and the 13 November 2016 MW 7.8 Kaikoura, New Zealand, earthquake.

Ionospheric disturbances generated by all four earthquakes have been observed by looking at the estimated vertical TEC (VTEC) and residual VTEC values. The results generated from these case studies are similar to those of published studies and validate the integrity of the in-house software.

2 Data Cleaning and Data Processing Methodology

Determining the absolute VTEC values are useful in order to understand the background ionospheric conditions when looking at the TEC perturbations, however small-scale variations in electron density are of primary interest. Quality checking processed GNSS data, applying carrier phase leveling to the measurements, and comparing the TEC perturbations with a polynomial fit creating residual plots are discussed in this section.

Time delay and phase advance observables can be measured from dual-frequency GNSS receivers to produce TEC data. Using data retrieved from the Center of Orbit Determination in Europe (CODE) site (ftp://ftp.unibe.ch/aiub/CODE), the differential code biases are subtracted from the ionospheric observables.

2.1 Determining VTEC: Thin Shell Mapping Function

The ionospheric shell height, H, used in ionosphere modeling has been open to debate for many years and typically ranges from 300 – 400 km, which corresponds to the maximum electron density within the ionosphere. The mapping function compensates for the increased path length traversed by the signal within the ionosphere. Figure 1 demonstrates the impact of varying the IPP height on the TEC values.

Figure 1 Impact on TEC values from varying IPP heights. The height of the thin shell, H, is increased in 50km increments from 300 to 500 km.

2.2 Phase Smoothing

For dual-frequency GNSS users TEC values can be retrieved with the use of dual-frequency measurements by applying calculations. Calculation of TEC for pseudorange measurements in practice produces a noisy outcome and so the relative phase delay between two carrier frequencies – which produces a more precise representation of TEC fluctuations – is preferred. To circumvent the effect of pseudorange noise on TEC data, GNSS pseudorange measurements can be smoothed by carrier phase measurements, with the use of the carrier phase smoothing technique, which is often referred to as carrier phase leveling.

Figure 2 Phase smoothed code differential delay

2.3 Residual Determination

For the purpose of this study the monitoring of small-scale variations in ionospheric electron density from the ionospheric observables are of particular interest. Longer period variations can be associated with diurnal alterations, and changes in the receiver- satellite elevation angles. In order to remove these longer period variations in the TEC time series as well as to monitor more closely the small-scale variations in ionospheric electron density, a higher-order polynomial is fitted to the TEC time series. This higher-order polynomial fit is then subtracted from the observed TEC values resulting in the residuals. The variation of TEC due to the TID perturbation are thus represented by the residuals. For this report the polynomial order applied was typically greater than 4, and was chosen to emulate the nature of the arc for that particular time series. The order number selected is dependent on the nature of arcs displayed upon calculating the VTEC values after an initial inspection of the VTEC plots.

3 Results

3.1 Tohoku Earthquake

For this particular report, the sampled data focused on what was retrieved from the IGS station, MIZU, located at Mizusawa, Japan. The MIZU site is 39N 08′ 06.61″ and 141E 07′ 58.18″. The location of the data collection site, MIZU, and the earthquake epicenter can be seen in Figure 3.

Figure 3 MIZU IGS station and Tohoku earthquake epicenter [generated using the Python library, Basemap]

Figure 4 displays the ionospheric delay in terms of vertical TEC (VTEC), in units of TECU (1 TECU = 1016 el m-2). The plot is split into two smaller subplots, the upper section displaying the ionospheric delay (VTEC) in units of TECU, the lower displaying the residuals. The vertical grey-dashed lined corresponds to the epoch of the earthquake at 05:46:23 UT (2:46:23 PM local time) on March 11 2011. In the upper section of the plot, the blue line corresponds to the absolute VTEC value calculated from the observations, in this case L1 and L2 on GPS, whereby the carrier phase leveling technique was applied to the data set. The VTEC values are mapped from the STEC values which are calculated from the LOS between MIZU and the GPS satellite PRN18 (on Figure 4 denoted G18). For this particular data set as seen in Figure 4, a polynomial fit of five degrees was applied, which corresponds to the red-dashed line. As an alternative to polynomial fitting, band-pass filtering can be employed when TEC perturbations are desired. However for the scope of this report polynomial fitting to the time series of TEC data was the only method used. In the lower section of Figure 4 the residuals are plotted. The residuals are simply the phase smoothed delay values (the blue line) minus the polynomial fit line (the red-dashed line). All ionosphere delay plots follow the same layout pattern and all time data is represented in UT (UT = GPS – 15 leap seconds, whereby 15 leap seconds correspond to the amount of leap seconds at the time of the seismic event). The time series shown for the ionosphere delay plots are given in terms of decimal of the hour, so that the format follows hh.hh.

Figure 4 VTEC and residual plot for G18 at MIZU on March 11 2011

3.2 South Georgia Earthquake

In the South Georgia Island region located in the North Scotia Ridge Transform (NSRT) plate boundary between the South American and Scotia plates on 19 August 2016, a magnitude of 7.4 MW earthquake struck at 7:32:22 UT. This subsection analyses the data retrieved from KEPA and KRSA. As well as computing the GPS and GLONASS TEC values, four Galileo satellites (E08, E14, E26, E28) are also analysed. Figure 5 demonstrates the TEC perturbations as computed for the Galileo L1 and L5 carrier frequencies.

Figure 5 VTEC and residual plots at KRSA on 19 August 2016. The plots are from the perspective of the GNSS receiver at KRSA, for four Galileo satellites (a) E08; (b) E14; (c) E24; (d) E26. The y-axes and x-axes in all plots do not conform with one another but are adjusted to fit the data. The y-axes for the residual section of each plot is consistent with one another.

Figure 6 Geometry of the Galileo (E08, E14, E24 and E26) satellites’ projected ground track whereby the IPP is set to 300km altitude. The orange lines correspond to tectonic plate boundaries.

4 Conclusion

The proximity of the MIZU site and magnitude of the Tohoku event has provided a remarkable – albeit a poignant – opportunity to analyse the ocean-ionospheric coupling aftermath of a deep submarine seismic event. The Tohoku event has also enabled the observation of the origin and nature of the TIDs generated by both a major earthquake and tsunami in close proximity to the epicenter. Further, the Python software developed is more than capable of providing this functionality, by drawing on its mathematical packages, such as NumPy, Pandas, SciPy, and Matplotlib, as well as employing the cartographic toolkit provided from the Basemap package, and finally by utilizing the focal mechanism generation library, Obspy.

Pre-seismic cursors have been investigated in the past and strongly advocated in particular by Kosuke Heki. The topic of pre-seismic ionospheric disturbances remains somewhat controversial. A potential future study area could be the utilization of the Python program – along with algorithmic amendments – to verify the existence of this phenomenon. Such work would heavily involve the use of Scikit-Learn in order to ascertain the existence of any pre-cursors.

Finally, the code developed is still retained privately and as of yet not launched to any particular platform, such as GitHub. More detailed information on this report can be obtained here:

Download as PDF