3

I wish to calculate a rather dense antisymmetric matrix quite fast based on a random distribution of points (and their radii) in the three dimensional space as follows:

ClearAll["Global`*"]; 
n = 12; size = n*n*n;
Print["The size of the matrix is = ", size];
grid = RandomReal[{-10, 10}, {size, 3}];
epsilon = 5.2; cl = 6.;
xgrid2 = Map[First, grid];
ygrid2 = Map[(#[[2]]) &, grid];
zgrid2 = Map[Last, grid];
rad[k_, l_] := EuclideanDistance[grid[[k]], grid[[l]]]
mat1 = ParallelTable[With[{radial = rad[i, j]},
     If[radial <= epsilon, -((
       56 Max[cl - radial, 
         0]^5 (cl + 5 radial) (xgrid2[[i]] - xgrid2[[j]]))/cl^8), 0]]
    , {i, 1, size}, {j, 1, size}]; // AbsoluteTiming
AntisymmetricMatrixQ[mat1]

Unfortunately, although I used ParallelTable, it seems it requires a considerable time when the size of the matrix is high, e.g., when $n=25$ and the size would be $15625$!

So, I would be thankful if someone could provide a couple of comments in order to accelerate such a process.

M.J.2
  • 491
  • 2
  • 8
  • 1
    First comment after a first glance: built-in functions are usually faster. Instead defining rad[], you could employ EuclideanDistance. – corey979 Sep 10 '16 at 15:35
  • Thanks. This has been done. But the computational time is still high. – M.J.2 Sep 10 '16 at 15:38
  • grid are points within a unit cube, so the maximal distance that can be achieved is $\sqrt{3}\approx 1.73$; why do you test whether radial <= epsilon? It always will. – corey979 Sep 10 '16 at 15:48
  • I revised the sample code for the general case. In general, the bounds for the cube could be any number. I think the anti-symmetric feature or an application of Compile[] can help us. – M.J.2 Sep 10 '16 at 15:52
  • Employing Nearest can substantially speed it up; see the beginning of this. – corey979 Sep 10 '16 at 15:57
  • 2
    How do you want to apply Nearest for constructing the matrix here? Please write your comment in a piece of implementation. – Faz Sep 10 '16 at 16:06
  • I'm attempting to find only elements that are not zero, and define the lacking ones as zeros. – corey979 Sep 10 '16 at 16:14
  • It seems that indeed my idea to use Nearest does not make the code faster, or at least I cannot employ it efficiently. So my last remarks are: (1) get rid of xgrid2, ygrid2, zgrid2; (2) replace (xgrid2[[i]] - xgrid2[[j]]) with grid[[i, 1]] - grid[[j, 1]]), and (3) compute the ParallelTable only for {j,i+1,size}, pad the rest with zeros and add to it a transposed and multiplied by $-1$ matrix. This should speed up the computation by a factor of two (the whole code for n=25 run 6 mins, so it easily be reduced to 3 mins). – corey979 Sep 10 '16 at 16:55
  • Please submit your third comment as an answer. – M.J.2 Sep 10 '16 at 17:02

2 Answers2

5

We can speed this up a lot by using that radial <= epsilon means only a small fraction of all those matrix elements have to be considered. Here's how I would do that:

computeMatrix[grid_, xgrid2_, size_] := Module[{dist, diff, eliminate, pos, values, gr = ConstantArray[xgrid2, size]},
  dist = DistanceMatrix[grid, DistanceFunction -> EuclideanDistance];
  diff = gr - Transpose[gr];
  eliminate = UnitStep[ConstantArray[epsilon, {size, size}] - dist];
  pos = SparseArray[eliminate]["NonzeroPositions"];
  dist = Extract[dist, pos];
  diff = Extract[diff, pos];
  values = ((56 Clip[cl - dist, {0, \[Infinity]}]^5 (cl + 5 dist) diff)/cl^8);
  SparseArray[pos -> values]
]

mat1 == computeMatrix[grid, xgrid2, size]
(* True *)

This evaluates in approximately 0.1 seconds on my computer, more than 100 times faster than the original solution.

Marius provided the vectorized function in the second to last line in a comment below. (If there is interest, one can see this solution without with a non-vectorized function in the edit history of this post.)

Note: This is a very memory hungry function for large $n$. A typical symptom for lack of memory is that the operating system is sluggish or non-responsive while the computation is ongoing. I might update this later with an improvement directed at memory efficiency, but if memory is a problem is might also be a good idea to ask a question about how this can be optimized with regards to memory. It also appears that the speed up with the vectorized approach as compared to my previous approach is lost when the computation is memory constrained.

Community
  • 1
C. E.
  • 70,533
  • 6
  • 140
  • 264
  • An amazing response. Thanks. Just one query. When I choose $n=25$, your code (without compiling) makes my laptop to get frozen! May you please report the whole time required for the case $n=25$ in your system? – M.J.2 Sep 10 '16 at 18:07
  • @M.J.2 It takes about 75 seconds on my computer. If the laptop freezes it may be due to lack of memory, so you may have to take that into consideration as well. – C. E. Sep 10 '16 at 18:16
  • 1
    How about using a vectorized function? Just uncompiled values = ((56 Clip[cl - dist, {0, \[Infinity]}]^5 (cl + 5 dist) diff)/cl^8) gives the same result as your MapThread and takes 0.024 seconds. Your's give 0.7 on my machine as well. – Marius Ladegård Meyer Sep 10 '16 at 23:14
  • @MariusLadegårdMeyer Yup! This is what I wanted to do, never got to it. I updated the post to reflect this great improvement. – C. E. Sep 11 '16 at 00:04
  • @M.J.2 With the new function it took me 106 seconds to compute for $n=25$ and I maxed out my 16 GB memory. If you have less than 16 GB memory, I'm not sure if you can run it and if you can run it the speed will suffer I guess. – C. E. Sep 11 '16 at 00:30
  • This modified response is really great. Yes, I am using 8 GB of RAM in my machine and I think because of that my machine freezes for the cases higher than $n=23$. The most memory-consuming portion of the code is gr = ConstantArray[xgrid2, size] and gr - Transpose[gr] at the beginning. – M.J.2 Sep 11 '16 at 12:12
1

While C. E.'s answer is outperforming my solution, I thought I'll elaborate on my last comment to the question for the sake of completeness.

First, xgrid2, ygrid2, zgrid2 are not needed, and the differences xgrid2[[i]] - xgrid2[[j]] might be written as grid[[i, 1]] - grid[[j, 1]]. Next, With is also not necessary. And finally, the main point is that the computed matrix is antisymmetric by construction, so one only needs to compute half of it:

mat1 = ParallelTable[
    If[rad[i, j] <= epsilon, 
        -((56 Max[cl - rad[i, j], 0]^5 (cl + 5 rad[i, j]) (grid[[i, 1]] - grid[[j, 1]]))/cl^8),0],
{i, 1, size}, {j, i + 1, size}]; // AbsoluteTiming

mat2 = PadLeft[mat1, {size, size}]; // AbsoluteTiming
mat3 = mat2 - Transpose@mat2; // AbsoluteTiming

The above is overall two times faster than the original code; mat1 takes 2/3 of the total time, while mat3 takes the remaining (the Transpose appears to be quite time-consuming; mat2 requires a negligible amount of time. So in fact the crucial information that is contained in mat1 can be obtained 3x faster than with the original code.

Community
  • 1
corey979
  • 23,947
  • 7
  • 58
  • 101
  • Why the output matrix mat3 in your response is different to the original output matrix mat1 in the question? – M.J.2 Sep 11 '16 at 12:18
  • 1
    It's not; I just run it to verify. Note that I also used mat1 to name a matrix, but it's different than yours. Maybe you're overriding definitions and in fact comparing my mat3 to my mat1, which are of course different. Change mat1 in your code to mat and run mat==mat3. – corey979 Sep 11 '16 at 12:23
  • Yes, I made a mistake. The matrices are the same. Do you know why Transpose is slow for large matrices! Furthermore, With[] is required, since in your way, rad[i, j]is computed more than once for each entry! Am I right? – M.J.2 Sep 11 '16 at 12:26
  • 1
    No and it puzzles me, especially when you take m = UpperTriangularize@RandomReal[1, {n^3, n^3}]; trm = m - Transpose@m; // AbsoluteTiming - so the same operation as in the code, the same size, reals as elements etc. - it works much faster for such random matrix. – corey979 Sep 11 '16 at 12:37