Optimizing matrix multiplication algorithm

Question

I have implemented my own version of algorithm to multiply two matrices together and I need someone who know what they are doing to see if there is anything that can be done in a more optimal way. I am also keen on understanding how I can make it robust such that it doesn't crash when non-rectangular matrices are passed as argument i.e. matrix that has uneven number of elements in different Vds it holds.

I am particularly interested in function matrixDot in the code bellow, everything else is only to show how I am using it in my project).

#include "iostream"
#include <vector>
#define LOG(m) std::cout << m << std::endl
struct Vd
{
 std::vector<double> v;
};
struct Md
{
 std::vector<Vd> m;
 //fill matrix with num
 void fill(unsigned const int rows, unsigned const int cols, const double num)
 {
 m.clear();
 for (unsigned int i = 0; i < rows; i++)
 {
 Vd tempVec;
 for (unsigned int j = 0; j < cols; j++)
 {
 tempVec.v.push_back(num);
 }
 m.push_back(tempVec);
 }
 }
 friend std::ostream& operator << (std::ostream& out, const Md& mat)
 {
 out << "[" << std::endl << std::endl;
 for (unsigned int i = 0; i < mat.m.size(); i++)
 {
 out << "[";
 for (unsigned int j = 0; j < mat.m[i].v.size(); j++)
 {
 if (j % mat.m[i].v.size() == mat.m[i].v.size() - 1)
 out << mat.m[i].v[j] << "]" << std::endl << std::endl;
 else
 out << mat.m[i].v[j] << ", ";
 }
 }
 out << "]" << std::endl;
 return out;
 }
};
inline void matrixDot(const Md& m1, const Md& m2, Md& outm)
{
 if (m1.m[0].v.size() && m2.m.size())
 if (m1.m[0].v.size() != m2.m.size())
 {
 LOG("Shape mismatch: " << "matrix1 columns: " << m1.m[0].v.size() << ", " << "matrix2 rows: " << m2.m.size());
 throw std::exception();
 }
 unsigned int m1x = 0; unsigned int m1y = 0; unsigned int m2y = 0; //m2x = m1y
 while (outm.m.size() < m1.m.size())
 {
 Vd tempv;
 while (tempv.v.size() < m2.m[0].v.size())
 {
 double total = 0.0;
 while (m1x < m1.m[0].v.size())
 {
 total += m1.m[m1y].v[m1x] * m2.m[m1x].v[m2y];
 m1x++;
 }
 tempv.v.push_back(total);
 m1x = 0;
 m2y < m2.m[0].v.size() - 1 ? m2y++ : m2y = 0;
 }
 m1y < m1.m.size() - 1 ? m1y++ : m1y = 0;
 outm.m.push_back(tempv);
 }
}
int main()
{
 Md mat1;
 mat1.fill(5, 2, 1.0);
 Md mat2;
 mat2.fill(2, 6, 2.0);
 Md mat3;
 matrixDot(mat1, mat2, mat3);
 std::cout << mat3;
}

• Loki Astari: I heard that some of the hard core matrix libraries don't actually do any arithmetic until the values are needed. That way if there are operations you can combine (addition/subtraction) or operations that can be removed (multiply by identity) these can be done before the hard work is done. Or if the value is never used the arithmatic is never actually done.

Answer 1

I see some things that you might want to use to improve your code.

Use using where appropriate

The code currently contains this:

struct Vd
{
 std::vector<double> v;
};

This is probably better expressed like this instead:

using Vd = std::vector<double>;

So now, instead of writing this:

out << mat.m[i].v[j] << ", ";

We can use this cleaner syntax:

out << mat.m[i][j] << ", ";

Understand header include paths

There is a subtle difference between #include "iostream" and #include <iostream>. Although it's implementation defined, most compiler implementations is that the quotation mark form searches locally first (e.g. the current directory) and then searches system include directories if that fails. The angle bracket form usually searches in system include directories. See this question for more. For that reason, this code should probably use #include <iostream>.

Don't use std::endl if you don't really need it

The difference betweeen std::endl and '\n' is that '\n' just emits a newline character, while std::endl actually flushes the stream. This can be time-consuming in a program with a lot of I/O and is rarely actually needed. It's best to only use std::endl when you have some good reason to flush the stream and it's not very often needed for simple programs such as this one. Avoiding the habit of using std::endl when '\n' will do will pay dividends in the future as you write more complex programs with more I/O and where performance needs to be maximized.

Use standard functions where appropriate

Rather than writing the ostream& operator<< as you have, a neat alternative would be to use std::copy and std::experimental::ostream_joiner like this:

friend std::ostream& operator << (std::ostream& out, const Md& mat)
{
 out << "[\n";
 for (const auto& row : mat.m) {
 out << "\t[";
 std::copy(row.begin(), row.end(), std::experimental::make_ostream_joiner(out, ", "));
 out << "]\n";
 }
 return out << "]\n";
}

Prefer return values to output parameters

It seems much more logical to have matrixDot return a new matrix, rather than using the third parameter as an output parameter. See F.20 for more details.

Consider an alternative representation

This code is somewhat fragile in that both the Md and Vd are implemented as struct, with all members public. Worse, is that we could have a jagged array in which each row does not have the same number of items. This is probably not going to result in anything nice. For both of those reasons, I'd suggest using a class and using a single vector to hold all items. See this question for some ideas and advice on how to do that. You might also look at std::valarray as an underlying type.

Provide a full class implementation

In addition to a constructor taking std::initializer_list arguments, I'd also suggest some of the other operators such as operator== for this class.

• Apple_Banana: "Rather than writing the ostream& operator<<" was that a typo and you meant "Rather than writing the ` out <<`" ? or am I missing something here?
• Edward: No, I meant the function itself, the full name of which is std::ostream& operator<<(std::ostream& out, const Md& mat). I was just trying to use a slightly shorter version of that.

Answer 2

I have to admit I'm a bit confused, you did the hard parts and have a question on the easy parts? I might be misunderstanding your question.

I am also keen on understanding how I can make it robust such that it doesn't crash when non-rectangular matrices are passed as argument i.e. matrix that has uneven number of elements in different Vds it holds.

You could validate that you have a well formed matrix by

inline bool isValid(const Md& mat)
{
 if (mat.m.size())
 {
 int size = mat.m[0].v.size();
 for (int i = 1; i < mat.m.size(); i++) {
 if (size != mat.m[i].v.size())
 {
 return false;
 }
 } 
 }
 
 return true;
}

and incorporating it into the matrixDot function for validation similar to the shape validation you have right now

 if (m1.m[0].v.size() && m2.m.size())
 if (m1.m[0].v.size() != m2.m.size())
 {
 LOG("Shape mismatch: " << "matrix1 columns: " << m1.m[0].v.size() << ", " << "matrix2 rows: " << m2.m.size());
 throw std::exception();
 }
 
 if (!isValid(m1))
 {
 LOG("Invalid matrix :: " << std::endl);
 std::cout << m1;
 throw std::exception();
 }
 
 if (!isValid(m2))
 {
 LOG("Invalid matrix :: " << std::endl);
 std::cout << m2;
 throw std::exception();
 }

Any optimizations I can think of are around utilizing std::array instead of std::vector as you are aware of row and column lengths at time of creation.

• Apple_Banana: I did think of loops but then I left like it would be not be worth it to sacrifice performance just to check if data is valid when I most likely would be passing the right data anyway. I was think more of a trick or something to validate without a loop.

Answer 3

Personally I'd extend the Md struct (class) so that it encapsulates the matrix fully. It would have:

--> Member variables for:

 The number of rows and columns. 
 Vector to hold the data. (Consider one dimensional array here).

--> Constructor that would allow a matrix of the right size to be created (rows,cols).

 This would allow you to use vector resize and reserve which will give you 
 a more performant implementation, especially if the matrices are re-used.
 So avoid using push_back and instead set values directly.

--> Get functions to get number of rows/columns

--> Get/Set functions to get/set matrix data values.

 Get/Set functions implement bounds checking.

I'd then derive a mathMatrix class which would add the matrix multiplication functionality. This would necessitate replacing most of the direct access to size functions/data items etc with calls from the above which would make it easier to read and hence maintain.

Then as the earlier poster suggested adding isValid or canMultiply functions would help with making the solution more robust.