A library to do maths with matrices written from scratch

Question

Background

Some time ago I've encountered some very good articles about neural networks that represented an ANN as a set of matrices, so everything was done using matrix operations. These articles show code written in Python, and since I know Python well, I decided to translate it to C++ (which I'm currently learning). Very soon I understood that it's too difficult to represent a matrix as a vector of vectors and mess with many of these vectors. That's why I decided to write my own library to do maths with matrices. Now it's a project on GitHub called Matrix.

Matrix

I know there are some libraries to work with matrices out there, like BLAS, Eigen, etc. Matrix was originally written for educational purpose, but now I'm also aiming for speed.

You can do almost any mathematical operation with an object of type Matrix as described in the Matrix Wiki.

Internally, all the data is stored in a vector of vectors of double. My attempts to make Matrix a template class were unsuccessful.

---

Code

This is the code responsible for matrix-by-matrix multiplication

Matrix Matrix::operator*(const Matrix& right) const {
    if (cols != right.rows) {
        std::string msg=std::string("Size mismatch while multiplying matrices: ").append(to_string(rows).append(std::string("X")).append(to_string(cols)));
        msg.append(std::string(" vs ").append(to_string(right.rows)).append(std::string("X")).append(to_string(right.cols)));
        throw SizeException(msg);
    }

    if (right.IsNum())
        return this->operator*(right.M[0][0]);


    size_t a, b, c;

    Matrix res(rows, right.cols);


    if (right.IsCol()) {
        for (a = 0; a < cols; ++a)
            res.M[0][0] += M[0][a] * right.M[a][0];
        return res;
    } else if (this->IsSquare(2) && right.IsSquare(2)) {
        // loop unrolling for 2x2 matrices
        res.M[0][0] = M[0][0] * right.M[0][0] + M[0][1] * right.M[1][0],
        res.M[0][1] = M[0][0] * right.M[0][1] + M[0][1] * right.M[1][1],
        res.M[1][0] = M[1][0] * right.M[0][0] + M[1][1] * right.M[1][0],
        res.M[1][1] = M[1][0] * right.M[0][1] + M[1][1] * right.M[1][1];
    
        return res;
    }


    for (a = 0; a < rows; ++a) {
        for (b = 0; b < right.cols; ++b) {
            double tmp;
            for (c = 0, tmp = 0; c < cols; ++c) tmp += M[a][c] * right.M[c][b];
            res.M[a][b] = tmp;
        }
    }

    return res;
}

Here M is a vector containing all the data of a matrix, IsNum and IsCol determine whether a matrix contains only one number and whether it consists of only one column.

---

This code provides operator[] to get a row of a matrix or its single number.

Matrix& Matrix::operator[](const long i) const {

    if (i < 0 || i == rows)
        throw SizeException("Index out of range");

    static Matrix ret;

    if (rows != 1) {
        ret.Reshape(1, cols);
        long a;

        for (a = 0; a < cols; ++a) ret.M[0][a] = this->M[i][a];
    } else {
        ret.Reshape(1, 1);
        ret.M[0][0] = this->M[0][i];
    }

    return ret;
}

Reshape(long rows, long cols) changes the size of a matrix to (rows, columns). I'm wondering whether this method could be made more effective.

---

This is a constructor that accepts two arguments: number of rows and columns in a matrix.

Matrix::Matrix(long rows, long cols) {
    long a;
    this->rows = rows, this->cols = cols;
    this->M.resize(rows);
    for (a = 0; a < rows; ++a) this->M[a].resize(cols);
    this->prettified=false;
}

The problem here is that the profiler shows that a lot of time is spent on resizing the M vector. Is it possible to avoid this?

Questions

1. I'm very new to C++, so my code could look like messed up garbage (although I've put lots of effort to make it look beautiful). What can be done to make my style better?

2. Compared to Python's NumPy, my matrix multiplication is quite a bit slower. Can it be made faster without using sophisticated matrix multiplication algorithms like Strassen's?

3. What else could be done to make Matrix more efficient?

---

Note:

I'm unable to copy and paste the whole code here since it's pretty huge, so please refer to Matrix GitHub repo.

Answer 1

Compared to Python's NumPy, my matrix multiplication is quite a bit slower. Can it be made faster without using sophisticated matrix multiplication algorithms like Strassen's?

One issue is that you use a vector of vectors - your memory accesses aren't going to be as contiguous. Typically one represents a matrix as a single really long vector, and then accesses into it like this matrix[rowNum*numColumns + colNum]. This will greatly improve your memory access patterns (MAPs) and caching behavior. With good enough caching you'll be able to start actually fighting against the memory-boundedness of your program. You'll also have to decide if you want to represent the matrices as row-major or column-major.

Secondly, think about the way matrix multiplication happens - \$row\ x\ column\$. I haven't had a chance to look at your GitHub repo so I don't know if you're representing a Matrix as row-major or column-major. The best way to multiply matrices on a CPU is with a row-major left-matrix and a column-major right-matrix. Obviously you can't always predict when something is going to be used, so you can use a tiled algorithm - essentially you multiply small internal matrices together and then put them all together. This is a more complicated algorithm, but you'll get much better caching behavior for large matrices. It'll look something like this (no guarantees that this is exactly right, I can never remember this off the top of my head). Note that this assumes the tile's size is a factor of the matrices' size, and that the matrices' side lengths are the same, however those can be adjusted for (or you can pad the matrices with 0's).

for tile in result matrix
    tile = 0 // zero out the result matrix, if they aren't initialized to 0
    for tileNum in tilesPerSide
        // These will both have to be adjusted - you'll have to
        // do some computation to actually get these tiles
        leftTile = leftMatrix.getTile(tileNum)
        rightTile = rightMatrix.getTile(tileNum)
        tile += leftTile * rightTile

The real algorithm is quite a bit more complicated than that, and involves an unfortunate number of nested loops, but it'll be faster (on a CPU).

You can also look at things like SIMD vectorization, prefetching, avoiding malloc, and instruction-level parallelism.

SIMD vectorization (or just vectorization) is when you perform a Single Instruction on Multiple Data - basically you perform the same instruction on a couple of (hopefully contiguous in memory) pieces of data. This can improve speed by (usually) a factor of 4 or 8, depending on the size of you computer's vector registers and if you're using floats or doubles. I recommend using Agner Fog's library instead of intrinsics - it's much cleaner, and more OS independent.

Prefetching is just what it sounds like - you ask for data from memory before you actually need it. Cache misses are expensive, so it can be helpful to request memory before you need it. This generally needs some tuning to find a good prefetch distance, and I won't go into all of the specifics here - there are decent tutorials online.

You also want to avoid malloc - this is generally a rather expensive operation. While you may not explicitly malloc or new anything, things like std::vector and std::string use malloc under the hood. This can end up in an operating system call (super slow) or using previously mallocd and freed memory (less slow, but not fantastic). If, for example, you know that messages will be below a certain size, consider using a statically allocated char array instead of a std::string - they can be a little trickier (don't forget the null \0 character) but you won't need to allocate as much on the heap.

Lastly, instruction level parallelism (I've always called this bucketing, but I don't know that anyone else does most people call this loop unrolling) takes advantage of your computer's pipeline and out of order execution to do things faster. If you have multiple subsequent and independent operations, then you can do them like that.

This is how you would normally sum a vector (well, now it isn't, but bear with me):

double result = 0;
for (unsigned i = 0; i < vector.size(); ++i) {
    result += vector[i];
}

Each iteration of the loop is independent of the other, so you can do something like this:

/* 
   Each r# is a "bucket", and the number you can have depends on your machine
   I've never actually been able to figure out how many a given architecture
   supports, however generally 4-8 is pretty safe. If you want to do some
   macro voodoo you can set this with a #define and compile differently 
   depending on architecture
*/
double r1 = 0; double r2 = 0; double r3 = 0; double r4 = 0;
unsigned int i = 0;
/*
   This upper bound is nasty, but it is safer (what if you have a size of 3
   and subtract 4?) and marginally (read as - probably not measurably) 
   faster. You can replace it with vector.size() - 4 if you'd rather.
*/
for (; i < (vector.size() & ~3); i+=4) {
    r1 += vector[i];
    r2 += vector[i+1];
    r3 += vector[i+2];
    r4 += vector[i+3];
}
for (; i < vector.size(); ++i) {
    r1 += vector[i];
}

double result = r1 + r2 + r3 + r4;

This is ugly, but it'll generally be faster due to out of order execution.

Combining a lot of these strategies should give you a substantial speedup.

There are also other parallelization options, GPUs, distributed systems, etc, but I'd argue that isn't really worth worrying about at the moment. If you want to, though, read about OpenMP, OpenACC, Kokkos, CUDA, OpenCL, and MPI and you can take advantage of those things as well.

Answer 2

for loops in C++

In C++ you can declare the loop variable in the statement instead, like this:

for (long a = 0; a < cols; ++a) {
    res.M[0][0] += M[0][a] * right.M[a][0];
}

This is a very nice feature, because the variable a in this example only exists in the scope of the loop, so it cannot be accidentally misused outside of the loop. So if you only need the index variable inside the loop, you should declare in the for statement itself.

Variable scope and live time

Related to the previous remark about the for loop, as a general rule, try to limit variables to the smallest scope possible. This also implies to declare variables when you need them. So instead of this:

size_t a, b, c;

It would be better to declare these variables right before you need them.

Readability

This snippet is hard to read because too many statements are jammed on a single line:

if (cols != right.rows) {
    std::string msg=std::string("Size mismatch while multiplying matrices: ").append(to_string(rows).append(std::string("X")).append(to_string(cols)));
    msg.append(std::string(" vs ").append(to_string(right.rows)).append(std::string("X")).append(to_string(right.cols)));
    throw SizeException(msg);
}

A general recommendation is to have one statement per line. Also, there are several redundant std::string calls there that can be simplified:

if (cols != right.rows) {
    std::string msg = std::string("Size mismatch while multiplying matrices: ")
        .append(to_string(rows))
        .append("X")
        .append(to_string(cols))
    );
    msg.append(" vs ")
        .append(to_string(right.rows))
        .append("X")
        .append(to_string(right.cols))
    );
    throw SizeException(msg);
}

Answer 3

Regarding the first question

What can be done to make my style better?

I've noticed in your Matrix.cpp file the class member variables are referred to inconsistently throughout the different methods.

For example, in Matrix::Zeros, member M is referred to as is, while in Matrix::Ones, you use this->M. There are similar examples in the other methods. Although you seem to correctly use the this pointer to disambiguate the class member variables from the method arguments, many methods are using the this pointer where no such ambiguity exists.

You might also consider making the sqr and exp functions static methods of the Matrix class, that way they don't have to be declared friend. If not, I would at least rename them so they aren't confused with functions from the cmath header.

The Matrix::operator== and Matrix::operator!= have different argument names.

It could be useful to have the name of the file that failed to open when a FileException is thrown.

The std::vector::operator[] takes an unsigned integral type as argument, so you might want to used an unsigned type (size_t) in your for loops when iterating the underlying vector of Matrix class. std::vector::operator[] doesn't check bounds, and if you accidentally pass it a negative number it will likely be converted to a large positive number.

As an aside, if you ever use unsigned types to iterate in reverse (starting variable at array/vector size and decreasing counter), be careful not to subtract one from the counter variable once it reaches zero (the counter variable will likely rollover to the largest possible value).

As far as your Matrix::operator[], I would be careful using that depending on what you intend to do with the returned row. The returned Matrix is a reference to an object static to the method. If you get a row with Matrix::operator[] and then mutate your Matrix, the returned row may no longer match up with your object. You will also mutate any prior Matrix references you retrieved with a call to Matrix::operator[]. So if you ever need multiple rows from a Matrix at the same time, you'll need to change the method to return a new Matrix, or copy the returned row before calling the method again.

Also, as @janos answer points out, the variables only used in the for loops can be declared in the initializer of the loop so they are scoped.