Message Boards Message Boards

Reduce my computing time in matrix multiplication?

Posted 4 years ago

Hello. I needed to do the calculation:

p = DiagonalMatrix[{Exp[-I \[Alpha]], 1, Exp[I \[Beta]]}]; ppc = 
 DiagonalMatrix[{Exp[I \[Alpha]p], 1, Exp[-I \[Beta]p]}];
phase = p.Transpose[ppc];
vckm = {};
Do[{
   Do[{
     vckm =  Append[vckm, Simplify[Transpose[ou[[i]]]. phase.od[[j]]]]}, {j, 1, Length[od]}]}, {i, 1, Length[ou]}] // AbsoluteTiming

Since the length of my matrices, ou and od are of the order of 4000 each, the computation time is huge. Is there any alternative to do this calculation, s.t. my computation time is shortened to a feasible interval?

POSTED BY: Nikhila Nikhila
5 Replies
Posted 4 years ago

To optimize the time I first try to estimate where much of the time is being spent.

I see two things and want to know how much each contributes to the time.

You are doing Simplify 4000^2 times. I would like to know how many seconds the calculation takes with the Simplify in your code and how many seconds the calculation takes if you remove the Simplify.

You are also implementing your own iteration to do the 4000^2 steps. I wonder how many seconds with your method and with a similar built in iteration, also with and without using Simplify.

vckm=Outer[(Transpose[#1].phase.#2)&,ou,od,1];

You have not included small examples of your ou and od data. Thus I have not been able to verify that I am producing exactly the same results as your code does. So you should not spend great amounts of time calculating results until we know that I am producing compatible results. I have assumed that ou and od are lists of 3x3 matricies as a test and the results of my calculations are similar to, but not exactly the same as your code.

If the use of Simplify is necessary and is taking a large amount of the time then I do not see any obvious way to avoid that.

If the use of Outer, instead of nested Do loops and subscription, is enough faster then that may be worth working on until we get exactly the calculated results that you want.

You might Transpose each element of ou before you begin instead of doing that 4000^2 times.

I believe your phase matrix is diagonal and if so then that might make it possible to substitute some calculation taking less time than a full matrix multiply 4000^2 times.

Your code seems simple enough that I do not see other opportunities to greatly speed it up.

POSTED BY: Bill Nelson

In addition to tests and improvements suggested by Bill Nelson, it is important to replace the Append with Sow/Reap, since that former has quadratic complexity. If the Simplify can be delayed until the end, get rid of the Do[Do[...Append...]] entirely, as below. I remark that I lacked patience to sort out dimensions needed to handle the inner Transpose so I made an example that does not use it.

SeedRandom[1234];
n = 100;
ou = RandomReal[1, {n, 3}];
od = RandomReal[1, {n, 3}];

vckm = {};
Do[{Do[{vckm = Append[vckm, ou[[i]].phase.od[[j]]]}, {j, 1, 
     Length[od]}]}, {i, 1, Length[ou]}] // AbsoluteTiming

(* Out[177]= {6.44489, Null} *)

pd = Diagonal[phase];
vckm2 = Flatten[ou.Transpose[Map[pd*# &, od]]]; // AbsoluteTiming

(* Out[179]= {0.027209, Null} *)

They are substantially the same results. Now double the size. The second computation will grow quadratically, so a factor of 4 more time is used. The first goes up quartically. That's the quadratic effect of Append compounding the quadratic size growth.

SeedRandom[1234];
n = 200;
ou = RandomReal[1, {n, 3}];
od = RandomReal[1, {n, 3}];
vckm = {};
Do[{Do[{vckm = Append[vckm, ou[[i]].phase.od[[j]]]}, {j, 1, 
     Length[od]}]}, {i, 1, Length[ou]}] // AbsoluteTiming
pd = Diagonal[phase];
vckm2 = Flatten[ou.Transpose[Map[pd*# &, od]]]; // AbsoluteTiming
Max[Abs[Expand[vckm - vckm2]]]

(* Out[195]= {126.234, Null}

Out[197]= {0.099383, Null}

Out[198]= 0. *)

To ramp to your desired size, scale that first timing by a factor of 20^4, vs. 20^2 for the second one. Upshot: even if Simplify doesn't drag out your computation, Append will.

POSTED BY: Daniel Lichtblau

Thank you so much. It was really helpful. Your approach is much more technical and substantially decreased the computation time. I hope I am able to use Mathematica to a maximum extent some day, just as you did. I also needed to shorten my computation time in the following calculation. I hope you can help. Here lis1, lis2 etc are of the type lis. Most of the time is spent in the nested Do loops below.

lis

{0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 2., 3., 4., 5., \
6., 7., 8., 9., 10.}

m1 = .003; m2 = 0.4; m3 = 120.; mdiag = 
 DiagonalMatrix[{m1, -m2, m3}]; mr = {Range[3], Range[3], Range[3]}; 
mr[[1, 1]] = lis; mr[[1, 2]] = lis1; mr[[1, 3]] = lis2; 
mr[[2, 2]] = lis4; mr[[2, 3]] = lis5; mr[[3, 3]] = lis8;

count1 = 0; tr1 = 0; ou = {}; mtr = {}; fa1 = 0;
Do[{eu = mr[[1, 1, i1]],
   Do[{au = mr[[1, 2, i2]],
     Do[{fu = mr[[1, 3, i3]],
       Do[{du = mr[[2, 2, i4]],
         Do[{bu = mr[[2, 3, i5]],

           Do[{cu = mr[[3, 3, i6]], 
             m = {{eu, au, fu}, {au, du, bu}, {fu, bu, cu}},
             If [Det[m - mdiag] == 0,
              {eig = Eigenvectors[m],

               If[{Normalize[eig[[3]]], Normalize[eig[[2]]], 
                   Normalize[
                    eig[[1]]]}.Transpose[{Normalize[eig[[3]]], 
                    Normalize[eig[[2]]], Normalize[eig[[1]]]}] == 1,
                {mtr = Append[mtr, m], 
                 ou = Append[ou, 
                   Transpose[{Normalize[eig[[3]]], 
                    Normalize[eig[[2]]], Normalize[eig[[1]]]}]], 
                 tr1 = tr1 + 1},
                fa1 = fa1 + 1]},
              count1 = count1 + 1]}, {i6, 1, n}]}, {i5, 1, n}]}, {i4, 
         1, n}]}, {i3, 1, n}]}, {i2, 1, n}]}, {i1, 1, 
   n}] // AbsoluteTiming
POSTED BY: Nikhila Nikhila
Posted 4 years ago

Since the code you showed is only a part of a larger program I can only offer a few small ideas which may not make enough difference to measure.

You might test replacing

eig = Eigenvectors[m],
If[{Normalize[eig[[3]]], Normalize[eig[[2]]], Normalize[eig[[1]]]}.
    Transpose[{Normalize[eig[[3]]], Normalize[eig[[2]]], Normalize[eig[[1]]]}] == 1

with

eig = Reverse[Map[Normalize,Eigenvectors[m]]],
If[eig.Transpose[eig] == 1

and replacing

ou=Append[ou,Transpose[{Normalize[eig[[3]]],Normalize[eig[[2]]],Normalize[eig[[1]]]}]],

with

ou = Append[ou, Transpose[eig]]

If the value of n in your program were small enough and you had enough memory then it might be possible to generate your matrices and produce your results without using Do loops. This might be faster, but I cannot test this idea with the information that I have available.

Please test this carefully to make certain I have made no mistakes.

POSTED BY: Bill Nelson

Thank you so much for your reply! Although Simplify was necessary for my calculation, I guess I can put it in a different step as removing Simplify decreases my time to one fourth. I did not know about the Outer command. It seems much simpler than my nested Do loops. And yes, I limit my sample size to maybe 10 before doing it for the complete sample. The matrices ou and od are 3x3 matrices, as you guessed. I forgot to mention it earlier. The phase matrix is diagonal. I'll put it in that way. Thank you so much.

POSTED BY: Nikhila Nikhila
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract