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"