fastest way to perform the following summation

Question

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:

MapThread[f, {listA, listX, listY}] 

Thanks in advance!

Answer 1

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

Answer 2

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.

Answer 3

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}
*)