Neural net in C++

Question

I wrote a Matrix struct and a neural net that uses it. Why is this slow? Gprof blames Matrix::operator()(int, int) which I suppose is my constructor, and Matrix::operator*(Matrix) and vector<float>::operator[].

#include <iostream>
//#include <C:/Users/David/Google Drive/Coding/neural_c/Launcher.h>
#include <fstream>
#include <assert.h>
#include <functional>
#include <vector>
#include <math.h>
#include <algorithm>
#include <random>
#include <time.h>
using namespace std;
#define debug(x) cout << #x << " = " << x << endl
#define debug_(x) cout << #x << " = " << x << " : "
#define linearize(r, c) r + rows*c
default_random_engine generator(time(NULL));
normal_distribution<float> _rng(0, 1);
auto rng = [](){return _rng(generator);};
auto logistic = [](float x){return 1/(1+exp(-x));};
auto square = [](float x){return x*x;};
struct Matrix{
 int rows, cols, size;
 int label=-1;
 vector<float> data;
 Matrix(int r=1, int c=1){ //column major
 rows = r;
 cols = c;
 size = r*c;
 data.resize(size);
 }
 Matrix one_to_many(){
 Matrix m(10);
 m.data[label] = 1;
 return m;
 }
 //
 // Matrix(const Matrix &m){
 // rows = m.rows;
 // cols = m.cols;
 // size = m.size;
 // data.resize(size);
 // copy(m.data.begin(), m.data.end(), data.begin());
 // }
 void randomize(){
 for (int i=0; i<size; i++)
 data[i] = rng();
 }
 float &operator()(int r, int c) {
 //assert(0<=r && r<rows && 0<=c && c<cols);
 return data[linearize(r,c)];
 //return data[r+rows*c];
 }
 float &operator[](int i){//linear access
 assert(0<=i && i<size);
 return data[i];
 }
 Matrix operator*(Matrix m) {
 //assert(cols==m.rows);
 Matrix product(rows, m.cols);
 for (int r=0; r<product.rows; r++)
 for (int c=0; c<product.cols; c++)
 for (int i=0; i<cols; i++)
 product(r,c) += (*this)(r, i) * m(i, c);
 return product;
 }
 Matrix operator*(float f){
 Matrix product(rows, cols);
 #pragma omp parallel for
 for (int i=0; i<size; i++)
 product[i] = (*this)[i]*f;
 return product;
 }
 float sum(){
 float sum=0;
 //#pragma omp parallel for reduction(+:sum)
 for (int i=0; i<size; i++)
 sum += (*this)[i];
 return sum;
 }
 float abs_sum(){
 float sum=0;
 //#pragma omp parallel for reduction(+:sum)
 for (int i=0; i<size; i++)
 sum += abs((*this)[i]);
 return sum;
 }
 Matrix operator+(Matrix &m){
 Matrix sum(rows, cols);
 #pragma omp parallel for
 for (int i=0; i<size; i++)
 sum[i] = (*this)[i] + m[i];
 return sum;
 }
 Matrix operator-(Matrix &m){
 //return m;
 Matrix negated=m*-1;
 return (*this) + negated;
 }
 bool operator==(Matrix &m){
 for (int i=0; i<size; i++)
 if (data[i] != m.data[i])
 return false;
 return true;
 }
 template<typename Func>
 void apply(Func f){
 //#pragma omp parallel for
 for (int i=0; i<size; i++)
 (*this)[i] = f((*this)[i]);
 }
 void print(){
 for (int r=0; r<rows; r++){
 for (int c=0; c<cols; c++)
 cout << (*this)(r, c) << " ";
 cout << endl;
 }
 cout << endl;
 }
 void square_print(){
 Matrix square(sqrt(size), sqrt(size));
 square.data = data;
 square.print();
 }
};
struct Net{
 vector<Matrix> weights;
 Net(vector<int> sizes){
 for (int i=0; i<sizes.size()-1; i++){
 weights.push_back(Matrix(sizes[i+1], sizes[i]));
 }
 for (Matrix& weight: weights)
 weight.randomize();
 }
 Matrix operator()(const Matrix &x){
 //return x;
 //for (int i=0; i<10000; i++){x[0]+=.00001;}
 Matrix activation = x;
 for (int i=0; i<weights.size()-1; i++){
 activation = weights[0]*activation;
 activation.apply(logistic);
 }
 return weights[weights.size()-1]*activation;
 }
 float difference_squared(Matrix x, Matrix t){
 Matrix y=(*this)(x); // 5/13
 Matrix difference_squared = t-y; // 4/12
 difference_squared.apply(square);
 return difference_squared.sum();
 }
 float difference_squared(vector<Matrix> X, vector<Matrix> T){
 float total=0;
 for (int i=0; i<X.size(); i++)
 total += difference_squared(X[i], T[i]);
 return total;
 }
 float error_normalized(Matrix x, Matrix t){
 float dif_sq = difference_squared(x, t);
 float t_average_magnitude = t.abs_sum() / t.size;
 return pow(dif_sq, .5) / t_average_magnitude;
 }
 float error_normalized(vector<Matrix> X, vector<Matrix> T){
 float total=0;
 for (int i=0; i<X.size(); i++)
 total += error_normalized(X[i], T[i]);
 return total/X.size();
 }
 float learn(Matrix x, Matrix t, float learning_rate=.0001, float dw=.001){
 float error, error_shifted, partial_derivative;
 for (int layer=0; layer<weights.size(); layer++)
 for (int i=0; i<weights[layer].size; i++){
 error = difference_squared(x, t);
 weights[layer].data[i] += dw;
 error_shifted = difference_squared(x, t);
 weights[layer].data[i] -= dw;
 partial_derivative = (error_shifted - error) / dw;
 weights[layer].data[i] -= learning_rate * partial_derivative;
 }
 return error;
 }
 void print(){
 for (Matrix m: weights){
 m.print();
 cout << "\n";
 }
 }
};
vector<Matrix> read_data(int images=1000){
 ifstream in("mnist_train.csv");
 vector<Matrix> data;
 for (int image=0; image<images; image++){
 Matrix datum(28*28);
 int temp=-1;
 in >> datum.label;
 for(int i=0; i<28*28; i++)
 in >> datum.data[i];
 data.push_back(datum);
 }
 return data;
}
int main(){
 vector<Matrix> images = read_data(10);
 vector<Matrix> labels;
 for (Matrix image: images){
 labels.push_back(image.one_to_many());
 assert(labels[labels.size()-1].abs_sum()==1);
 // image.square_print();
 //labels[labels.size()-1].print();
 }
 printf("data read\n");
 Net brain({28*28, 3, 10});
 printf("brain created!\n");
 //
 int epochs=3;
 debug(brain.difference_squared(images, labels));
 for (int epoch=0; epoch<epochs; epoch++){
 for (int i=0; i<images.size(); i++)
 brain.learn(images[i], labels[i], .01);
 debug_(epoch);
 debug(brain.difference_squared(images, labels));
 //debug(brain.error_normalized(images, labels));
 }
 printf("trained\n");
 printf("done\n");
}
//alias c="g++ neural.cpp -o neural.exe -O3 -fopenmp --std=c++11;time ./neural.exe"

• user1118321: Gprof blames Matrix::operator()(int, int) which I suppose is my constructor - no, it's your access operator. This one: float &operator()(int r, int c)

Answer 1

You are using a naive matrix product, this is slow. There are faster ways to perform a matrix product. See Wikipedia:Matrix Multiplication.

In general you should not write your own math primitives. I would recommend that you use for example Eigen or Armadillo or any of the umpteen linear algebra libraries.