# Fast approach to do Summation in Compile[]?

My code does massive Summation and Matrix multiplication

Compile[] has boosted it distinctly. But I read some literatures related to my program, It seems there are approches to make it even faster. Maybe it can be improved from optimizing MMA language or algorithm itself.

My code is below.

MomentComputing =
Compile[{{Mmax, _Integer}, {k, _Integer}, {image, _Real,
2}, {xLegendreP, _Real, 2}, {yLegendreP, _Real, 2}},
Block[{m, n, width, momentMatrix},
width = Length[image];
momentMatrix = Table[0., {Mmax + 1}, {Mmax + 1}];
Do[
momentMatrix[[m + 1, n + 1]] = ((2. (m - n) + 1.) (2. n + 1.)/((k k)*width width)) xLegendreP[[
m - n + 1]].image.yLegendreP[[n + 1]], {m, 0, Mmax}, {n, 0,
m}];

momentMatrix], CompilationTarget -> "C",
RuntimeAttributes -> {Listable}, Parallelization -> True,
RuntimeOptions -> "Speed"]


It should be better if I don't use any loop operations. But I can not figure out any other approaches. Probably matrix vector multiplication should be time-consuming as well.

• As a general comment, matrix-vector multiplication is likely to be most efficient. Unless you discard part of the result. It's built on LAPACK and highly optimized. – Michael E2 Jul 26 at 17:23
• What happens if you omit momentMatrix & its initialization and change the Do loop to return Table[((2. (m - n) + 1.) (2. n + 1.)/((k k)*width width)) xLegendreP[[m - n + 1]].image.yLegendreP[[n + 1]], {m, 0, Mmax}, {n, 0, m}]? Ditto for reconstructedImage. – Michael E2 Jul 26 at 17:26
• Including working test code would help you get better answers. – Michael E2 Jul 26 at 17:28
• @MichaelE2 I just tried the your suggestion. Theorietically it should be faster, but it didn't make differences well. Maybe it is fast enough and can not be boosted any more. Thanks anyway. – PalvinWang Jul 26 at 18:41
• Might be faster to do xLegendreP.image.Transpose[yLegendreP] and transpose, then shift row j by j-1 places rightward, then transpose back. Compute those multipliers separately (precomputing k^2*widh^2) and multiply matrices component-wise. – Daniel Lichtblau Jul 26 at 20:30

First we run the original with an example (which should have been in the post).

Mmax = 400;
W = 1024;
deltaX = .1;
SeedRandom[1111];
lambdaMatrix = RandomReal[1, {W, W}];
lPoly = DeveloperToPackedArray[RandomReal[{-10, 10}, {Mmax + 1, W}]];

AbsoluteTiming[
res = MomentComputing[Mmax, 5, lambdaMatrix, lPoly, lPoly];]

(* Out[52]= {12.1522, Null} *)


What I had in mind is this (I could combine some steps and maybe improve speed a bit further).

MomentComputing2 =
Compile[{{Mmax, _Integer}, {k, _Integer}, {image, _Real,
2}, {xLegendreP, _Real, 2}, {yLegendreP, _Real, 2}},
Block[{m, n, width, mult, mat1, momentMatrix},
width = Length[image];
mult = 1./(k^2*width^2);
mat1 =
Table[(2. (m - n) + 1.) (2. n + 1.), {m, 0, Mmax}, {n, 0, Mmax}]*
mult;
momentMatrix = Transpose[xLegendreP.image.Transpose[yLegendreP]];
momentMatrix =
Transpose[
MapIndexed[(PadLeft[Drop[#1, -(#2[[1]] - 1)], Mmax + 1, 0.]) &,
momentMatrix]];
mat1*momentMatrix
], CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];


Repeat the example.

AbsoluteTiming[
res2 = MomentComputing2[Mmax, 5, lambdaMatrix, lPoly, lPoly];]

(* Out[54]= {0.023139, Null} *)


Check that the results are within numeric fuzz:

In[57]:= Max[Abs[res - res2]]

(* Out[57]= 3.69482*10^-13 *)

Upshot: machine precision matrix products will beat iterated vector products by a wide margin. To the point where it did not matter that I discard nearly half the computed elements.

• Sorry to reply late. Your solution is amazing. But I dont understand the meaning of Drop[#1, -(#2[[1]] - 1)] inside your code. Could you give explanation on it? – PalvinWang Jul 31 at 15:48
• For MapIndexed over a list, the list position is given by that #2[[1]]. We use the list position to determine how much to shift rightward and replace by zeros (the Drop and PadLeft` do this). – Daniel Lichtblau Jul 31 at 15:54