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;

    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)  {
        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;

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)

    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;
        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
        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){
            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];
    return data;

int main(){
    vector<Matrix> images = read_data(10);
    vector<Matrix> labels;
    for (Matrix image: images){
        //  image.square_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(brain.difference_squared(images, labels));
        //debug(brain.error_normalized(images, labels));


//alias c="g++ neural.cpp -o neural.exe -O3 -fopenmp --std=c++11;time ./neural.exe"
  • 1
    \$\begingroup\$ 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) \$\endgroup\$ Jun 30, 2016 at 4:54

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.