3
$\begingroup$

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.

Improve this question
$\endgroup$
2
$\begingroup$

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 *)
Improve this answer
$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Sumit Jun 28 '15 at 7:44
  • $\begingroup$ @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. $\endgroup$ – C. E. Jun 28 '15 at 11:23
  • $\begingroup$ 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. $\endgroup$ – Sumit Jun 28 '15 at 12:14
  • $\begingroup$ @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). $\endgroup$ – C. E. Jun 28 '15 at 13:30
  • $\begingroup$ 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. $\endgroup$ – Sumit Jul 18 '15 at 12:11
1
$\begingroup$

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...

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.