# Quickly sum over vector outer products

I am trying to do a simple summation over outer products. I wrote the code as straightforwardly as I could, but perhaps there is a more efficient way to do it.

{mat1,mat2} =
Import["https://dl.dropboxusercontent.com/s/d9m7mgd4gio78t4/12SeptemberData.nc",
{"Datasets",{"matrix1","matrix2"}}];
{sdim,ldim}=Dimensions[mat1];
dim = sdim + ldim;

DM = ConstantArray[0, {dim,dim}];
kcounter=0;
Do[
If[kcounter != k, Print[k]; kcounter += 1;];
(*The evaluation seems to take forever,
so this just spits out the current value of the k iterator
so I can keep track of what is going on.*)
If[b != a && Chop[mat1[[k, a]] mat1[[k, b]]] != 0 ,
DM += mat1[[k, a]] mat1[[k, b]] Outer[Times, mat2[[a+sdim]], mat2[[b + sdim]]]];

,{k, sdim}, {a, ldim}, {b, ldim}];


Before I added the term where it checks to make sure the necessary term is not zero, it was taking a lot longer. But still, I know there are a lot of Mathematica functions I don't know about, and if there were a way to speed up this computation that would be fantastic.

• As a side note, you can use Monitor[Do[..], k] to monitor k. – Michael E2 Sep 12 '13 at 23:02
• Are you aware that the first sdim rows of your matrix mat2 are never used anywhere? – bill s Sep 12 '13 at 23:05
• @MichaelE2, that is great to know! – Jason B. Sep 12 '13 at 23:28
• @bill_s, those rows get used in the next calculation, which is very similar to this one. – Jason B. Sep 12 '13 at 23:28

If I understand your code correctly than it can be reduced to

mat3 = mat2[[sdim + 1 ;; -1]];
DM = Transpose[mat3].Transpose[mat1].mat1.mat3;


Or in a more efficient way

DM = Transpose[#].# &[mat1.mat3];


Update:

You can take into account a!=b in the following way

DM = Transpose[mat3].(# - DiagonalMatrix@Diagonal[#]) &[Transpose[mat1].mat1].mat3;

• This almost does it! If I comment out the b!=a part of the If statement, then this matches my result exactly, and does it 20,000 times faster. – Jason B. Sep 12 '13 at 23:30
• Maybe using this won't be too slow: DM=(1 - IdentityMatrix[dim])Transpose[#].#&[mat1.mat3] – Timothy Wofford Sep 13 '13 at 6:42
• @TimothyWofford Unfortunately it doesn't match the original problem. See my update. – ybeltukov Sep 13 '13 at 12:17
• @JasonCasamir Oh, I missed this condition, see my update! – ybeltukov Sep 13 '13 at 12:19
• @JasonCasamir Matrix multiplications are not needed to be compiled. They are already implemented as fast as possible. – ybeltukov Sep 13 '13 at 20:41