Collecting algebraic performance tricks

General Usage Performance
1

There are many ways to dramatically improve performance by using some algebraic equivalencies, for example :

tr(A*B) == dot(A,B') 
diag(A*B) == vec(sum(A.*B',dims=2))

Do you know if there are ressources listing this kind of useful equivalencies?
And if not how about starting an organized collaborative document where everyone can share his knowledge?

Edit : corrected formula for diag

2 
@article{petersen2008matrix,
  title={The matrix cookbook},
  author={Petersen, Kaare Brandt and Pedersen, Michael Syskind and others},
  journal={Technical University of Denmark},
  volume={7},
  number={15},
  pages={510},
  year={2008}
}

You can find some updates online up to 2012 I think.

3 Likes
3

The Matrix Cookbook is indeed the most lifesaving ressource for matrix calculus!
However when it comes to implementation it does not help so much. I was thinking of something more concrete and code oriented.

4

Math textbooks? Typically, this is the job of mathematicians, so the language for identities is usually in written in pure math, and not a specific programming language. This is why you won’t find many textbooks that give identities in a specific computer programming language.

However, you should also be able to find a whole bunch of identities like this in any sort of computer algebra term rewrite system.

2 Likes
5

sum(A .* B', dims = 2)

2 Likes
6

I think that in the context of Julia the best approach is to define methods to do what you mean, without worrying about the implementation details, and then depending on the details, define optimized methods when it makes sense.

This approach establishes a clean separation between intent and implementation, and you can also use the generic method to test the specific one.

Regarding matrix indentities: the optimal approach can depend on structure (sense, n-diagonal, triangular, sparse, block, …). Also, sometimes the a highly optimized generic BLAS approach will beat a theoretically superior native one — one should always benchmark.

1 Like
7 
tr(A*B) == dot(A,B') 
diag(A*B) == vec(sum(A,B',dims=2))

This sort of matlab-style trickery might not be that useful in Julia where you can just write loops. Exceptions are when you can recognize a BLAS2 or BLAS3 operation.

8
9

I get a “not found” for this.

10

what if you click the pdf link here ? The Matrix Cookbook

11

Yes, this one works, thanks.

This is the version I had, and the latest I know of. It is a pity that they did not upload the source to a repo somewhere with an open license.