# fastest way to perform the following summation

I have a vector $A$ with components $A_{i}$, and two matrices $X$ and $Y$ with components $X_{ij}$, and $Y_{ij}$ respectively. What is the fastest way of performing the following summation

$$A_{i}A_{j}X_{ik}X_{jl}Y_{kl}$$

where I would like to sum over the repeated indices?

The problem I have is that I need to perform this summation thousands of times for vectors of length 100, matrices 100×100 etc. So I want to avoid using Sum. I think once I know how to do this bit I am happy with how to then apply this to a list of such objects - I guess I just create a compiled function f that performs the above summation then use this with MapThread. i.e.: Say I have lists of the above objects listA, listX and listY, with length 1000. Then provided I can define a function f that performs the summation above, all I would need to do is:

• Can you show the (slow) code you have for the Sum so that we can understand the exact calculation you want to make? Jan 31, 2015 at 19:56

Since your central question was about speed (or time complexity), you might wish to know an important result from elementary theory of algorithms and computational complexity, which is that the time and space complexity of matrix multiplication depends upon the order of such multiplication (see Introduction to Algorithms by Cormen, Leiserson and Rivest). As such, your timing will depend upon the order of multiplication (specified by grouping) and the dimensions of your vectors and matrices. Here is an example:

A = Table[RandomReal[], {3000}];
X = Table[RandomReal[], {3000}, {8000}];
Y = Table[RandomReal[], {8000}, {8000}];

Timing[A.X.Y.X\[Transpose].A]

(* {0.956274, 1.79343*10^13} *)

Timing[(A.X).Y.(X\[Transpose].A)]

(* {0.544384, 1.79343*10^13} *)

Timing[A.(X.Y).(X\[Transpose].A)]

(* {36.569151, 1.79343*10^13} *)

So whereas all these give the same answer, one is 1/67 times as fast as another, which in turn is twice as fast (for these dimensions) than the solution provided by @DumpsterDoofus. If you know your matrix sizes ahead of time, you can choose the computation order optimally, as given in Theory of Algorithms. @DumpsterDoofus is correct, though, that for dimensions around 1000, all the results should be fast enough for your purposes.

Note too, that if you have very large matrices, you may want to Parallelize[] your computation and here too, the groupings can affect the timing (and space complexity, or memory usage) significantly.

## Why does the order of multiplication matter?

Because several comments expressed an interest in understanding the roots of the importance of proper grouping of matrices during multiplication, let me give a simple example to illustrate. Let $A_{p \times q}$, $B_{q \times r}$ and $C_{p \times r}$ be three matrices with the dimensions given, and we seek to compute $ABC$. Consider $A_{10 \times 100}$, $B_{100 \times 5}$ and $C_{5 \times 50}$.

If we multiply according to $(AB)C$ we perform $10 \cdot 5 \cdot 5 = 5000$ scalar multiplications for $(AB)$, then an additional $10 \cdot 5 \cdot 50$ to then multiply by $C$: total = 7500 multiplications.

If we multiply according to $A(BC)$ we perform $100 \cdot 5 \cdot 50 = 25000$ scalar multiplications for $(BC)$ then an additional $10 \cdot 100 \cdot 50 = 50000$ to then multiply by $A$: total = 75000 multiplications. A factor of 10 higher than in the other case.

There are smart linear programming methods to determine the optimal (i.e., least complex) groupings, but that would take us far afield here...

• Nice. Couple comments: on my computer and using your definitions of A, X and Y, Timing[A.X.Y.X\[Transpose].A] is actually slightly faster, taking 0.12 seconds compared to 0.15 seconds for Timing[(A.X).Y.(X\[Transpose].A)]. I'm not sure why this is the case, seeing as you found that the former was twice as slow as the latter. Second, good job pointing out that grouping order matters! I hadn't considered it at all. Third, Parallelize will give no speedup, since LAPACK calls are automatically parallelized across available CPU cores, and Dot calls LAPACK. Jan 31, 2015 at 21:31
• @DumpsterDoofus: The above timing results were on my MacBook Air laptop at home, not my powerful Mac desktop at my research lab. I may re-test when I return to a real work machine, but I believe that smart parallelization could possibly speed up the overall computation by splitting the multiplications optimally. If the matrix operations are $A.B.C.D.E.F.G.H$ I suspect (but do not know) that the proper partitioning (e.g., $(A.B).(C.D).(E.F).(G.H)$ might lead to faster execution on a multi-core machine. Feb 1, 2015 at 0:13
• This question turned out a great deal more interresting than I was expecting. Thanks a lot for this! Feb 2, 2015 at 14:25

In Einstein notation your initial expression is $$A_{i}A_{j}X_{ik}X_{jl}Y_{kl}\\=A_iX_{ik}Y_{kl}X^\mathsf{T}_{lj}A_j\\=AXYX^\mathsf{T}A$$

which can be easily entered into Mathematica as

A.X.Y.X\[Transpose].A

It's unlikely you'll be able to devise anything faster than this, but this should be very fast. On my machine, this executes in 0.062 seconds for $1000\times1000$ random real matrices, so performance will not be a concern.

Here's another way that's quite a bit faster on David Stork's example: #.Y.# &[A.X]

SeedRandom[0];
A = Table[RandomReal[], {3000}];
X = Table[RandomReal[], {3000}, {8000}];
Y = Table[RandomReal[], {8000}, {8000}];

AbsoluteTiming[A.X.Y.X\[Transpose].A]
AbsoluteTiming[(A.X).Y.(X\[Transpose].A)]
AbsoluteTiming[A.(X.Y).(X\[Transpose].A)]

AbsoluteTiming[#.Y.# &[A.X]]
(*
{0.198840, 1.76089*10^13}
{0.180308, 1.76089*10^13}
{5.247255, 1.76089*10^13}

{0.038420, 1.76089*10^13}
*)
• Wow! Do you mind explaining to me why this last one is so fantastically fast? Feb 2, 2015 at 14:23
• @user12876 I don't really know. Note AbsoluteTiming[(A.X).Y.(A.X)] is nearly as fast. My code does A.X only once (instead of twice), and I expected that to be the only speed-up. It could be that vector.matrix.vector is highly optimized/parallelized, but then the second option in my answer (from David Stork) should be fast, since X\[Transpose].A == A.X is True. Dot is Flat, but I don't think that is relevant. (I mention it in case someone can point out that my assumption is wrong.) Feb 2, 2015 at 15:10