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} = 
dim = sdim + ldim;

DM = ConstantArray[0, {dim,dim}];
   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.

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

1 Answer 1


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];


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

DM = Transpose[mat3].(# - DiagonalMatrix@Diagonal[#]) &[Transpose[mat1].mat1].mat3;
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jason B.
    Sep 12, 2013 at 23:30
  • $\begingroup$ Maybe using this won't be too slow: DM=(1 - IdentityMatrix[dim])Transpose[#].#&[mat1.mat3] $\endgroup$ Sep 13, 2013 at 6:42
  • $\begingroup$ @TimothyWofford Unfortunately it doesn't match the original problem. See my update. $\endgroup$
    – ybeltukov
    Sep 13, 2013 at 12:17
  • $\begingroup$ @JasonCasamir Oh, I missed this condition, see my update! $\endgroup$
    – ybeltukov
    Sep 13, 2013 at 12:19
  • 2
    $\begingroup$ @JasonCasamir Matrix multiplications are not needed to be compiled. They are already implemented as fast as possible. $\endgroup$
    – ybeltukov
    Sep 13, 2013 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.