3
$\begingroup$

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.

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ First comment after a first glance: built-in functions are usually faster. Instead defining rad[], you could employ EuclideanDistance. $\endgroup$
    – corey979
    Commented Sep 10, 2016 at 15:35
  • $\begingroup$ Thanks. This has been done. But the computational time is still high. $\endgroup$
    – M.J.2
    Commented Sep 10, 2016 at 15:38
  • $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Sep 10, 2016 at 15:48
  • $\begingroup$ 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. $\endgroup$
    – M.J.2
    Commented Sep 10, 2016 at 15:52
  • 2
    $\begingroup$ How do you want to apply Nearest for constructing the matrix here? Please write your comment in a piece of implementation. $\endgroup$
    – Faz
    Commented Sep 10, 2016 at 16:06

2 Answers 2

Reset to default
5
$\begingroup$

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.

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – M.J.2
    Commented Sep 10, 2016 at 18:07
  • $\begingroup$ @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. $\endgroup$
    – C. E.
    Commented Sep 10, 2016 at 18:16
  • 1
    $\begingroup$ 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. $\endgroup$
    – Marius Ladegård Meyer
    Commented Sep 10, 2016 at 23:14
  • $\begingroup$ @MariusLadegårdMeyer Yup! This is what I wanted to do, never got to it. I updated the post to reflect this great improvement. $\endgroup$
    – C. E.
    Commented Sep 11, 2016 at 0:04
  • $\begingroup$ @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. $\endgroup$
    – C. E.
    Commented Sep 11, 2016 at 0:30
1
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Why the output matrix mat3 in your response is different to the original output matrix mat1 in the question? $\endgroup$
    – M.J.2
    Commented Sep 11, 2016 at 12:18
  • 1
    $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Sep 11, 2016 at 12:23
  • $\begingroup$ 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? $\endgroup$
    – M.J.2
    Commented Sep 11, 2016 at 12:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – corey979
    Commented Sep 11, 2016 at 12:37

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.