2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ As a side note, you can use Monitor[Do[..], k] to monitor k. $\endgroup$ – Michael E2 Sep 12 '13 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 '13 at 23:05
  • $\begingroup$ @MichaelE2, that is great to know! $\endgroup$ – Jason B. Sep 12 '13 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 '13 at 23:28
8
$\begingroup$

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;
$\endgroup$
  • 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 '13 at 23:30
  • $\begingroup$ Maybe using this won't be too slow: DM=(1 - IdentityMatrix[dim])Transpose[#].#&[mat1.mat3] $\endgroup$ – Timothy Wofford Sep 13 '13 at 6:42
  • $\begingroup$ @TimothyWofford Unfortunately it doesn't match the original problem. See my update. $\endgroup$ – ybeltukov Sep 13 '13 at 12:17
  • $\begingroup$ @JasonCasamir Oh, I missed this condition, see my update! $\endgroup$ – ybeltukov Sep 13 '13 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 '13 at 20:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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