A Neural Network in 10 lines of CUDA C++ Code

Purpose: For education purposes only. The code demonstrates supervised learning task using a very simple neural network.

Reference: inspired by Andrew Trask‘s post.

Here is a follow-up post featuring a little bit more complicated code:

Neural Network in C++ (Part 2: MNIST Handwritten Digits Dataset)

The core component of the code, the learning algorithm, is only 10 lines:

	__global__ void kFit( const float* X, const int X_w, const int X_h, const float* y, const int y_w, float* l1, const int l1_w, float* l_1_d, float* pred, float* pred_d, float* W0, float* W1, float* buffer) {
	for (unsigned i = 0; i < 50; ++i) {
	dSigmoid(dDot(X, W0, l1, X_h, X_w, l1_w), l1, X_h, l1_w);
	dSigmoid(dDot(l1, W1, pred, X_h, l1_w, y_w), pred, X_h, y_w);
	dMartixByMatrixElementwise(dMartixSubstractMatrix(y, pred, pred_d, X_h, y_w), dSigmoid_d(pred, buffer, X_h, y_w), pred_d, X_h, y_w );
	dMartixByMatrixElementwise(dDot_m1_m2T(pred_d, W1, l_1_d, X_h, y_w, l1_w), dSigmoid_d(l1, buffer, X_h, l1_w), l_1_d, X_h, l1_w);
	dDot_m1T_m2( l1, pred_d, W1, X_h, l1_w, y_w );
	dDot_m1T_m2( X, l_1_d, W0, X_h, X_w, l1_w );
	}
	}

view raw learn.cu hosted with ❤ by GitHub

The loop above runs for 50 iterations (epochs) and fits the vector of attributes X to the vector of classes y. I am going to use 4 records from Iris flower dataset. The attributes (X) are sepal length, sepal width, petal length, and petal width. In my example, I have 2 (Iris Setosa (0) and Iris Virginica (1)) of 3 classes you can find in the original dataset. Predictions are stored in vector pred.

Neural network architecture. Values of vectors W0, W1, layer_1 and pred change over the course of training the network, while vectors X and y must not be changed:

	X W0 layer_1 W1 pred y
	5.1 3.5 1.4 0.2 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.05 0
	4.9 3.0 1.4 0.2 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.06 0
	6.2 3.4 5.4 2.3 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.93 1
	5.9 3.0 5.1 1.8 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.92 1
	0.5
	0.5
	0.5
	0.5

view raw init.cu hosted with ❤ by GitHub

The size of matrix X is the size of the batch by the number of attributes.

Line 3. Finding the values of the hidden layer:

dSigmoid(dDot(X, W0, l1, X_h, X_w, l1_w), l1, X_h, l1_w);

view raw layer1.cu hosted with ❤ by GitHub

In order to calculate the hidden layer, first of all, we will need to multiply a 4 x 4 matrix X by a 4 x 4 matrix W0. Then, we will need to apply an activation function; in this case, we will use a sigmoid function.

A subroutine for matrix multiplication:

	__global__ void kDot(const float m1, const float m2, float *output, const int m1_rows , const int m1_columns, const int m2_columns ){
	/* Computes the product of two matrices: m1 x m2.
	Inputs:
	m1: array, left matrix of size m1_rows x m1_columns
	m2: array, right matrix of size m1_columns x m2_columns (the number of rows in the right matrix
	must be equal to the number of the columns in the left one)
	output: array, the results of the computation are to be stored here:
	m1 * m2, product of two arrays m1 and m2, a matrix of size m1_rows x m2_columns
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_columns: int, number of columns in the right matrix m2
	*/
	const int id = blockIdx.x * blockDim.x + threadIdx.x;
	const int r = (int)id / m2_columns;
	const int c = id % m2_columns;
	float t_output = 0.f;

	for( int k = 0; k < m1_columns; ++k ) {
	t_output += m1[ r * m1_columns + k ] * m2[ k * m2_columns + c ];
	}

	output[ id ] = t_output;
	}

	__device__ float* dDot(const float m1, const float m2, float *output, const int m1_rows , const int m1_columns, const int m2_columns )
	{
	kDot <<< m1_rows, m2_columns >>> (m1, m2, output, m1_rows , m1_columns, m2_columns );
	cudaDeviceSynchronize();
	return output;
	}

view raw dot.cu hosted with ❤ by GitHub

A subroutine for the sigmoid function:

	__global__ void kSigmoid(float const input, float output) {
	/* Computes the value of the sigmoid function f(x) = 1/(1 + e^-x).
	Inputs:
	input: array
	output: array, the results of the computation are to be stored here
	*/

	const int id = blockIdx.x * blockDim.x + threadIdx.x;

	output[id] = 1.0 / (1.0 + std::exp(-input[id]));
	}

	__device__ void dSigmoid(float const input, float output, const int height, const int width){
	kSigmoid <<< height, width >>> (input, output);
	cudaDeviceSynchronize();
	}

view raw sigmoid.cu hosted with ❤ by GitHub

Sigmoid function (red) and its first derivative (blue graph):

Line 4. Finding the matrix with predictions pred. In order to do so, we will need to multiply a 4 x 8 matrix l1 by a 8 x 1 matrix W1. Then, we will need to apply an activation function:

dSigmoid(dDot(l1, W1, pred, X_h, l1_w, y_w), pred, X_h, y_w);

view raw pred.cu hosted with ❤ by GitHub

Line 5. Determine the vector of prediction errors pred_d. First, subtract pred from y. Then, calculate sigmoid( pred ) and, finally, multiply (elementwise) the result of these two operations.

dMartixByMatrixElementwise(dMartixSubstractMatrix(y, pred, pred_d, X_h, y_w), dSigmoid_d(pred, buffer, X_h, y_w), pred_d, X_h, y_w );

view raw pred_d.cu hosted with ❤ by GitHub

CUDA kernel which finds the difference between two matrices:

	__global__ void kMartixSubstractMatrix(const float m1, const float m2, float *output) {
	/* Computes the (elementwise) difference between two arrays
	Inputs:
	m1: array
	m2: array
	output: array,the results of the computation are to be stored here
	*/

	const int id = blockIdx.x * blockDim.x + threadIdx.x;

	output[id] = m1[id] - m2[id];
	}

	__device__ float* dMartixSubstractMatrix(const float m1, const float m2, float *output, const int width, const int height){

	kMartixSubstractMatrix <<< width, height >>> ( m1, m2, output );
	cudaDeviceSynchronize();
	return output;
	}

view raw MartixSubstractMatrix.cu hosted with ❤ by GitHub

Elemetwise multiplicaton of two vectors:

	__global__ void kMartixByMatrixElementwise(const float m1, const float m2, float *output) {
	/* Computes the product of two arrays (elementwise multiplication).
	Inputs:
	m1: array
	m2: array
	output: array,the results of the multiplication are to be stored here
	*/

	const int id = blockIdx.x * blockDim.x + threadIdx.x;

	output[id] = m1[id] * m2[id];
	}

	__device__ float* dMartixByMatrixElementwise(const float m1, const float m2, float *output, const int width, const int height){

	kMartixByMatrixElementwise <<< width, height >>> ( m1, m2, output );
	cudaDeviceSynchronize();
	return output;
	}

view raw MartixByMatrixElementwise.cu hosted with ❤ by GitHub

Line 6. Back propagate the prediction errors to l_1_d. First, multiply pred_d by transposed W1. Then, calculate sigmoid( l1 ) and, finally, multiply (elementwise) the result of these two operations.

dMartixByMatrixElementwise(dDot_m1_m2T(pred_d, W1, l_1_d, X_h, y_w, l1_w), dSigmoid_d(l1, buffer, X_h, l1_w), l_1_d, X_h, l1_w);

view raw l_1_d.cu hosted with ❤ by GitHub

A subroutine that multiplies matrix by transposed matrix:

	__global__ void kDot_m1_m2T(const float m1, const float m2, float *output, const int m1_columns, const int m2_rows ){
	/* Updates the output matrix with the product of two matrices: m1 and m2 transposed.
	Inputs:
	m1: array, left matrix of size m1_rows x m1_columns
	m2: array, right matrix of size m2_rows x m1_columns (m2 transposed will be of size m1_columns x m2_rows)
	output: array, the results of the computation are to be stored here:
	m1 * m2, product of two arrays m1 and m2, a matrix of size m1_rows x m2_rows
	m1_columns: int, number of columns in the left matrix m1
	m2_rows: int, number of rows in the left matrix m2
	*/

	const int id = blockIdx.x * blockDim.x + threadIdx.x;
	const int r = (int)id / m2_rows;
	const int c = id % m2_rows;
	float t_output = 0.0;
	int id_T;

	for( int k = 0; k < m1_columns; ++k ) {
	id_T = c * m1_columns + k;
	t_output += m1[ r * m1_columns + k ] * m2[ id_T ];
	}

	output[ id ] = t_output;
	}

	__device__ float* dDot_m1_m2T(const float m1, const float m2, float *output, const int m1_rows , const int m1_columns, const int m2_rows )
	{
	kDot_m1_m2T <<< m1_rows, m2_rows >>> ( m1, m2, output, m1_columns, m2_rows );
	cudaDeviceSynchronize();
	return output;

view raw dot_m1_m2T.cu hosted with ❤ by GitHub

Line 7. Update weights W1 with the result of matrix multiplication of transposed l1 and pred_d:

This line computes weight updates. In order to do that, we need to perform matrix multiplication of transposed matrix X by matrix pred_delta.

vector W_delta = dot(transpose( &X[0], 4, 4 ), pred_delta, 4, 4, 1);

view raw w_delta.cpp hosted with ❤ by GitHub

Line 8. Update weights W0 with the result of matrix multiplication of transposed X and l_1_d:

dDot_m1T_m2( X, l_1_d, W0, X_h, X_w, l1_w );

view raw W0.cu hosted with ❤ by GitHub

Complete code:

	//
	// onehiddenlayerperceptron.cu
	// onehiddenlayerperceptron
	//
	// Created by Sergei Bugrov on 8/21/17.
	// Copyright © 2017 Sergei Bugrov. All rights reserved.
	//

	#include <stdio.h>

	__global__ void kMartixByMatrixElementwise(const int nThreads, const float m1, const float m2, float *output) {
	/* Computes the product of two arrays (elementwise multiplication).
	Inputs:
	m1: array
	m2: array
	output: array,the results of the multiplication are to be stored here
	*/
	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	output[i] = m1[i] * m2[i];
	}
	}

	__device__ float* dMartixByMatrixElementwise(const float m1, const float m2, float *output, const int width, const int height){

	kMartixByMatrixElementwise <<< width, height >>> ( width * height, m1, m2, output );
	cudaDeviceSynchronize();
	return output;
	}

	__global__ void kMartixSubstractMatrix(const int nThreads, const float m1, const float m2, float *output) {
	/* Computes the (elementwise) difference between two arrays
	Inputs:
	m1: array
	m2: array
	output: array,the results of the computation are to be stored here
	*/

	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	output[i] = m1[i] - m2[i];
	}
	}

	__device__ float* dMartixSubstractMatrix(const float m1, const float m2, float *output, const int width, const int height){

	kMartixSubstractMatrix <<< width, height >>> ( width * height, m1, m2, output );
	cudaDeviceSynchronize();
	return output;
	}

	__global__ void kSigmoid(const int nThreads, float const input, float output){
	/* Computes the value of the sigmoid function f(x) = 1/(1 + e^-x).
	Inputs:
	input: array
	output: array, the results of the computation are to be stored here
	*/

	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	output[i] = 1.0 / (1.0 + std::exp(-input[i]));
	}
	}

	__device__ void dSigmoid(float const input, float output, const int height, const int width){

	kSigmoid <<< height, width >>> (height * width, input, output);
	cudaDeviceSynchronize();
	}

	__global__ void kSigmoid_d(const int nThreads, float const input, float output) {
	/* Computes the value of the sigmoid function derivative f'(x) = f(x)(1 - f(x)),
	where f(x) is sigmoid function.
	Inputs:
	input: array
	output: array, the results of the computation are to be stored here:
	x(1 - x) for every element of the input matrix m1.
	*/

	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	output[i] = input[i] * (1 - input[i]);
	}
	}

	__device__ float* dSigmoid_d(float const input, float output, const int rows, const int columns){
	kSigmoid_d <<< rows, columns >>> (rows*columns, input, output);
	cudaDeviceSynchronize();
	return output;
	}

	__global__ void kDot(const int nThreads, const float m1, const float m2, float *output, const int m1_rows , const int m1_columns, const int m2_columns ){
	/* Computes the product of two matrices: m1 x m2.
	Inputs:
	m1: array, left matrix of size m1_rows x m1_columns
	m2: array, right matrix of size m1_columns x m2_columns (the number of rows in the right matrix
	must be equal to the number of the columns in the left one)
	output: array, the results of the computation are to be stored here:
	m1 * m2, product of two arrays m1 and m2, a matrix of size m1_rows x m2_columns
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_columns: int, number of columns in the right matrix m2
	*/

	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	int r = (int)i / m2_columns;
	int c = i % m2_columns;
	float t_output = 0.f;

	for( int k = 0; k < m1_columns; ++k ) {
	t_output += m1[ r * m1_columns + k ] * m2[ k * m2_columns + c ];
	}

	output[i] = t_output;
	}
	}

	__device__ float* dDot(const float m1, const float m2, float *output, const int m1_rows , const int m1_columns, const int m2_columns ){

	kDot <<< m1_rows, m2_columns >>> (m1_rows * m2_columns, m1, m2, output, m1_rows , m1_columns, m2_columns );
	cudaDeviceSynchronize();
	return output;
	}

	__global__ void kDot_m1_m2T(const int nThreads, const float m1, const float m2, float *output, const int m1_columns, const int m2_rows ){
	/* Updates the output matrix with the product of two matrices: m1 and m2 transposed.
	Inputs:
	m1: array, left matrix of size m1_rows x m1_columns
	m2: array, right matrix of size m2_rows x m1_columns (m2 transposed will be of size m1_columns x m2_rows)
	output: array, the results of the computation are to be stored here:
	m1 * m2, product of two arrays m1 and m2, a matrix of size m1_rows x m2_rows
	m1_columns: int, number of columns in the left matrix m1
	m2_rows: int, number of rows in the left matrix m2
	*/

	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	int r = (int)i / m2_rows;
	int c = i % m2_rows;
	float t_output = 0.0;
	int id_T;

	for( int k = 0; k < m1_columns; ++k ) {
	id_T = c * m1_columns + k;
	t_output += m1[ r * m1_columns + k ] * m2[ id_T ];
	}

	output[i] = t_output;
	}
	}

	__device__ float* dDot_m1_m2T(const float m1, const float m2, float *output, const int m1_rows , const int m1_columns, const int m2_rows )
	{
	kDot_m1_m2T <<< m1_rows, m2_rows >>> ( m1_rows * m2_rows, m1, m2, output, m1_columns, m2_rows );
	cudaDeviceSynchronize();
	return output;
	}

	__global__ void kDot_m1T_m2(const int nThreads, const float m1, const float m2, float *output, const int m1_rows,
	const int m1_columns, const int m2_columns ){
	/* Increments the output matrix with the product of two matrices: m1 transposed and m2.
	Inputs:
	m1: array, left matrix of size m1_rows x m1_columns (m1 transposed will be of size m1_columns x m1_rows)
	m2: array, right matrix of size m1_rows x m2_columns
	output: array, the results of the computation are to be stored here:
	m1 * m2, product of two arrays m1 and m2, a matrix of size m1_columns x m2_columns
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_rows: int, number of rows in the left matrix m2
	*/

	for (int i = blockIdx.x * blockDim.x + threadIdx.x;
	i < nThreads;
	i += blockDim.x * gridDim.x)
	{
	int r = (int)i / m2_columns;
	int c = i % m2_columns;
	int id_T;
	float t_output = 0.0;

	for( int k = 0; k < m1_rows; ++k ) {
	id_T = k * m1_columns + r;
	t_output += m1[ id_T ] * m2[ k * m2_columns + c ];
	}

	output[i] += t_output;
	}
	}

	__device__ void dDot_m1T_m2(const float m1, const float m2, float *output, const int m1_height , const int m1_width, const int m2_width )
	{
	kDot_m1T_m2 <<< m1_width, m2_width >>> (m1_width * m2_width, m1, m2, output, m1_height, m1_width, m2_width );
	cudaDeviceSynchronize();
	}

	__device__ void kPrintMatrix (const float* M, int h, int w) {
	/* Prints out the input array as h x w matrix.
	Inputs:
	m: vector, matrix of size n_rows x n_columns
	h: int, number of rows in the matrix M
	w: int, number of columns in the matrix M
	*/
	for (int i = 0; i < h; i++){
	for (int j = 0; j < w; j++){
	printf("%f ", M[i*w+j]);
	}
	printf("\n");
	}
	printf("\n");
	}

	__global__ void kFit( const float* X, const int X_w, const int X_h,
	const float* y, const int y_w,
	float* l1, const int l1_w, float* l_1_d,
	float* pred, float* pred_d,
	float* W0,
	float* W1,
	float* buffer
	)
	{
	for (unsigned i = 0; i < 50; ++i) {

	dSigmoid(dDot(X, W0, l1, X_h, X_w, l1_w), l1, X_h, l1_w);
	dSigmoid(dDot(l1, W1, pred, X_h, l1_w, y_w), pred, X_h, y_w);
	dMartixByMatrixElementwise(dMartixSubstractMatrix(y, pred, pred_d, X_h, y_w), dSigmoid_d(pred, buffer, X_h, y_w), pred_d, X_h, y_w );
	dMartixByMatrixElementwise(dDot_m1_m2T(pred_d, W1, l_1_d, X_h, y_w, l1_w), dSigmoid_d(l1, buffer, X_h, l1_w), l_1_d, X_h, l1_w);
	dDot_m1T_m2( l1, pred_d, W1, X_h, l1_w, y_w );
	dDot_m1T_m2( X, l_1_d, W0, X_h, X_w, l1_w );
	}
	}

	int main(void){

	const int TRAINING_SIZE = 4;
	const int TRAINING_DIM = 4;
	const int L1_SIZE = 8;

	// X, the first 4 lines from Iris dataset
	float h_X[TRAINING_SIZE*TRAINING_DIM] = { 5.1, 3.5, 1.4, 0.2,
	4.9, 3.0, 1.4, 0.2,
	6.2, 3.4, 5.4, 2.3,
	5.9, 3.0, 5.1, 1.8 };

	const signed int X_size = sizeof(h_X);

	float *d_X;
	cudaMalloc(&d_X, X_size);
	cudaMemcpy(d_X, h_X, X_size, cudaMemcpyHostToDevice);

	//WEIGHTS_0
	const long signed int W0_size = L1_SIZETRAINING_DIMsizeof(float);
	float h_W0 = (float)malloc(W0_size);
	for (int i = 0; i < L1_SIZE*TRAINING_DIM; i++){
	h_W0[i] = 0.1 * (2.0*rand()/RAND_MAX-1.0);
	}

	float *d_W0;
	cudaMalloc(&d_W0, W0_size);
	cudaMemcpy(d_W0, h_W0, W0_size, cudaMemcpyHostToDevice);

	//LAYER_1, LAYER_1_DELTA AND BUFFER OF LAYER 1 SIZE
	const long signed int L1_size = L1_SIZETRAINING_SIZEsizeof(float);

	float* h_layer_1 = (float*)malloc(L1_size);
	float* h_layer_1_delta = (float*)malloc(L1_size);
	float* h_buffer = (float*)malloc(L1_size);

	for (int i = 0; i < L1_SIZE*TRAINING_SIZE; i++){
	h_layer_1[i] = 0.0;
	h_buffer[i] = 0.0;
	h_layer_1_delta[i] = 0.0;
	}

	float *d_layer_1;
	cudaMalloc(&d_layer_1, L1_size);
	cudaMemcpy(d_layer_1, h_layer_1, L1_size, cudaMemcpyHostToDevice);

	float *d_buffer;
	cudaMalloc(&d_buffer, L1_size);
	cudaMemcpy(d_buffer, h_buffer, L1_size, cudaMemcpyHostToDevice);

	float *d_layer_1_delta;
	cudaMalloc(&d_layer_1_delta, L1_size);
	cudaMemcpy(d_layer_1_delta, h_layer_1_delta, L1_size, cudaMemcpyHostToDevice);

	//WEIGHTS_1
	const long signed int W1_size = L1_SIZE*sizeof(float);
	float h_W1 = (float)malloc(W1_size);
	for (int i = 0; i < L1_SIZE; i++){
	h_W1[i] = 0.1* (2.0*rand()/RAND_MAX-1.0);
	}

	float *d_W1;
	cudaMalloc(&d_W1, W1_size);
	cudaMemcpy(d_W1, h_W1, W1_size, cudaMemcpyHostToDevice);

	//Y
	float h_y[4] = { 0,
	0,
	1,
	1 };
	const signed int y_size = sizeof(h_y);
	float *d_y;
	cudaMalloc(&d_y, y_size);
	cudaMemcpy(d_y, h_y, y_size, cudaMemcpyHostToDevice);

	//PRED AND PRED_DELTA
	float* h_pred = (float*)malloc(y_size);
	float* h_pred_delta = (float*)malloc(y_size);
	for (int i = 0; i < TRAINING_SIZE; i++){
	h_pred[i] = 0.0;
	h_pred_delta[i] = 0.0;
	}

	float *d_pred;
	cudaMalloc(&d_pred, y_size);
	cudaMemcpy(d_pred, h_pred, y_size, cudaMemcpyHostToDevice);

	float *d_pred_delta;
	cudaMalloc(&d_pred_delta, y_size);
	cudaMemcpy(d_pred_delta, h_pred_delta, y_size, cudaMemcpyHostToDevice);

	kFit <<< 1, 1 >>> ( d_X, TRAINING_DIM, TRAINING_SIZE,
	d_y, 1,
	d_layer_1, L1_SIZE, d_layer_1_delta,
	d_pred,
	d_pred_delta,
	d_W0,
	d_W1,
	d_buffer);

	cudaMemcpy(h_pred, d_pred, y_size, cudaMemcpyDeviceToHost);

	cudaFree(d_pred);
	cudaFree(d_X);
	cudaFree(d_y);
	cudaFree(d_layer_1_delta);
	cudaFree(d_pred_delta);
	cudaFree(d_W0);
	cudaFree(d_W1);
	cudaFree(d_buffer);

	free(h_layer_1_delta);
	free(h_pred_delta);
	free(h_W0);
	free(h_W1);
	free(h_buffer);

	for (int i = 0; i < TRAINING_SIZE; i++){
	printf("Prediction[%i] : %f True Value[%i] : %f Error[%i] : %f\n", i, h_pred[i], i, h_y[i], i, h_pred[i] - h_y[i]);
	}

	free(h_pred);
	}

view raw onehiddenlayerperceptron.cu hosted with ❤ by GitHub

Output:

	Prediction[0] : 0.060997 True Value[0] : 0.000000 Error[0] : 0.060997
	Prediction[1] : 0.076193 True Value[1] : 0.000000 Error[1] : 0.076193
	Prediction[2] : 0.927551 True Value[2] : 1.000000 Error[2] : -0.072449
	Prediction[3] : 0.918263 True Value[3] : 1.000000 Error[3] : -0.081737

view raw out.cu hosted with ❤ by GitHub

Compile…

nvcc -arch=sm_50 -rdc=true -lcudadevrt onehiddenlayerperceptron.cu -o perceptron

view raw compile.sh hosted with ❤ by GitHub

… and run

./perceptron

view raw run.sh hosted with ❤ by GitHub

Neural Network in C++ (Part 2: MNIST Handwritten Digits Dataset)

In this post, I’ll describe how a neural network with two hidden layers works. The code is highly unoptimized to make it as simple to understand as possible. I’ll train the model on a part of MNIST dataset. So, you will need to download this file containing both the labels (1st column) and the variables. Size of y is 42000×1, and the size of X is 42000×784. Every line of X is a 28×28 grayscale picture of a handwritten number. Every element of y is a number from 0 to 9.

The whole code is here with the explanations following after it:

	//
	// nn.cpp
	//
	// To compile: g++ -o nn nn.cpp -std=c++11
	// To run: ./nn
	// Created by Sergei Bugrov on 4/20/18.
	// Copyright © 2017 Sergei Bugrov. All rights reserved.
	// Download dataset from: https://drive.google.com/file/d/1OdtwXHf_-2T0aS9HLBnxU3o-72mklCZY/view?usp=sharing

	#include <iostream>
	#include <vector>
	#include <math.h>
	#include <fstream>
	#include <sstream>
	#include <string>
	#include <random>

	using namespace std;

	void print ( const vector <float>& m, int n_rows, int n_columns ) {

	/* "Couts" the input vector as n_rows x n_columns matrix.
	Inputs:
	m: vector, matrix of size n_rows x n_columns
	n_rows: int, number of rows in the left matrix m1
	n_columns: int, number of columns in the left matrix m1
	*/

	for( int i = 0; i != n_rows; ++i ) {
	for( int j = 0; j != n_columns; ++j ) {
	cout << m[ i * n_columns + j ] << " ";
	}
	cout << '\n';
	}
	cout << endl;
	}

	int argmax ( const vector <float>& m ) {

	return distance(m.begin(), max_element(m.begin(), m.end()));
	}

	vector <float> relu(const vector <float>& z){
	int size = z.size();
	vector <float> output;
	for( int i = 0; i < size; ++i ) {
	if (z[i] < 0){
	output.push_back(0.0);
	}
	else output.push_back(z[i]);
	}
	return output;
	}

	vector <float> reluPrime (const vector <float>& z) {
	int size = z.size();
	vector <float> output;
	for( int i = 0; i < size; ++i ) {
	if (z[i] <= 0){
	output.push_back(0.0);
	}
	else output.push_back(1.0);
	}
	return output;
	}

	static vector<float> random_vector(const int size)
	{
	random_device rd;
	mt19937 gen(rd());
	uniform_real_distribution<> distribution(0.0, 0.05);
	static default_random_engine generator;

	vector<float> data(size);
	generate(data.begin(), data.end(), [&]() { return distribution(generator); });
	return data;
	}

	vector <float> softmax (const vector <float>& z, const int dim) {

	const int zsize = static_cast<int>(z.size());
	vector <float> out;

	for (unsigned i = 0; i != zsize; i += dim) {
	vector <float> foo;
	for (unsigned j = 0; j != dim; ++j) {
	foo.push_back(z[i + j]);
	}

	float max_foo = *max_element(foo.begin(), foo.end());

	for (unsigned j = 0; j != dim; ++j) {
	foo[j] = exp(foo[j] - max_foo);
	}

	float sum_of_elems = 0.0;
	for (unsigned j = 0; j != dim; ++j) {
	sum_of_elems = sum_of_elems + foo[j];
	}

	for (unsigned j = 0; j != dim; ++j) {
	out.push_back(foo[j]/sum_of_elems);
	}
	}
	return out;
	}

	vector <float> sigmoid_d (const vector <float>& m1) {

	/* Returns the value of the sigmoid function derivative f'(x) = f(x)(1 - f(x)),
	where f(x) is sigmoid function.
	Input: m1, a vector.
	Output: x(1 - x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = m1[ i ] * (1 - m1[ i ]);
	}

	return output;
	}

	vector <float> sigmoid (const vector <float>& m1) {

	/* Returns the value of the sigmoid function f(x) = 1/(1 + e^-x).
	Input: m1, a vector.
	Output: 1/(1 + e^-x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = 1 / (1 + exp(-m1[ i ]));
	}

	return output;
	}

	vector <float> operator+(const vector <float>& m1, const vector <float>& m2){

	/* Returns the elementwise sum of two vectors.
	Inputs:
	m1: a vector
	m2: a vector
	Output: a vector, sum of the vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> sum (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	sum[i] = m1[i] + m2[i];
	};

	return sum;
	}

	vector <float> operator-(const vector <float>& m1, const vector <float>& m2){

	/* Returns the difference between two vectors.
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 - m2, difference between two vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> difference (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	difference[i] = m1[i] - m2[i];
	};

	return difference;
	}

	vector <float> operator*(const vector <float>& m1, const vector <float>& m2){

	/* Returns the product of two vectors (elementwise multiplication).
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 * m2, product of two vectors m1 and m2
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m1[i] * m2[i];
	};

	return product;
	}

	vector <float> operator*(const float m1, const vector <float>& m2){

	/* Returns the product of a float and a vectors (elementwise multiplication).
	Inputs:
	m1: float
	m2: vector
	Output: vector, m1 * m2, product of two vectors m1 and m2
	*/

	const unsigned long VECTOR_SIZE = m2.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m1 * m2[i];
	};

	return product;
	}

	vector <float> operator/(const vector <float>& m2, const float m1){

	/* Returns the product of a float and a vectors (elementwise multiplication).
	Inputs:
	m1: float
	m2: vector
	Output: vector, m1 * m2, product of two vectors m1 and m2
	*/

	const unsigned long VECTOR_SIZE = m2.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m2[i] / m1;
	};

	return product;
	}

	vector <float> transpose (float *m, const int C, const int R) {

	/* Returns a transpose matrix of input matrix.
	Inputs:
	m: vector, input matrix
	C: int, number of columns in the input matrix
	R: int, number of rows in the input matrix
	Output: vector, transpose matrix mT of input matrix m
	*/

	vector <float> mT (C*R);

	for(unsigned n = 0; n != C*R; n++) {
	unsigned i = n/C;
	unsigned j = n%C;
	mT[n] = m[R*j + i];
	}

	return mT;
	}

	vector <float> dot (const vector <float>& m1, const vector <float>& m2, const int m1_rows, const int m1_columns, const int m2_columns) {

	/* Returns the product of two matrices: m1 x m2.
	Inputs:
	m1: vector, left matrix of size m1_rows x m1_columns
	m2: vector, right matrix of size m1_columns x m2_columns (the number of rows in the right matrix
	must be equal to the number of the columns in the left one)
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_columns: int, number of columns in the right matrix m2
	Output: vector, m1 * m2, product of two vectors m1 and m2, a matrix of size m1_rows x m2_columns
	*/

	vector <float> output (m1_rows*m2_columns);

	for( int row = 0; row != m1_rows; ++row ) {
	for( int col = 0; col != m2_columns; ++col ) {
	output[ row * m2_columns + col ] = 0.f;
	for( int k = 0; k != m1_columns; ++k ) {
	output[ row * m2_columns + col ] += m1[ row * m1_columns + k ] * m2[ k * m2_columns + col ];
	}
	}
	}

	return output;
	}

	vector<string> split(const string &s, char delim) {
	stringstream ss(s);
	string item;
	vector<string> tokens;
	while (getline(ss, item, delim)) {
	tokens.push_back(item);
	}
	return tokens;
	}

	int main(int argc, const char * argv[]) {

	string line;
	vector<string> line_v;

	cout << "Loading data ...\n";
	vector<float> X_train;
	vector<float> y_train;
	ifstream myfile ("train.txt");
	if (myfile.is_open())
	{
	while ( getline (myfile,line) )
	{
	line_v = split(line, '\t');
	int digit = strtof((line_v[0]).c_str(),0);
	for (unsigned i = 0; i < 10; ++i) {
	if (i == digit)
	{
	y_train.push_back(1.);
	}
	else y_train.push_back(0.);
	}

	int size = static_cast<int>(line_v.size());
	for (unsigned i = 1; i < size; ++i) {
	X_train.push_back(strtof((line_v[i]).c_str(),0));
	}
	}
	X_train = X_train/255.0;
	myfile.close();
	}

	else cout << "Unable to open file" << '\n';

	int xsize = static_cast<int>(X_train.size());
	int ysize = static_cast<int>(y_train.size());

	// Some hyperparameters for the NN
	int BATCH_SIZE = 256;
	float lr = .01/BATCH_SIZE;

	// Random initialization of the weights
	vector <float> W1 = random_vector(784*128);
	vector <float> W2 = random_vector(128*64);
	vector <float> W3 = random_vector(64*10);

	cout << "Training the model ...\n";
	for (unsigned i = 0; i < 10000; ++i) {

	// Building batches of input variables (X) and labels (y)
	int randindx = rand() % (42000-BATCH_SIZE);
	vector<float> b_X;
	vector<float> b_y;
	for (unsigned j = randindx784; j < (randindx+BATCH_SIZE)784; ++j){
	b_X.push_back(X_train[j]);
	}
	for (unsigned k = randindx10; k < (randindx+BATCH_SIZE)10; ++k){
	b_y.push_back(y_train[k]);
	}

	// Feed forward
	vector<float> a1 = relu(dot( b_X, W1, BATCH_SIZE, 784, 128 ));
	vector<float> a2 = relu(dot( a1, W2, BATCH_SIZE, 128, 64 ));
	vector<float> yhat = softmax(dot( a2, W3, BATCH_SIZE, 64, 10 ), 10);

	// Back propagation
	vector<float> dyhat = (yhat - b_y);
	// dW3 = a2.T * dyhat
	vector<float> dW3 = dot(transpose( &a2[0], BATCH_SIZE, 64 ), dyhat, 64, BATCH_SIZE, 10);
	// dz2 = dyhat * W3.T * relu'(a2)
	vector<float> dz2 = dot(dyhat, transpose( &W3[0], 64, 10 ), BATCH_SIZE, 10, 64) * reluPrime(a2);
	// dW2 = a1.T * dz2
	vector<float> dW2 = dot(transpose( &a1[0], BATCH_SIZE, 128 ), dz2, 128, BATCH_SIZE, 64);
	// dz1 = dz2 * W2.T * relu'(a1)
	vector<float> dz1 = dot(dz2, transpose( &W2[0], 128, 64 ), BATCH_SIZE, 64, 128) * reluPrime(a1);
	// dW1 = X.T * dz1
	vector<float> dW1 = dot(transpose( &b_X[0], BATCH_SIZE, 784 ), dz1, 784, BATCH_SIZE, 128);

	// Updating the parameters
	W3 = W3 - lr * dW3;
	W2 = W2 - lr * dW2;
	W1 = W1 - lr * dW1;


	if ((i+1) % 100 == 0){
	cout << "-----------------------------------------------Epoch " << i+1 << "--------------------------------------------------" <<"\n";
	cout << "Predictions:" << "\n";
	print ( yhat, 10, 10 );
	cout << "Ground truth:" << "\n";
	print ( b_y, 10, 10 );
	vector<float> loss_m = yhat - b_y;
	float loss = 0.0;
	for (unsigned k = 0; k < BATCH_SIZE*10; ++k){
	loss += loss_m[k]*loss_m[k];
	}
	cout << " Loss " << loss/BATCH_SIZE <<"\n";
	cout << "--------------------------------------------End of Epoch :(------------------------------------------------" <<"\n";
	};
	};

	return 0;
	}

view raw nn.cpp hosted with ❤ by GitHub

The neural network part is pretty short:

	// # Feed forward
	vector<float> a1 = relu(dot( b_X, W1, BATCH_SIZE, 784, 128 ));
	vector<float> a2 = relu(dot( a1, W2, BATCH_SIZE, 128, 64 ));
	vector<float> yhat = softmax(dot( a2, W3, BATCH_SIZE, 64, 10 ), 10);

	// # Back propagation
	vector<float> dyhat = yhat - b_y;
	vector<float> dW3 = dot(transpose( &a2[0], BATCH_SIZE, 64 ), dyhat, 64, BATCH_SIZE, 10);
	vector<float> dz2 = dot(dyhat, transpose( &W3[0], 64, 10 ), BATCH_SIZE, 10, 64) * reluPrime(a2);
	vector<float> dW2 = dot(transpose( &a1[0], BATCH_SIZE, 128 ), dz2, 128, BATCH_SIZE, 64);
	vector<float> dz1 = dot(dz2, transpose( &W2[0], 128, 64 ), BATCH_SIZE, 64, 128) * reluPrime(a1);
	vector<float> dW1 = dot(transpose( &b_X[0], BATCH_SIZE, 784 ), dz1, 784, BATCH_SIZE, 128);

	// # Updating the parameters
	W3 = W3 - lr * dW3;
	W2 = W2 - lr * dW2;
	W1 = W1 - lr * dW1;

view raw ff_bp.cpp hosted with ❤ by GitHub

The most interesting is, probably, backpropagation:
Step 1. Calculate the loss and its derivative

The first thing you need to know here is the loss function you are going to use. Here I use Squared Error (more precisely it is 1/2 SE).

In order to propagate the loss, first of all, we need to calculate the derivative of the loss w.r.t the prediction vector yhat.

vector<float> dyhat = yhat - b_y;

view raw bp_1.cpp hosted with ❤ by GitHub

If you use any othe loss function, you need to find its derivative w.r.t yhat.

Step 2. Calculate the gradient of the matrix of parameters W3

vector<float> dW3 = dot(transpose( &a2[0], BATCH_SIZE, 64 ), dyhat, 64, BATCH_SIZE, 10);

view raw bp_2.cpp hosted with ❤ by GitHub

Step 3. Calculate the gradient of the second hidden layer

	// dz2 = dyhat * W3.T * relu'(a2)
	vector<float> dz2 = dot(dyhat, transpose( &W3[0], 64, 10 ), BATCH_SIZE, 10, 64) * reluPrime(a2);

view raw dz2.cpp hosted with ❤ by GitHub

Step 4. Calculate the gradient of the matrix of parameters W2

	// dW2 = a1.T * dz2
	vector<float> dW2 = dot(transpose( &a1[0], BATCH_SIZE, 128 ), dz2, 128, BATCH_SIZE, 64);

view raw dW2.cpp hosted with ❤ by GitHub

Step 5. Calculate the gradient of the first hidden layer

	// dz1 = dz2 * W2.T * relu'(a1)
	vector<float> dz1 = dot(dz2, transpose( &W2[0], 128, 64 ), BATCH_SIZE, 64, 128) * reluPrime(a1);

view raw dz1.cpp hosted with ❤ by GitHub

Step 6. Calculate the gradient of the matrix of parameters W1

	// dW1 = X.T * dz1
	vector<float> dW1 = dot(transpose( &b_X[0], BATCH_SIZE, 784 ), dz1, 784, BATCH_SIZE, 128);

view raw dW1.cpp hosted with ❤ by GitHub

Step 7. Update parameters W1, W2, and W3

	// Updating the parameters
	W3 = W3 - lr * dW3;
	W2 = W2 - lr * dW2;
	W1 = W1 - lr * dW1;

view raw update_w1w2w3.cpp hosted with ❤ by GitHub

The output will look like the following:

Presenting the Importance of Random Initialization of the Weights

The problem of weights initialization is explained here.

“This turns out to be a mistake, because if every neuron in the network computes the same output, then they will also all compute the same gradients during backpropagation and undergo the exact same parameter updates. In other words, there is no source of asymmetry between neurons if their weights are initialized to be the same.”

Basically, if done improperly, it would result in serious problems with learning features. This post is intended to provide some simple evidence of the importance of the asymmetry in weights initialization.

Configuration of the neural network:

	// XOR Dataset
	vector<float> X {
	0.0, 0.0,
	0.0, 1.0,
	1.0, 0.0,
	1.0, 1.0};

	// Quasi random numbers
	vector<float> W0 {
	-0.07555777, -0.04661271, -0.0982434, 0.01800294,
	-0.03213882, -0.09211904, -0.0674924, 0.04203922};

	vector<float> W1 {
	0.03445321,
	0.07976875,
	-0.0502343,
	-0.0995293};

	// All 0.1s
	vector<float> W0_b {
	0.1, 0.1, 0.1, 0.1,
	0.1, 0.1, 0.1, 0.1};

	vector<float> W1_b {
	0.1,
	0.1,
	0.1,
	0.1};

	vector<float> y {
	0.0,
	1.0,
	1.0,
	0.0 };

view raw XOR_init.cpp hosted with ❤ by GitHub

Learning loop:

	int main(int argc, const char * argv[]) {
	for (unsigned i = 0; i != 10000; ++i) {
	vector<float> layer_1 = sigmoid(dot(X, W0, 4, 2, 4 ) );
	vector<float> layer_2 = sigmoid(dot(layer_1, W1, 4, 4, 1 ) );
	vector<float> layer_2_delta = (y - layer_2) * sigmoid_d(layer_2);
	vector<float> layer_1_delta = dot(layer_2_delta, transpose( &W1[0], 4, 1 ), 4, 1, 4) * sigmoid_d(layer_1);
	W1 = W1 + dot(transpose( &layer_1[0], 4, 4 ), layer_2_delta, 4, 4, 1);
	W0 = W0 + dot(transpose( &X[0], 4, 2 ), layer_1_delta, 2, 4, 4);
	};
	return 0;
	}

view raw main_loop.cpp hosted with ❤ by GitHub

Predictions if initialized with assymetrical weights:

	0.0314169
	0.97234
	0.972556
	0.0295431

view raw layer_2_randon_weights.cpp hosted with ❤ by GitHub

Predictions if all weights are initialized with 0.1s:

	0.500027
	0.500029
	0.500029
	0.50003

view raw layer_2_all_0.1_init.cpp hosted with ❤ by GitHub

After 10,000 iterations the network failed to solve simple XOR problem — embarrasing, kind of.

Complete code:

	//
	// main.cpp
	// mlperceptron
	//
	// Created by Sergei Bugrov on 7/11/17.
	// Copyright © 2017 Sergei Bugrov. All rights reserved.
	//

	#include <iostream>
	#include <vector>
	#include <math.h>

	using std::vector;
	using std::cout;
	using std::endl;

	// XOR Dataset
	vector<float> X {
	0.0, 0.0,
	0.0, 1.0,
	1.0, 0.0,
	1.0, 1.0};

	// Quasi random numbers
	vector<float> W0 {
	-0.07555777, -0.04661271, -0.0982434, 0.01800294,
	-0.03213882, -0.09211904, -0.0674924, 0.04203922};

	vector<float> W1 {
	0.03445321,
	0.07976875,
	-0.0502343,
	-0.0995293};

	// All 0.1s
	vector<float> W0_b {
	0.1, 0.1, 0.1, 0.1,
	0.1, 0.1, 0.1, 0.1};

	vector<float> W1_b {
	0.1,
	0.1,
	0.1,
	0.1};

	vector<float> y {
	0.0,
	1.0,
	1.0,
	0.0 };

	vector <float> sigmoid_d (const vector <float>& m1) {

	/* Returns the value of the sigmoid function derivative f'(x) = f(x)(1 - f(x)),
	where f(x) is sigmoid function.
	Input: m1, a vector.
	Output: x(1 - x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = m1[ i ] * (1 - m1[ i ]);
	}

	return output;
	}

	vector <float> sigmoid (const vector <float>& m1) {

	/* Returns the value of the sigmoid function f(x) = 1/(1 + e^-x).
	Input: m1, a vector.
	Output: 1/(1 + e^-x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	if (fabs(m1[ i ]) < 6)
	{
	output[ i ] = 1 / (1 + exp(-m1[ i ]));
	} else {
	output[ i ] = 0.0;
	}
	}

	return output;
	}

	vector <float> operator+(const vector <float>& m1, const vector <float>& m2){

	/* Returns the elementwise sum of two vectors.
	Inputs:
	m1: a vector
	m2: a vector
	Output: a vector, sum of the vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> sum (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	sum[i] = m1[i] + m2[i];
	};

	return sum;
	}

	vector <float> operator-(const vector <float>& m1, const vector <float>& m2){

	/* Returns the difference between two vectors.
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 - m2, difference between two vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> difference (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	difference[i] = m1[i] - m2[i];
	};

	return difference;
	}

	vector <float> operator*(const vector <float>& m1, const vector <float>& m2){

	/* Returns the product of two vectors (elementwise multiplication).
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 * m2, product of two vectors m1 and m2
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m1[i] * m2[i];
	};

	return product;
	}

	vector <float> operator*(const vector <float>& m1, const float s){

	/* Returns the product of a vector and a number (elementwise multiplication).
	Inputs:
	m1: vector
	s: float
	Output: vector, m1 * s, product of a vector and a number
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m1[i] * s;
	};

	return product;
	}

	vector <float> transpose (float *m, const int C, const int R) {

	/* Returns a transpose matrix of input matrix.
	Inputs:
	m: vector, input matrix
	C: int, number of columns in the input matrix
	R: int, number of rows in the input matrix
	Output: vector, transpose matrix mT of input matrix m
	*/

	vector <float> mT (C*R);

	for(unsigned n = 0; n != C*R; n++) {
	unsigned i = n/C;
	unsigned j = n%C;
	mT[n] = m[R*j + i];
	}

	return mT;
	}

	vector <float> dot (const vector <float>& m1, const vector <float>& m2, const int m1_rows, const int m1_columns, const int m2_columns) {

	/* Returns the product of two matrices: m1 x m2.
	Inputs:
	m1: vector, left matrix of size m1_rows x m1_columns
	m2: vector, right matrix of size m1_columns x m2_columns (the number of rows in the right matrix
	must be equal to the number of the columns in the left one)
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_columns: int, number of columns in the right matrix m2
	Output: vector, m1 * m2, product of two vectors m1 and m2, a matrix of size m1_rows x m2_columns
	*/

	vector <float> output (m1_rows*m2_columns);

	for( int row = 0; row != m1_rows; ++row ) {
	for( int col = 0; col != m2_columns; ++col ) {
	output[ row * m2_columns + col ] = 0.f;
	for( int k = 0; k != m1_columns; ++k ) {
	output[ row * m2_columns + col ] += m1[ row * m1_columns + k ] * m2[ k * m2_columns + col ];
	}
	}
	}

	return output;
	}

	void print ( const vector <float>& m, int n_rows, int n_columns ) {

	/* "Couts" the input vector as n_rows x n_columns matrix.
	Inputs:
	m: vector, matrix of size n_rows x n_columns
	n_rows: int, number of rows in the matrix m
	n_columns: int, number of columns in the matrix m
	*/

	for( int i = 0; i != n_rows; ++i ) {
	for( int j = 0; j != n_columns; ++j ) {
	cout << m[ i * n_columns + j ] << " ";
	}
	cout << '\n';
	}
	cout << endl;
	}

	int main(int argc, const char * argv[]) {

	for (unsigned i = 0; i != 10000; ++i) {

	vector<float> layer_1 = sigmoid(dot(X, W0, 4, 2, 4 ) );

	vector<float> layer_2 = sigmoid(dot(layer_1, W1, 4, 4, 1 ) );

	vector<float> layer_2_delta = (y - layer_2) * sigmoid_d(layer_2);

	vector<float> layer_1_delta = dot(layer_2_delta, transpose( &W1[0], 4, 1 ), 4, 1, 4) * sigmoid_d(layer_1);

	W1 = W1 + dot(transpose( &layer_1[0], 4, 4 ), layer_2_delta, 4, 4, 1);

	W0 = W0 + dot(transpose( &X[0], 4, 2 ), layer_1_delta, 2, 4, 4);

	if (i == 9999){
	print ( layer_2, 4, 1 );
	};
	};

	for (unsigned i = 0; i != 10000; ++i) {

	vector<float> layer_1 = sigmoid(dot(X, W0_b, 4, 2, 4 ) );

	vector<float> layer_2 = sigmoid(dot(layer_1, W1_b, 4, 4, 1 ) );

	vector<float> layer_2_delta = (y - layer_2) * sigmoid_d(layer_2);

	vector<float> layer_1_delta = dot(layer_2_delta, transpose( &W1_b[0], 4, 1 ), 4, 1, 4) * sigmoid_d(layer_1);

	W1_b = W1_b + dot(transpose( &layer_1[0], 4, 4 ), layer_2_delta, 4, 4, 1);

	W0_b = W0_b + dot(transpose( &X[0], 4, 2 ), layer_1_delta, 2, 4, 4);

	if (i == 9999){
	print ( layer_2, 4, 1 );
	};
	};

	return 0;
	}

view raw 1hlperceptron.cpp hosted with ❤ by GitHub

Downloading more than 20 years of The New York Times

Articles for the period from 1987 to present are available without subscription. Their copyright notice is web scraping friendly:

“… you may download material from The New York Times on the Web (one machine readable copy and one print copy per page) for your personal, noncommercial use only.”

Why waste the opportunity to download these articles then?

Please read their terms of service here.
Please subscribe to The New York Times here.

	# -- coding: utf-8 --

	# Sergei Bugrov
	# 7-10-17
	#
	# Downloads all available articles from https://www.nytimes.com
	#
	# usage : python nytimes.py
	#
	# python version : 3.6.1

	import requests, bs4, os, errno, time, datetime, re

	def download_page(url):
	try:
	page = requests.get(url, timeout=10.0)
	except requests.exceptions.Timeout:
	print('Timeout\n')
	return None
	except requests.exceptions.ConnectionError:
	print('ConnectionError\n')
	time.sleep(120)
	return None
	except requests.exceptions.HTTPError:
	print('HTTPError\n')
	return None
	except requests.exceptions.TooManyRedirects:
	print('TooManyRedirects\n')
	return None
	else:
	return page


	def main():

	max_attempts = 10

	r_unwanted = re.compile('[\n\t\r]')

	urls_to_articles = []

	if not os.path.exists('articles/'):
	try:
	os.makedirs('articles/')
	except OSError as e:
	if e.errno != errno.EEXIST:
	raise

	# STEP 1. BUILD THE LIST OF URLS TO ARTICLES
	if not os.path.exists('urls_to_articles.txt'):

	for year in range(1987, datetime.datetime.now().year + 1):

	catalog_page_by_years = 'http://spiderbites.nytimes.com/free_%s/index.html' % (year)

	links_to_parts = []

	attempts = 0

	print('Year: ', year)

	with open('logfile.log', 'w') as f:
	f.write('STEP 1. Year: ' + str(year) + '\n')

	catalog_page = download_page(catalog_page_by_years)

	while not (catalog_page or attempts > max_attempts):
	catalog_page = download_page(catalog_page_by_years)
	attempts += 1

	if catalog_page:
	catalog_page = bs4.BeautifulSoup(catalog_page.text, "lxml")
	if year > 1995:
	links_to_parts.append(['http://spiderbites.nytimes.com%s' % (el.get('href')) for el in catalog_page.select('body > div > div > div > div > div > div > ul > li > a')])
	else:
	links_to_parts.append(['http://spiderbites.nytimes.com/free_%s/%s' % (year, el.get('href')) for el in catalog_page.select('body > div > div > div > div > div > div > ul > li > a')])

	links_to_parts = [item for sublist in links_to_parts for item in sublist]

	for link_to_parts in links_to_parts:

	attempts = 0

	parts_page = download_page(link_to_parts)

	while not (parts_page or attempts > max_attempts):
	parts_page = download_page(link_to_parts)
	attempts += 1

	if parts_page:
	parts_page = bs4.BeautifulSoup(parts_page.text, "lxml")
	urls_to_articles.append([el.get('href') for el in parts_page.select('body > div > div > div > div > ul > li > a')])

	urls_to_articles = [item for sublist in urls_to_articles for item in sublist]

	# Backing up the list of URLs
	with open('urls_to_articles.txt', 'w') as output:
	for u in urls_to_articles:
	output.write('%s\n' % (u.strip()))

	# STEP 2. DOWNLOAD ARTICLES
	# If, at some point, Step 2 is interrupted due to unforeseen
	# circumstances (power outage, loss of internet connection), replace the number
	# (value of the variable url_num) below with the one you will find in the logfile.log
	url_num = 0

	if os.path.exists('urls_to_articles.txt') and len(urls_to_articles) == 0:
	with open('urls_to_articles.txt', 'r') as f:
	urls_to_articles = f.read().splitlines()

	print('Number of articles that are about to be downloaded: ', len(urls_to_articles))

	for url in urls_to_articles[url_num:]:

	if len(url) > 34:

	attempts = 0

	if url_num % 1000 == 0:

	print('Downloading article #', url_num, ' from ', url)

	with open('logfile.log', 'w') as f:
	f.write('STEP 2. Downloading article #' + str(url_num) + ' from ' + url + '\n')

	article_page = download_page(url)

	while not (article_page or attempts > max_attempts):
	article_page = download_page(url)
	attempts += 1

	if article_page:
	article_page = bs4.BeautifulSoup(article_page.text, "lxml")

	title = [el.getText() for el in article_page.find_all(class_="articleHeadline")]
	if len(title) > 0:
	title = title[0]
	else:
	title = [el.getText() for el in article_page.find_all(class_="headline")]

	if len(title) > 0:
	title = title[0]
	else:
	title = ""

	dateline = [el.getText() for el in article_page.find_all(class_="dateline")]
	if len(dateline) > 0:
	dateline = dateline[0]
	else:
	dateline = ""

	byline = [el.getText().strip() for el in article_page.find_all(class_="byline")]
	if len(byline) > 0:
	byline = ' '.join(byline)
	else:
	byline = ""

	body = [el.getText() for el in article_page.find_all(class_="articleBody")]
	if len(body) > 0:
	body = '\n'.join(body)
	body = r_unwanted.sub("", body)
	body = re.sub(' +', ' ', body)

	with open('articles/' + str(url_num) + url.split('/')[-1] + '.txt', 'w') as output:
	output.write('(c) ' + str(datetime.datetime.now().year) + ' The New York Times Company\n')
	output.write(url + '\n')
	output.write(title + '\n')
	output.write(dateline + '\n')
	output.write(byline + '\n')
	output.write('\n' + body)
	else:

	body = [el.getText() for el in article_page.find_all(class_="story-body-text")]

	if len(body) > 0:
	body = '\n'.join(body)
	body = r_unwanted.sub("", body)
	body = re.sub(' +', ' ', body)

	with open('articles/' + str(url_num) + url.split('/')[-1] + '.txt', 'w') as output:
	output.write('(c) ' + str(datetime.datetime.now().year) + ' The New York Times Company\n')
	output.write(url + '\n')
	output.write(title + '\n')
	output.write(dateline + '\n')
	output.write(byline + '\n')
	output.write('\n' + body)
	url_num += 1


	if __name__ == '__main__':
	"""
	The main function is called when nytimes.py is run from the command line
	"""

	main()

view raw nytimes.py hosted with ❤ by GitHub

Next time, I’ll modify the code so you can download articles from some other major online newspaper.

Downloading all English books from gutenberg.org with Python

Project Gutenberg (PG) is probably second most popular source (after Wikipedia: here you will find a torrent file for the latest Wikipedia dump btw) of text corpora for NLP. The code below will download all available books in .txt format in the English language. It consists of two steps: (1) first, it collects all direct URLs to the books and (2) then, it downloads them one by one, extracts text files from archives and, then, deletes .zip files.

After you run the code, you will get approximately 16,486,020,098 bytes (16.57 GB on disk) for 41,599 items.

	# -- coding: utf-8 --

	# Sergei Bugrov
	# 7-9-17
	#
	# Downloads all available books in English language in .txt format from http://www.gutenberg.org,
	# unpacks them from .zip archives, saves them to ../books/ folder, and deletes .zip files.
	#
	# usage : python gutenberg.py
	#
	# python version : 3.6.1

	import requests, bs4, os, errno, zipfile, glob
	from urllib.request import urlretrieve


	def main():

	if not os.path.exists('books/'):
	try:
	os.makedirs('books/')
	except OSError as e:
	if e.errno != errno.EEXIST:
	raise

	# STEP 1. BUILD A LIST OF URLS

	urls_to_books = []

	if not os.path.exists('urls_to_books.txt'):

	page_w_books_url = 'http://www.gutenberg.org/robot/harvest?filetypes%5B%5D=txt&langs%5B%5D=en'

	while 1 == 1:

	is_last_page = False

	print('Reading page: ' + page_w_books_url)

	page_w_books = requests.get(page_w_books_url, timeout=20.0)

	if page_w_books:
	page_w_books = bs4.BeautifulSoup(page_w_books.text, "lxml")
	urls = [el.get('href') for el in page_w_books.select('body > p > a[href^="http://aleph.gutenberg.org/"%5D')]
	url_to_next_page = page_w_books.find_all('a', string='Next Page')

	if len(urls) > 0:
	urls_to_books.append(urls)

	if url_to_next_page[0]:
	page_w_books_url = "http://www.gutenberg.org/robot/" + url_to_next_page[0].get('href')
	else:
	is_last_page = True

	if is_last_page:
	break

	urls_to_books = [item for sublist in urls_to_books for item in sublist]

	# Backing up the list of URLs
	with open('urls_to_books.txt', 'w') as output:
	for u in urls_to_books:
	output.write('%s\n' % u)

	# STEP 2. DOWNLOAD BOOKS

	# If, at some point, Step 2 is interrupted due to unforeseen
	# circumstances (power outage, lost of internet connection), replace the number
	# (value of the variable url_num) below with the one you will find in the logfile.log
	# Example
	# logfile.log : Unzipping file #99 books/10020.zip
	# the number : 99
	url_num = 0

	if os.path.exists('urls_to_books.txt') and len(urls_to_books) == 0:
	with open('urls_to_books.txt', 'r') as f:
	urls_to_books = f.read().splitlines()

	for url in urls_to_books[url_num:]:

	dst = 'books/' + url.split('/')[-1].split('.')[0].split('-')[0]

	with open('logfile.log', 'w') as f:
	f.write('Unzipping file #' + str(url_num) + ' ' + dst + '.zip' + '\n')

	if len(glob.glob(dst + '*')) == 0:
	urlretrieve(url, dst + '.zip')

	with zipfile.ZipFile(dst + '.zip', "r") as zip_ref:
	try:
	zip_ref.extractall("books/")
	print(str(url_num) + ' ' + dst + '.zip ' + 'unzipped successfully!')
	except NotImplementedError:
	print(str(url_num) + ' Cannot unzip file:', dst)

	os.remove(dst + '.zip')

	url_num += 1


	if __name__ == '__main__':
	"""
	The main function is called when gutenberg.py is run from the command line
	"""

	main()

view raw gutenberg.py hosted with ❤ by GitHub

Next time, I will build word embeddings using word2vec model based on the PG text corpus.

A Neural Network in 10 lines of C++ Code

Purpose: For education purposes only. The code demonstrates supervised learning task using a very simple neural network. In my next post, I am going to replace the vast majority of subroutines with CUDA kernels.

Reference: Andrew Trask‘s post.

The core component of the code, the learning algorithm, is only 10 lines:

	int main(int argc, const char * argv[]) {
	for (unsigned i = 0; i != 50; ++i) {
	vector<float> pred = sigmoid(dot(X, W, 4, 4, 1 ) );
	vector<float> pred_error = y - pred;
	vector<float> pred_delta = pred_error * sigmoid_d(pred);
	vector<float> W_delta = dot(transpose( &X[0], 4, 4 ), pred_delta, 4, 4, 1);
	W = W + W_delta;
	};
	return 0;
	}

view raw main_loop.cpp hosted with ❤ by GitHub

The loop above runs for 50 iterations (epochs) and fits the vector of attributes X to the vector of classes y through the vector of weights W. I am going to use 4 records from Iris flower dataset. The attributes (X) are sepal length, sepal width, petal length, and petal width. In my example, I have 2 (Iris Setosa (0) and Iris Virginica (1)) of 3 classes you can find in the original dataset. Predictions are stored in vector pred.

Neural network architecture. Values of vectors W and pred change over the course of training the network, while vectors X and y must not be changed:

	X W pred y
	5.1 3.5 1.4 0.2 0.5 0.00 0
	4.9 3.0 1.4 0.2 0.5 0.00 0
	6.2 3.4 5.4 2.3 0.5 0.99 1
	5.9 3.0 5.1 1.8 0.5 0.99 1

view raw nn_arch.cpp hosted with ❤ by GitHub

The size of matrix X is the size of the batch by the number of attributes.

Line 3. Make predictions:

vector pred = sigmoid(dot(X, W, 4, 4, 1 ) );

view raw pred.cpp hosted with ❤ by GitHub

In order to calculate predictions, first of all, we will need to multiply a 4 x 4 matrix X by a 4 x 1 matrix W. Then, we will need to apply an activation function; in this case, we will use a sigmoid function.

A subroutine for matrix multiplication:

	vector <float> dot (const vector <float>& m1, const vector <float>& m2,
	const int m1_rows, const int m1_columns, const int m2_columns) {

	/* Returns the product of two matrices: m1 x m2.
	Inputs:
	m1: vector, left matrix of size m1_rows x m1_columns
	m2: vector, right matrix of size m1_columns x m2_columns
	(the number of rows in the right matrix must be equal
	to the number of the columns in the left one)
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_columns: int, number of columns in the right matrix m2
	Output: vector, m1 * m2, product of two vectors m1 and m2,
	a matrix of size m1_rows x m2_columns
	*/

	vector <float> output (m1_rows*m2_columns);

	for( int row = 0; row != m1_rows; ++row ) {
	for( int col = 0; col != m2_columns; ++col ) {
	output[ row * m2_columns + col ] = 0.f;
	for( int k = 0; k != m1_columns; ++k ) {
	output[ row * m2_columns + col ] += m1[ row * m1_columns + k ] * m2[ k * m2_columns + col ];
	}
	}
	}

	return output;
	}

view raw dot.cpp hosted with ❤ by GitHub

A subroutine for the sigmoid function:

	vector <float> sigmoid (const vector <float>& m1) {

	/* Returns the value of the sigmoid function f(x) = 1/(1 + e^-x).
	Input: m1, a vector.
	Output: 1/(1 + e^-x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = 1 / (1 + exp(-m1[ i ]));
	}

	return output;
	}

view raw sigmoid.cpp hosted with ❤ by GitHub

Sigmoid function (red) and its first derivative (blue graph):

Line 4. Calculate pred_error, it is simply a difference between the predictions and the truth:

vector<float> pred_error = y - pred;

view raw pred_error.cpp hosted with ❤ by GitHub

In order to subtract one vector from another, we will need to overload the “-” operator:

	vector <float> operator-(const vector <float>& m1, const vector <float>& m2){

	/* Returns the difference between two vectors.
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 - m2, difference between two vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> difference (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	difference[i] = m1[i] - m2[i];
	};

	return difference;
	}

view raw operator-.cpp hosted with ❤ by GitHub

Line 5. Determine the vector of deltas pred_delta:

vector<float> pred_delta = pred_error * sigmoid_d(pred);

view raw pred_delta.cpp hosted with ❤ by GitHub

In order to perform elemetwise multiplicaton of two vectors, we will need to overload the “*” operator:

	vector <float> operator*(const vector <float>& m1, const vector <float>& m2){

	/* Returns the product of two vectors (elementwise multiplication).
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 * m2, product of two vectors m1 and m2
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m1[i] * m2[i];
	};

	return product;
	}

view raw operator*.cpp hosted with ❤ by GitHub

A subroutine for the derivative of the sigmoid function (d_sigmoid):

Basically, we use the first derivative to find the slope of the line tangent to the graph of the sigmoid function. At x = 0 the slope equals to 0.25. The further the prediction is from 0, the closer the slope is to 0: at x = ±10 the slope equals to 0.000045. Hence, the deltas will be small if either the error is small or the network is very confident about its prediction (i.e. abs(x) is greater than 4).

	vector <float> sigmoid_d (const vector <float>& m1) {

	/* Returns the value of the sigmoid function derivative f'(x) = f(x)(1 - f(x)),
	where f(x) is sigmoid function.
	Input: m1, a vector.
	Output: x(1 - x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = m1[ i ] * (1 - m1[ i ]);
	}

	return output;
	}

view raw sigmoid_d.cpp hosted with ❤ by GitHub

Line 6. Calculate W_delta:

This line computes weight updates. In order to do that, we need to perform matrix multiplication of transposed matrix X by matrix pred_delta.

vector W_delta = dot(transpose( &X[0], 4, 4 ), pred_delta, 4, 4, 1);

view raw w_delta.cpp hosted with ❤ by GitHub

The subroutine that transposes matrices:

	vector <float> transpose (float *m, const int C, const int R) {

	/* Returns a transpose matrix of input matrix.
	Inputs:
	m: vector, input matrix
	C: int, number of columns in the input matrix
	R: int, number of rows in the input matrix
	Output: vector, transpose matrix mT of input matrix m
	*/

	vector <float> mT (C*R);

	for(int n = 0; n!=C*R; n++) {
	int i = n/C;
	int j = n%C;
	mT[n] = m[R*j + i];
	}

	return mT;
	}

view raw transpose.cpp hosted with ❤ by GitHub

Line 7. Update the weights W:

W = W + W_delta;

view raw w.cpp hosted with ❤ by GitHub

In order to perform matrix addition operation, we need to overload the “+” operator:

	vector <float> operator+(const vector <float>& m1, const vector <float>& m2){

	/* Returns the elementwise sum of two vectors.
	Inputs:
	m1: a vector
	m2: a vector
	Output: a vector, sum of the vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> sum (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	sum[i] = m1[i] + m2[i];
	};

	return sum;
	}

view raw operator+.cpp hosted with ❤ by GitHub

Complete code:

	//
	// main.cpp
	// mlperceptron
	//
	// Created by Sergei Bugrov on 7/1/17.
	// Copyright © 2017 Sergei Bugrov. All rights reserved.
	//

	#include <iostream>
	#include <vector>
	#include <math.h>

	using std::vector;
	using std::cout;
	using std::endl;

	vector<float> X {
	5.1, 3.5, 1.4, 0.2,
	4.9, 3.0, 1.4, 0.2,
	6.2, 3.4, 5.4, 2.3,
	5.9, 3.0, 5.1, 1.8
	};

	vector<float> y {
	0,
	0,
	1,
	1 };

	vector<float> W {
	0.5,
	0.5,
	0.5,
	0.5};

	vector <float> sigmoid_d (const vector <float>& m1) {

	/* Returns the value of the sigmoid function derivative f'(x) = f(x)(1 - f(x)),
	where f(x) is sigmoid function.
	Input: m1, a vector.
	Output: x(1 - x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = m1[ i ] * (1 - m1[ i ]);
	}

	return output;
	}

	vector <float> sigmoid (const vector <float>& m1) {

	/* Returns the value of the sigmoid function f(x) = 1/(1 + e^-x).
	Input: m1, a vector.
	Output: 1/(1 + e^-x) for every element of the input matrix m1.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> output (VECTOR_SIZE);


	for( unsigned i = 0; i != VECTOR_SIZE; ++i ) {
	output[ i ] = 1 / (1 + exp(-m1[ i ]));
	}

	return output;
	}

	vector <float> operator+(const vector <float>& m1, const vector <float>& m2){

	/* Returns the elementwise sum of two vectors.
	Inputs:
	m1: a vector
	m2: a vector
	Output: a vector, sum of the vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> sum (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	sum[i] = m1[i] + m2[i];
	};

	return sum;
	}

	vector <float> operator-(const vector <float>& m1, const vector <float>& m2){

	/* Returns the difference between two vectors.
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 - m2, difference between two vectors m1 and m2.
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> difference (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	difference[i] = m1[i] - m2[i];
	};

	return difference;
	}

	vector <float> operator*(const vector <float>& m1, const vector <float>& m2){

	/* Returns the product of two vectors (elementwise multiplication).
	Inputs:
	m1: vector
	m2: vector
	Output: vector, m1 * m2, product of two vectors m1 and m2
	*/

	const unsigned long VECTOR_SIZE = m1.size();
	vector <float> product (VECTOR_SIZE);

	for (unsigned i = 0; i != VECTOR_SIZE; ++i){
	product[i] = m1[i] * m2[i];
	};

	return product;
	}

	vector <float> transpose (float *m, const int C, const int R) {

	/* Returns a transpose matrix of input matrix.
	Inputs:
	m: vector, input matrix
	C: int, number of columns in the input matrix
	R: int, number of rows in the input matrix
	Output: vector, transpose matrix mT of input matrix m
	*/

	vector <float> mT (C*R);

	for(unsigned n = 0; n != C*R; n++) {
	unsigned i = n/C;
	unsigned j = n%C;
	mT[n] = m[R*j + i];
	}

	return mT;
	}

	vector <float> dot (const vector <float>& m1, const vector <float>& m2, const int m1_rows, const int m1_columns, const int m2_columns) {

	/* Returns the product of two matrices: m1 x m2.
	Inputs:
	m1: vector, left matrix of size m1_rows x m1_columns
	m2: vector, right matrix of size m1_columns x m2_columns (the number of rows in the right matrix
	must be equal to the number of the columns in the left one)
	m1_rows: int, number of rows in the left matrix m1
	m1_columns: int, number of columns in the left matrix m1
	m2_columns: int, number of columns in the right matrix m2
	Output: vector, m1 * m2, product of two vectors m1 and m2, a matrix of size m1_rows x m2_columns
	*/

	vector <float> output (m1_rows*m2_columns);

	for( int row = 0; row != m1_rows; ++row ) {
	for( int col = 0; col != m2_columns; ++col ) {
	output[ row * m2_columns + col ] = 0.f;
	for( int k = 0; k != m1_columns; ++k ) {
	output[ row * m2_columns + col ] += m1[ row * m1_columns + k ] * m2[ k * m2_columns + col ];
	}
	}
	}

	return output;
	}

	void print ( const vector <float>& m, int n_rows, int n_columns ) {

	/* "Couts" the input vector as n_rows x n_columns matrix.
	Inputs:
	m: vector, matrix of size n_rows x n_columns
	n_rows: int, number of rows in the left matrix m1
	n_columns: int, number of columns in the left matrix m1
	*/

	for( int i = 0; i != n_rows; ++i ) {
	for( int j = 0; j != n_columns; ++j ) {
	cout << m[ i * n_columns + j ] << " ";
	}
	cout << '\n';
	}
	cout << endl;
	}

	int main(int argc, const char * argv[]) {

	for (unsigned i = 0; i != 50; ++i) {

	vector<float> pred = sigmoid(dot(X, W, 4, 4, 1 ) );
	vector<float> pred_error = y - pred;
	vector<float> pred_delta = pred_error * sigmoid_d(pred);
	vector<float> W_delta = dot(transpose( &X[0], 4, 4 ), pred_delta, 4, 4, 1);
	W = W + W_delta;

	if (i == 49){
	print ( pred, 4, 1 );
	};
	};

	return 0;
	}

view raw main.cpp hosted with ❤ by GitHub

Output:

	0.0511965
	0.0696981
	0.931842
	0.899579

	Program ended with exit code: 0

view raw output.cpp hosted with ❤ by GitHub

How to install NVIDIA CUDA 8.0, cuDNN 5.1, TensorFlow, and Keras on Ubuntu 16.04

Please follow the instructions below and you will be rewarded with Keras with Tenserflow backend and, most importantly, GPU support.

The latest version of CUDA Toolkit you can download from here. It is also clear from that page which versions of Ubuntu are supported.

The latest version of cuDNN you can download from here. TensorFlow, however, requires cuDNN 5.1 and GPU card with CUDA Compute Capability 3.0 or higher.

Step 1. Linux

Update apt repositories and install the linux -image-extra-virtual package.
This package includes the kernel module that’s required by the NVIDIA drivers.

sudo apt-get update
sudo apt-get install -y linux-image-extra-virtual

Install the version of the headers that matches the freshly installed kernel from the previous step.

sudo apt-get install linux-source linux-headers-`uname-r`
sudo reboot

Step 2. Python

Download (from here) and Install Anaconda Python 3.6 64 bit

chmod +x Anaconda3-4.4.0-Linux-x86_64.sh
sudo ./Anaconda3-4.3.1-Linux-x86_64.sh

Step 3. NVIDIA Drivers and CUDA

Blacklist Noveau which has a conflict with the NVIDIA Drivers

echo -e "blacklist nouveau\nblacklist lbm-nouveau\noptions nouveau modeset=0\nalias nouveau off\nalias lbm-nouveau off\n" | sudo tee /etc/modprobe.d/blacklist-nouveau.conf

Disable the Kernel Nouveau

echo options nouveau modeset=0 | sudo tee -a /etc/modprobe.d/nouveau-kms.conf
sudo update-initramfs -u
sudo reboot

Download the Installer and make it executable

chmod +x cuda_8.0.61_375.26_linux.run

Hit Ctrl + Alt + F1

Kill X server

sudo systemctl stop lightdm.service
sudo init 3

Run the Installer and accept the license agreement and install samples

sudo sh cuda_8.0.61_375.26_linux.run

Enable NVIDIA Driver

sudo modprobe nvidia

Restart X server

sudo service lightdm restart

Compile and run the deviceQuery sample from the CUDA distribution to validate the NVIDIA driver installation was successful.

cd /home/evg/NVIDIA_CUDA-8.0_Samples/1_Utilities/deviceQuery/
make ./deviceQuery

Step 4. cuDNN v5.1 for CUDA 8.0

Download cuDNN

Unzip the .tar archive

tar -xzf cudnn-8.0-linux-x64-v5.1.tgz

Copy the cuDNN libraries and header file to the CUDA folders

sudo cp cuda/lib64/libcudnn* /usr/local/cuda/lib64
sudo cp cuda/include/cudnn.h /usr/local/cuda/include/
sudo chmod a+r /usr/local/cuda/lib64/libcudnn*

Add some environment variables

gedit ~/.bashrc

Insert the following lines and save the changes

export CUDA_HOME="/usr/local/cuda"
export LD_LIBRARY_PATH="/usr/local/cuda-8.0/lib64"
export PATH="/usr/local/cuda-8.0/bin:$PATH"

Enable changes in bashrc

source ~/.bashrc

Check if the environment variables contain the paths from the previous step

echo $CUDA_HOME
echo $PATH
echo $LD_LIBRARY_PATH

alias sudo='sudo env PATH=$PATH'

Step 5. Tensorflow

Create a conda environment named tensorflow to run a version of Python by invoking the following command:

conda create -n tensorflow

Activate the conda environment by issuing the following command:

source activate tensorflow

Issue a command of the following format to install TensorFlow inside your conda environment:

sudo pip install –ignore-installed –upgrade TF_PYTHON_URL where TF_PYTHON_URL is the URL of the TensorFlow Python package. For example, the following command installs the CPU-only version of TensorFlow for Python 3.6:

sudo pip install --ignore-installed --upgrade https://storage.googleapis.com/tensorflow/linux/gpu/tensorflow_gpu-1.1.0-cp36-cp36m-linux_x86_64.whl

Test Tensorflow

Invoke python from your shell as follows:

python

Enter the following short program inside the python interactive shell:

import tensorflow as tf
hello = tf.constant('Hello, TensorFlow!')
sess = tf.Session()
print(sess.run(hello))

Exit python interactive shell

exit()

Step 6. Keras

sudo pip install keras

Cognitive Demons

Hello World!

Menu

A Neural Network in 10 lines of CUDA C++ Code

Neural Network in C++ (Part 2: MNIST Handwritten Digits Dataset)

Presenting the Importance of Random Initialization of the Weights

Downloading more than 20 years of The New York Times

Downloading all English books from gutenberg.org with Python

A Neural Network in 10 lines of C++ Code

How to install NVIDIA CUDA 8.0, cuDNN 5.1, TensorFlow, and Keras on Ubuntu 16.04