Super C++ optimization of matrix multiplication with Armadillo

Question

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?

Answer 1

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.

Answer 2

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

Answer 3

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.

Answer 4

I'd like to refute your premise that this is “very, very inefficient”. Performance-wise, it's actually pretty irrelevant, because matrix multiplication is¹ $\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.

---

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

Answer 5

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.

Answer 6

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.