3

Here is the example/problem.

I have a array of $n$ random number.

n = 2000;
p = RandomReal[{0, 9}, n];

I have to create a $n\times n$ sparse matrix with $(i,j)$th element given by this condition

$m_{i,j} = \begin{array}{ll} a(p_i-p_j)&2<|p_i-p_j|<3 \\ b(p_i-p_j)&5<|p_i-p_j|<6 \end{array}$

A brutal way to do it is using If or Piecewise

f1[x_, y_] := If[2 < Abs[x - y] < 3, a[x-y], If[5 < Abs[x - y] < 6, b[x-y], 0]]
SparseArray@Table[f1[p[[i]], p[[j]]], {i, n}, {j, n}]; // AbsoluteTiming

{24.858066, Null}

f2[x_, y_] := Piecewise[{{a[x-y], 2 < Abs[x - y] < 3}, {b[x-y], 5 < Abs[x - y] < 6}}]
SparseArray@Table[f2[p[[i]], p[[j]]], {i, n}, {j, n}]; // AbsoluteTiming

{26.741343, Null}

And of course Map works much faster

SparseArray@Partition[
If[2 < # < 3, a[#], If[5 < # < 6, b[#], 0]] & /@ 
 Flatten[Abs[# - p] & /@ p], n]; // AbsoluteTiming

{7.797710, Null}

(I don't know how to use Map in this kind of 2D array so I combined Flatten and Partition.)

However I can Parallelize the Table and reduce the time to ~1/8th (I have to use the Workstation in my office) and apparently it looks like the If will be a winner. ParallelMap doesn't show similar reduction rate in time.

SparseArray@ParallelTable[f1[p[[i]], p[[j]]], {i, n}, {j, n}]; // AbsoluteTiming

{13.573756, Null}

SparseArray@Partition[
ParallelMap[If[2 < # < 3, a[#], If[5 < # < 6, b[#], 0]] &,
 Flatten[Abs[# - p] & /@ p]], n]; // AbsoluteTiming

{8.580731, Null}

So what would be the optimised (to minimise time) way to do the job?

For my actual problem n~10000 and a and b are combination of trigonometric and arithmetic functions.

Sumit
  • 15,912
  • 2
  • 31
  • 73

2 Answers2

2

I replaced a[#] with Abs[#] and b[#] with Abs[#] in your fastest solution. It took 7.8 seconds to evaluate. The following takes 0.085 seconds to evaluate and gives the same result:

m = Compile[{{p, _Real, 1}},
   Outer[Which[
      2 < Abs[# - #2] < 3, Abs[# - #2],
      5 < Abs[# - #2] < 6, Abs[# - #2],
      True, 0
      ] &, p, p], CompilationTarget -> "C", RuntimeOptions -> {"Speed"}];

m[p] // AbsoluteTiming // First
(* Out: 0.084618 *)
C. E.
  • 70,533
  • 6
  • 140
  • 264
  • Thanks @Pickett, the timing is really impressive. I never used Compile before, so thanks for that also. I was wondering if it is possible to store the array as functions? For my case the a, b functions are parameter dependent and I have to produce this matrix for a particular set p and for different parameter values. With your syntax I can afford to generate the matrix for each parameter value, but I think storing the matrix in a functional form can be a bit faster! – Sumit Jun 28 '15 at 07:44
  • @Sumit Could you elaborate on what you mean by "store the array as functions"? I will try to answer the question, but I don't quite understand it yet. – C. E. Jun 28 '15 at 11:23
  • Sorry @Pickett for not making it clear. Let say a[x_]=Sin[k1 x] and b[x_]=Cos[k2 x]. So the final matrix will be something like m1[k1_, k2_] = SparseArray[{{1, 3} -> Sin[2.345 k1], {4, 2} -> Cos[5.678 k2],...}] (*= m[p]*). Now I can simply give the value of k1 and k2 and generate the whole matrix (m1[k1,k2]). This needs to evaluate perhaps few thousands terms only. – Sumit Jun 28 '15 at 12:14
  • @Sumit I see, no that's not possible because the C language can't handle symbolic expressions. It has to receive and return numeric values (or another valid C data type). – C. E. Jun 28 '15 at 13:30
  • Hi @Pickett. Is it possible to generalise your method if each element of p is a vector? If you are interested, I post it as a new question How to use Compile to generate an n×n array using n vectors. – Sumit Jul 18 '15 at 12:11
1

Not particularly memory efficient, but ~20X faster than your first, ~5X faster than your non-parallel map (which, btw, does not produce the same results as your first) on n=1000...

oo = Flatten@Outer[Subtract, p, p];
aoo = Abs@oo;
rl = (Length@p)^2;
aa = Pick[Range@rl, Unitize@Clip[aoo, {2, 3}, {0, 0}], 1];
bb = Pick[Range@rl, Unitize@Clip[aoo, {5, 6}, {0, 0}], 1];
sa = Partition[SparseArray[Join[aa, bb] -> Join[a /@ oo[[aa]], b /@ oo[[bb]]], rl],Length@p];

Just a cobbled-up idea, perhaps something in there useful to you...

ciao
  • 25,774
  • 2
  • 58
  • 139