Summing matrix products

Question

I need to compute a double sum over a weighted matrix product: $L[M]=\sum_{i,j}^{N}\Lambda[[i]]\;\omega[[i]].M.\omega^{\dagger}[[j]]$.

$\Lambda$ is a list of with N complex values(weights) and $\omega$ is a list of N matrices, so that each different $i$ corresponds to different a matrix. I need to compute this function many times with different M matrices, so I need to optimise the execution time considerably. The final result i am after is a matrix with the elements

$\mathcal{L}_{j,l}=Tr[ M_{j}^{}\dagger]L[M_{l}]]$

The matrix $\mathcal{L}$ is large ($80^2\times 80^2$), and I have parallelised my code over $j$ and $l$, so I cannot parallelise the $L$ function. As I need to compute the value of $L[M_{i}]$ 6400 I really need a fast way to compute it. One thing which I hoped could speed things up is the fact that $M$ matrices are zero matrices with only one element different from 0, which is then 1. I have not succeeded in coding this an efficient way in mathematica as the following takes 2.5s to run on my current machine.

n = 90;
Λ = RandomReal[1, n];
ω = Table[RandomReal[1, {80, 80}], {i, 1, n}];
M = IdentityMatrix[80^2];
M1 = Partition[M[[All, 1]], 80];
M2 = Partition[M[[All, 2]], 80];
LM[M_] := Sum[Λ[[i]] ω[[i]].M.ConjugateTranspose[ω[[j]]], {i, 1, n}, {j, 1, n}]
AbsoluteTiming[LM[M1];]

Edit

As suggested by J.M. I have tried to "vectorise" the sum. In theory I think this should work: $A=(c_1 \omega_1.M\quad c_2\omega_2.M\quad c_3\omega_3.M\quad...)$

$B=\begin{bmatrix}\omega_1 &\omega_2&\omega_3&...\\ \omega_1 &\omega_2&\omega_3&...\\ \omega_1 &\omega_2&\omega_3&...\\ .&.&.&.\\ .&.&.&.\\ \end{bmatrix}$

$L=Total[A.B]$

But doing this the obvious way in mathematica gives me a wrong result. My implementation looks like this:

(*simple test arrays*)
w = {{{1, 2}, {3, 1}}, {{1, 1}, {1, 2}}};
c = {2, 3};
M = {{1, 0}, {0, 0}}; 

(*Correct results*)
Correct = Sum[c[[i]] w[[i]].M.w[[j]], {i, 1, Length[c]}, {j, 1, 2}]//MatrixForm

(*Vectorisation of double sum*)
A = Table[c[[i]] w[[i]].M, {i, 1, Length[c]}];
B = Table[w[[i]], {j, 1, 2}, {i, 1, Length[c]}];
VecRes = Total[A.B]

'VecRes' is a completely different format ({2,2,2,2}) than the expected output which should be a $2\times 2$ matrix and not an array of matrices. As Michael points out this is due to the indexing being wrong, so flatten could be used to fix.

Vectorisation Remarks

Huge thanks to Michael for the thorough description! Having incorporated Michael Weyrauch example for my entire problem have really shown the power of the vectorisation. What previously took 2400s with 4 cpu now take a single cpu 10s! The only drawback is that I do not know how to effectively parallelise this algorithm compared to the previous sum which was embarrassingly parallel. However Michaels algorithm is still much quicker, but if it can be parallelised efficiently it would be a really nice!

Memory effects

This procedure seems to take up a huge amount of memory. So much in fact that it easily uses all of the 65gb available on 16 core machine and thereby shutting it down. Any clever way of avoiding this? An example of a line that uses all the memory:

Lx1 = Total[Flatten[(La*om).M.
Conjugate[Flatten[om, {{3}, {1}, {2}}]], {{1}, {4}, {3}, {2}, {5}}], 2];

Where the dimensions are: Dim[La]={1561}, Dim[om]={1561,80,80}, Dim[M]={80,6400,80}.

So La is a list of complex numbers, om is a tensor with 1561 different 80x80 matrices, and similar for M except it consists of 6400 matrices.

Answer 1

Here's a redefinition of the LM[] function in the OP:

LM[M_] :=
   Total[DiagonalMatrix[Λ].Outer[Dot, ω.M, ConjugateTranspose[ω, {1, 3, 2}], 1], 2]

On my limited set of tests with much smaller matrices and smaller lists of matrices, this version runs slightly faster than the Sum[]-based solution. You will have to do your own tests on your larger matrices.

Answer 2

For such tasks you need to use the appropriate functions, which are optimized for such purposes. In particular, these are Dot[], Flatten[], and Total[]. Never use looping constructs such as Do or Sum for such purposes if performance is an issue.

So, I would create the the N complex valued matrices $\omega$ as well as the other objects as appropriately dimensioned tensors making sure that they are "packed arrays" for efficient numerics. Here I do this using RandomComplex, which automatically packs.

n = 10;  
om = RandomComplex[{-10. - 10. I, 10. + 10. I}, {n, 100, 100}];
bM = RandomComplex[{-10. - 10. I, 10. + 10. I}, {20, 100, 100}];
lam = RandomComplex[{-10. - 10. I, 10. + 10.*I}, {n}];

With these tensors one could program your problem as

zw = (lam*om).Flatten[bM, {{2}, {1}, {3}}].
              Flatten[om, {{2}, {1}, {3}}];
L = Total[Flatten[zw, {{1}, {4}, {2}, {3}, {5}}], 2];
res = Tr[Flatten[Conjugate[Flatten[bM, {{3}, {1}, {2}}]].L,
           {{1}, {4}, {2}, {3}}], Plus, 2];

which runs rather fast. res should contain the matrix you want to produce. The various Flatten are necessary in order to bring the indices of the tensors into the appropriate order. Note that Dot always contracts the last index into the first index.

For more details check the documentation and just take the proposed code as an example which you adapt as needed to your particular problem.

As the OP noted the above code takes too much memory and does not take into account that the bM tensors are very sparse. It was meant to illustrate the use of Dot[], Flatten[], and Total[].

The problem at hand may be implemented taking note of its special structure as follows (using SparseArrays):

(*z = 1561;*)  z = 30;
om = SparseArray[
RandomComplex[{-10. - 10. I, 10. + 10. I}, {z, 80, 80}]];
lam = SparseArray[RandomComplex[{-10. - 10. I, 10. + 10.*I}, z]];
omt = Total[om, 1];
lamomt = Total[lam*om, 1];

(*y = 6400;*) y = 20;
bm = SparseArray[
Table[SparseArray[{{RandomInteger[{1, 80}], 
      RandomInteger[{1, 80}]} -> 1}, {80, 80}], y]];
bm = Flatten[bm, {{2}, {1}, {3}}];

(I use z=30 and y=20 instead of the requested z = 1561 and y = 6400 because I work on a laptop with little memory.) Note that I separated the sums over $i$ and $j$ which appears to be possible here.

With these building blocks one can obtain the desired matrix as follows:

term = (Flatten[Conjugate[bm], {{3}, {2}, {1}}].lamomt).(bm.omt);
res = Total[Table[term[[i, All, All, i]], {i, 80}], 1];

Note that I avoided to use Tr, which would have required to reorganize the indices using Flatten, which would have taken a lot of memory.