5

I have a given SparseArray, let's call it A.

I want to construct another SparseArray, let's call it B, which has the same "shape" as A. Right now, I construct B by applying a certain defining function on the positions of the elements of A. So something like

B = SparseArray[ Map[ some_function[], A["ExplicitPositions"] ] ]

It works OK by it turns out to be quite slow for large SparseArray...

Is there another more efficient method maybe ?

Thanks in advance.

coussin
  • 361
  • 1
  • 8

2 Answers2

5

What makes building a SparseArray from triples {i,j,a} (explicit positions and explicit values) so expensive is the sorting of the indices. But once the matrix is constructed, then the indices are sorted. You can exploit this by using an undocumented constructor that allows to initialize a matrix only by the "RowPointers", "ColumnInidces" and , "ExplicitValues". See this example, where I use a Kirchoff matrix of a random matrix as initial sparse matrix.

A = KirchhoffMatrix[RandomGraph[{1000, 2000}]];

B = With[{ data = {Automatic, A["Dimensions"], A["Background"], {1, {A["RowPointers"], A["ColumnIndices"]}, A["ExplicitValues"]}} }, SparseArray @@ data ];

A == B

True

So in your application you might want to try

B = With[{
    data = {Automatic, A["Dimensions"], A["Background"],
      {1, {A["RowPointers"], A["ColumnIndices"]}, 
       Map[f, A["ExplicitPositions"]]}}
    },
   SparseArray @@ data
   ];
Henrik Schumacher
  • 106,770
  • 7
  • 179
  • 309
  • Thank you very much for your answer. In my case, all of the time was spend in the function defining the elements so there was only a marginal gain of time constructing the SparseArray from triples {i,j,a} or from {RowPointers,ColumnIndices}.

    But this discussion was nevertheless very useful for me :)

    – coussin Dec 19 '23 at 11:10
  • You're welcome! – Henrik Schumacher Dec 19 '23 at 12:01
3

Is this even necessary to improve? Let's define a sparse array $A$,

A = SparseArray[{3, 2} -> 6, {10, 20}];

which has rules

ArrayRules[A]
(*    {{3, 2} -> 6, {_, _} -> 0}    *)

Now define a function that we want to map on the elements of $A$, defining it so that it logs every call:

f[x_] := Echo[g[x]]

Map it onto the elements of $A$:

B = Map[f, A, {2}];

(* g[6] ) ( g[0] *)

Note that $f$ was only called twice. The result is

ArrayRules[B]
(*    {{3, 2} -> g[6], {_, _} -> g[0]}    *)

There's no indication that "too much work" was done in constructing $B$.

Roman
  • 47,322
  • 2
  • 55
  • 121
  • Yeah, when I say "it's quite slow" I guess it's just the time to call the defining function times the length of A["ExplicitPositions"]... I guess I was expecting some kind of internal vectorization by using another method which do not map[] blindly over the list of ExplicitPositions... – coussin Dec 18 '23 at 11:42
  • Can you be more specific about the kind of magic you expect? – Roman Dec 18 '23 at 12:19
  • Huh, what I rote this morning was totally unreadable. Must have been due to jetlag... So here the correct version: Roman, I think the issue here is that user coussin wants to compute new "ExplicitValues" from "ExplicitPositions". In your example, the new "ExplicitValues" are just a function of the original "ExplicitValues", so the function can be simply mapped over the whole SparseArray. – Henrik Schumacher Dec 18 '23 at 21:41