# Super C++ optimization of matrix multiplication with Armadillo

I'm using Armadillo to do very intensive matrix multiplications with side lengths $2^n$, where $n$ can be up to 20 or even more. I'm using Armadillo with OpenBLAS for matrix multiplication, which seems to be doing a very good job in parallel cores, except that I have a problem with the formalism of multiplication in Armadillo for super optimization of performance.

Say that I have a loop in the following form:

arma::cx_mat stateMatrix, evolutionMatrix; //armadillo complex matrix type
for(double t = t0; t < t1; t += 1/sampleRate)
{
...
stateMatrix = evolutionMatrix*stateMatrix;
...
}


In fundamental C++, I find the problem here is that C++ will allocate a new object of cx_mat to store evolutionMatrix*stateMatrix, and then copy the new object to stateMatrix with operator=(). This is very, very inefficient. It's well known that returning complex classes of large data types is a bad idea, right?

The way I see this going way more efficient is with a function doing the multiplication in the form:

void multiply(const cx_mat& mat1, const cx_mat& mat2, cx_mat& output)
{
... //multiplication of mat1 and mat2 and then store it in output
}


This way, One doesn't have to copy huge objects with return value, and the output doesn't have to be reallocated with every multiplication.

The Question: How can I find a compromise, in which I can use Armadillo for multiplication with its nice interface of BLAS, and do this efficiently without having to recreate matrix objects and copy them with each operation?

Isn't this an implementation problem in Armadillo?

• "Super optimization" is actually a thing, which you probably did not mean to refer to. It's a very old and advanced form of compile-time code specialization that still hasn't caught on yet. May 1 '15 at 7:43
• Most of the answers (and the question itself!) seem to miss the point that matrix multiplication is not something you do in place.
– user15247
May 1 '15 at 13:31
• @hurkyl what do you mean with "in place"? May 1 '15 at 14:20
• When you compute $A = A + B$, you modify $A$ "in place" in the sense that you leave the contents of $A$ where they are in memory and do all the work modifying that memory. $A = A * B$ or $A = B * A$ is not computed that way at all. There is no reasonable algorithm for multiplication that leaves $A$ where it is in memory and writes the output of the multiplication into the same memory as its being computed. The update has to be done out of place -- temporary memory has to be used in some fashion.
– user15247
May 1 '15 at 16:34
• Looking at Armadillo's source code, the expression stateMatrix = evolutionMatrix*stateMatrix won't do any copying whatsoever. Instead, Armadillo does a fancy memory pointer change. New memory for the result will still be allocated (there is no way around that), but instead of copying, the stateMatrix matrix will simply use the new memory and discard the old memory. Aug 4 '15 at 16:23

In fundamental C++, I find the problem here is that C++ will allocate a new object of cx_mat to store evolutionMatrix*stateMatrix, and then copy the new object to stateMatrix with operator=().

I think you're right that it's creating temporaries, which is too slow, but I think the reason for why it's doing that is wrong.

Armadillo, like any good C++ linear algebra library, uses expression templates to implement delayed evaluation of expressions. When you write down a matrix product like A*B, no temporaries are created, instead Armadillo makes a temporary object (x) that keeps references to A and B, and then later, given an expression like C = x, computes the matrix product storing the result directly in C, without creating any temporaries.

It also uses this optimization to handle expressions like A*B*C*D, where, depending on matrix sizes, certain orders of multiplication are more efficient than others.

Isn't this an implementation problem in Armadillo?

If Armadillo is not performing this optimization, that would be a bug in Armadillo that should be reported to developers.

However, in your case, there is another problem that is more important. In an expression like A=B*C the storage of A doesn't contain any input data if A doesn't alias B or C. In your case, A = A*B, writing anything to the output matrix would modify one of the input matrices as well.

Even given your suggested function

multiply(const cx_mat&, const cx_mat&, cx_mat&)


how exactly would that function help in the expression multiply(A, B, A)? For most ordinary implementations of that function, that would lead to a bug. It would need to use some temporary storage on its own, to make sure its input data is not corrupted. Your suggestion is pretty much how Armadillo implements matrix multiplication already, but I think it probably takes care to avoid situtations like multiply(A, B, A) by allocating a temporary.

The most likely explanation of why Armadillo isn't doing this optimization is that it would be incorrect to do that.

Finally, there is a much simpler way to do what you want, like this:

cx_mat *A, *Atemp, B;
for (;;) {
*Atemp = (*A)*B;
swap(A, Atemp);
}


This is identical to

cx_mat A, B;
for (;;) {
A = A*B;
}


but it allocates one temporary matrix, instead of one temporary matrix per iteration.

• That “much simpler way to do it” – apart from being obscure-looking (though yeah, swap-instead-of-copy is actually a C++ idiom, fortunately little needed since C++11), and crashing if you don't new-initialise Atemp – doesn't gain you anything at all: it still involves generating a fresh temporary matrix (*A)*B and copying it into *Atemp, unless RVO prevents it. Apr 30 '15 at 23:18
• @leftaroundabout No, if an extra temporary is created in my example, then that's an Armadillo bug. Linear algebra libraries that rely on expression templates explicitly avoid creating temporaries in intermediate results. The value of (*A)*B is not a temporary matrix, but an expression object that keeps track of the expression and its inputs. I tried to explain why this optimisation doesn't fire in the original example, and it's got nothing to do with RVO (or move semantics as in another answer). I skipped all initialisation code, it's not important in the example, I just showed the types. Apr 30 '15 at 23:23
• Ok, I see what you're getting at, but this still seems a very hackish, unreliable way to do it. If the designers had conveived the option to optimise destructive multiplication this way, they'd certainly have implemented it with a dedicated method, or at least provided a custom swap so you don't have to do this kind of pointer juggling. Apr 30 '15 at 23:29
• @leftaroundabout Also, the example doesn't swap matrices, it swaps pointers to matrices, to avoid any copying at all. There are two temporary matrices, and which one of them is considered the temporary one switches every iteration. Apr 30 '15 at 23:29
• @leftaroundabout: There's no memory management going on here with this use of pointers. It's just a small block of code where you have two objects and need to keep track of which object you're using for which purpose.
– user15247
May 1 '15 at 13:15

@BillGreene points to the "return value optimization" as a way around the fundamental problem, but this actually only helps for one half of it. Assume you have code of this form:

struct ExpensiveObject { ExpensiveObject(); ~ExpensiveObject(); };

ExpensiveObject operator+ (ExpensiveObject &obj1,
ExpensiveObject &obj2)
{
ExpensiveObject tmp;
...compute tmp based on obj1 and obj2...
return tmp;
}

void f() {
ExpensiveObject o1, o2, o3;
...initialize o1, o2...;
o3 = o1 + o2;
}


A naive compiler will

1. create a slot to store the result of the plus operation (one temporary),
2. call operator+,
3. create the 'tmp' object inside operator+ (a second temporary),
4. compute tmp,
5. copy tmp into the result slot,
6. destroy tmp,
7. copy the result object into o3
8. destroy the result object

Return value optimization can only unify the 'tmp' object and the 'result' slot, but not remove the need for a copy. So, you are still left with the creation of a temporary, the copy operation, and the destruction of a temporary.

The only way around this is operator+ does not return an object, but an object of some intermediate class that, when assigned to an ExpensiveObject, does the addition and copy operation in place. This is the approach typical used in expression template libraries.

• Thank you for this information. Could you provide an example that I can use with Armadillo to avoid this problem? Apr 30 '15 at 16:09
• And I would like to ask: This is an implementation problem in Armadillo, right? I mean it's not really so smart to do it this way... at least they have to give the result to reference option. Right? Apr 30 '15 at 16:12
• The key part of this answer is the end. Armadillo uses expression templates to evaluate expressions lazily when possible. That cuts down on the number of temporaries that are created. The main thing that the OP should keep in mind is to run a profiler to determine where slowdowns are occurring, then focus on optimizing those. Often, theories about code that "should be slow" don't turn out to be true. Apr 30 '15 at 16:59
• I don't believe any temporaries are created for this example when compiled with a modern C++ compiler. I have created a simple example that shows this and updated my post. I don't disagree with the value of the expression template technique, in general, but it is irrelevant for a simple, single-operator expression like the one shown above. Apr 30 '15 at 18:02
• @BillGreene: Create a class with a constructor, copy constructor, assignment operator and destructor and compile the example. You will see that a temporary is created. Also: needs to be created because the compiler can't eliminate it without merging copy operator, constructor and destructor. That's simply not possible for non-trivial operations such as memory allocation. May 1 '15 at 1:46

Stackoverflow (https://stackoverflow.com/) is probably a better discussion forum for this question. However, here is a short answer.

I doubt that the C++ compiler is generating code for this expression like you describe above. All modern C++ compilers implement an optimization called "return value optimization" (http://en.wikipedia.org/wiki/Return_value_optimization). With return value optimization the result of evolutionMatrix*stateMatrix is stored directly in stateMatrix; no copy is made.

There is obviously considerable confusion on this topic and that is one of the reasons I suggested Stackoverflow might be a better forum. There are many C++ "language lawyers" there while most of us here would rather spend our time on CSE. ;-)

I created the following simple example based on Professor Bangerth's post:

#ifndef NDEBUG
#include <iostream>

using namespace std;
#endif

class ExpensiveObject  {
public:
ExpensiveObject () {
#ifndef NDEBUG
cout << "ExpensiveObject  constructor called." << endl;
#endif
v = 0;
}
ExpensiveObject (int i) {
#ifndef NDEBUG
cout << "ExpensiveObject  constructor(int) called." << endl;
#endif
v = i;
}
ExpensiveObject (const ExpensiveObject  &a) {
v = a.v;
#ifndef NDEBUG
cout << "ExpensiveObject  copy constructor called." << endl;
#endif
}
~ExpensiveObject() {
#ifndef NDEBUG
cout << "ExpensiveObject  destructor called." << endl;
#endif
}
ExpensiveObject  operator=(const ExpensiveObject  &a) {
#ifndef NDEBUG
cout << "ExpensiveObject  assignment operator called." << endl;
#endif
if (this != &a) {
return ExpensiveObject (a);
}
}
void print() const {
#ifndef NDEBUG
cout << "v=" << v << endl;
#endif
}
int getV() const {
return v;
}
private:
int v;
};

ExpensiveObject  operator+(const ExpensiveObject  &a1, const ExpensiveObject  &a2) {
#ifndef NDEBUG
cout << "ExpensiveObject  operator+ called." << endl;
#endif
return ExpensiveObject (a1.getV() + a2.getV());
}

int main()
{
ExpensiveObject  a(2), b(3);
ExpensiveObject  c = a + b;
#ifndef NDEBUG
c.print();
#endif
}


It looks more complicated than it actually is because I wanted to completely remove all code for printing output when compiling in optimized mode. When I run the version compiled with a debug option, I get the following output:

ExpensiveObject  constructor(int) called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  operator+ called.
ExpensiveObject  constructor(int) called.
v=5
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.


The first thing to notice is that no temporaries are constructed-- only a, b, and c. The default constructor and the assignment operator are never called because they aren't needed in this example.

Professor Bangerth mentioned expression templates. Indeed, this optimization technique is very important in obtaining good performance in a matrix class library. But it is important only when the object expressions are more complicated than simply a + b. If, for example, my test was instead:

  ExpensiveObject  a(2), b(3), c(9);
ExpensiveObject  d = a + b + c;


I would get the following output:

ExpensiveObject  constructor(int) called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  operator+ called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  operator+ called.
ExpensiveObject  constructor(int) called.
ExpensiveObject  destructor called.
v=14
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.
ExpensiveObject  destructor called.


This case shows the undesirable construction of a temporary (5 calls to the constructor and two calls of operator +). Proper use of expression templates (a topic well beyond the scope of this forum) would prevent this temporary. (For the highly motivated, a particularly readable discussion of expression templates can be found in chapter 18 of http://www.amazon.com/C-Templates-The-Complete-Guide/dp/0201734842).

Finally, the real "proof" of what the compiler is actually doing comes from examining the assembly code output by the compiler. For the first example, when compiled in optimized mode, this code is astonishly simple. All function calls have been optimized away and the assembly code essentially loads 2 into one register, 3 into a second, and adds them.

• I was actually hesitating to put it here or on stackoverflow... I'm pretty sure if I had put it on stackoverflow, someone would've commented that I should've put it here :-) . Anyway; the return value optimization is good news and I didn't know it before (+1). Thanks for that. Unfortunately I don't know anything in assembly code, so that's not a check that I can do. Apr 30 '15 at 13:33
• If I am not mistaken, even considering return-value optimization the compiler works with three matrices in memory, not two. "Multiply A and B, and put the result in C" is a different function than "multiply A and B, and overwrite B with the result". Apr 30 '15 at 14:56
• Interesting point. I was focusing the the desire of the poster to have a matrix multiply implementation as efficient as his multiply() function but with the nice overloading of the multiply operator. Is there a way to implement a general matrix multiply without three matrices? RVO, of course, eliminates the need to have a copy of the output matrix. Apr 30 '15 at 15:34
• @BillGreene's reference to return value optimization only avoids the need for a second temporary, but one is still needed. I'll comment on this in another answer. Apr 30 '15 at 15:47
• @BillGreene: Your example is too simple. Optimizing away some of the assignments, the creation of temporaries, etc, is possible because there are no side effects the compiler has to accommodate. In essence, you're just working on a single scalar. Try an example where instead of a single scalar, the class requires allocating and deleting memory. In this case, you have to call malloc and free and the compiler cannot optimize away pairs of them without tripping up memory monitors etc. May 1 '15 at 2:06

I'd like to refute your premise that this is “very, very inefficient”. Performance-wise, it's actually pretty irrelevant, because matrix multiplication is1 $\mathcal{O}(n^{2.8})$, whereas copying one of those matrices is only $\mathcal{O}(n^2)$. So for $n$ as big as you have, the multiplication will probably completely outweigh the copying time.

That is, unless you incur an enormous constant into the copying – which is actually not so far-fetched, because the version with copying is much more expensive in another regard: it needs way more memory. So if you end up having to swap to and from hard disk, the copying might indeed become the bottleneck. However, even if you don't copy anything yourself, a strongly parallelised algorithm may well do some copies of its own. Really, the only way to make sure not too much memory will be used in each step is to split up the multiplication in columns of stateMatrix, so only small multiplications are done at a time. For instance, you can define

class HChunkMatrix // optimised for destructive left-multiplication updates
{
std::vector<arma::cx_mat> colChunks; // e.g. for an m×n matrix,
//  use √n chunks, each an m×√n matrix
public:
...

HChunkMatrix& operator *= (const arma::cx_mat& lhMult) {
for (&colChunk: colChunks) {
colChunk = lhMult * colChunk;
}
return *this;
}
}


You should also consider whether you even need to evolute that stateMatrix as one in the first place. If you basically just want independent time-evolution of n state kets, then you might as well evolute them one by one, which is much less memory-costly. In particular if evolutionMatrix is sparse, which you should definitely check out! For this is basically just a Hamiltonian, isn't it? Hamiltonians are often sparse or approximately sparse.

1Or $\mathcal{O}(n^{2.38})$, but that one isn't practically useful.

• This is the best answer; the others miss the important point that matrix multiplication is really not the sort of thing you do in-place.
– user15247
May 1 '15 at 13:28

Modern C++ has a solution for the problem by using "move constructors" and "rvalue references".

A "move constructor" is a constructor for a class, for example a matrix class, which takes another instance of the same class and moves the data from the other instance into the new instance, leaving the original instance empty. Typically, a matrix object will have two numbers for the size and a pointer to the data. Where a normal constructor would duplicate the data, a move constructor will only copy the two numbers and the pointer, so this is really fast.

An "rvalue reference", written for example as "matrix&&" instead of the usual "matrix&" is used for temporary variables. You would declare a matrix multiplication as returning a matrix&&. By doing that, the compiler will make sure that a very cheap move constructor will be used to get the result out of the function that calls it. So an expression like result = (a + b) * (c + d) where a, b, c, d are huge matrix objects, will happen without any copying.

Googling for "rvalue references and move constructors" will find examples and tutorials.

You might also want to benchmark the operations in question in Eigen. It is the only linear algebra library (that I know of) that even attempts to optimize matrix expressions involving more than two arrays. For example, I'm pretty sure if you do $v'*M*M*M*M*M$ or $M*M*M*M*M*v$ it will do a bunch of vector-matrix multiplies, rather than matrix-matrix multiplies. And if you do (v1 + v2 - v3) it will do the computation entry-wise in one pass with no intermediate, rather than two with an intermediate.

Then again, I gather OpenBLAS has a bigger collection of architecture-specific optimizations, so Eigen may or may not be a win for you. Unfortunately, there is no linear algebra library so awesome that you don't even have to consider the others when fighting for "the last 10%" of performance. Wrappers are not a 100% solution; most (all?) of them cannot take advantage of eigen's ability to merge computations this way.

• note, there are ~application specific libraries that do fancier stuff; I think Apple's API's for image compositing do things similar to what eigen does, plus mapping the computation onto the GPU... And I imagine audio stream libraries make similar optimizations... May 1 '15 at 7:38