10

I'm trying to use a SparseArray for some large homological algebra calculations. I'd like to find indices of a SparseArray whose entries are equal to 1. If I were operating naively on the ArrayRules version of the array, I'd just need do something like Cases[sparserules, HoldPattern[_List -> 1]] and that would give me what I want. (FWIW: here's why you need the HoldPattern.)

If it helps, I am also happy with just being able to pull off one such element at a time - Cases[sparserules, HoldPattern[_List -> 1], All, 1] does what I want on ArrayRules output. The order is not important - if this returns a random entry that's equal to 1, that's also fine by me.

Unfortunately, this doesn't seem to work on SparseArray objects, and I don't want to be packing and unpacking SparseArrays all the time. Is there an efficient way to get all (or just one) indices of a SparseArray whose entries are equal to a given element?

Some possible starting points: SparseArrays are atomic; using Pick on SparseArrays. Both of these are a bit beyond my level of sophistication, however.

Community
  • 1
dvitek
  • 295
  • 1
  • 5
  • 3
    Is the array numeric? Then (SparseArray[Unitize[#1 - #2], Automatic, 1]["NonzeroPositions"]) & should fit your needs - call with array and value to index. – ciao May 03 '16 at 05:52
  • 1
    Ciao ciao, thanks for the tip on "NonzeroPositions" - that could be useful. Unfortunately the array is not numeric - mostly polynomials, so Unitize won't work. I like the trick, though. – dvitek May 03 '16 at 07:10
  • 1
    @ciao Got your answer to work - I just used Map[Boole[#===1]&, sparsearray, {2}]["NonzeroPositions"], which fixed the non-numeric problem. If you put your comment as an answer I'd be happy to accept it. – dvitek May 03 '16 at 16:46
  • I presume you wrapped that Map in SparseArray? Otherwise, it's possibly not doing what you think... and even if it is, if you're going to map over the array, might as well just MapIndexed over it, and return the index for elements that match your tests. In any case, feel free to use the ideas in a self-answer... – ciao May 03 '16 at 23:32
  • In fact you don't have to wrap Map in SparseArray - it will automatically thread over SparseArray objects, whereas MapIndexed will convert them to lists. If you want a concrete demonstration, run the following: a = SparseArray[(# -> RandomInteger[{1, 200}]) & /@ RandomInteger[{1, 2000}, {4000, 2}]]; a1 = Map[#^2&, a, {2}]; a2 = MapIndexed[#1^2&, a, {2}]; If you do this, a1 will have head SparseArray, and the command will execute quickly. But a2 will take a while and have head List. See https://reference.wolfram.com/language/tutorial/LinearAlgebraSparseArrays.html. – dvitek May 04 '16 at 00:11
  • LOL - I'm not talking about the obvious - I'm talking about the results of ["NonzeroElements"] on resultants.... – ciao May 04 '16 at 00:32
  • OK, then I'd appreciate it if you could be more explicit about what you're talking about. I certainly don't understand what you mean by "the results of ["NonzeroElements"] on resultants". – dvitek May 04 '16 at 06:45

2 Answers2

5

Here's a fast way if you have general values, like @ciao's:

SeedRandom[0];
sa = SparseArray[(# -> RandomChoice[Join[{x^2 + y^2, x^3}, Range@200]]) & /@ 
    RandomInteger[{1, 2000}, {4000, 2}]];

Extract[sa["NonzeroPositions"], 
  Position[sa["NonzeroValues"], x^2 + y^2, 1]] // RepeatedTiming

{0.00015,
  {{22, 1644}, {165, 37}, {207, 910}, {332, 291}, {354, 1432}, {997, 552},
   {1211, 1944}, {1240, 1514}, {1505, 1421}, {1632, 899}, {1735, 702},
   {1808, 1690}, {1816, 1445}, {1931, 1856}, {1934, 1709}, {1988, 312}}}
*)

If all values are numbers, then this is faster:

SeedRandom[0];
sa2 = SparseArray[# -> RandomInteger[{1, 200}] & /@ 
    RandomInteger[{1, 2000}, {4000, 2}]];

Pick[sa2["NonzeroPositions"], sa2["NonzeroValues"], 1] // RepeatedTiming
(*
{0.0000261,
  {{261, 1553}, {295, 1413}, {392, 1762}, {400, 1357}, {418, 850},
   {580, 237}, {665, 1699}, {739, 1314}, {776, 1589}, {874, 206}, {932, 828},
   {1215, 730}, {1636, 150}, {1770, 214}, {1815, 1922}, {1928, 676}}}
*)
Community
  • 1
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Big +1 - that's the way to do it... – ciao May 04 '16 at 12:40
  • Ah, great, thanks a bunch! I didn't know about NonzeroValues either - this doesn't appear to be documented. And if I only want to pick off particular expressions one element at a time, it's just one more argument to Position! I'm still curious why the Pick method is not faster on non-numerics, though - interesting. – dvitek May 04 '16 at 15:40
  • @dvitek Pick only works when the expression tree of its second argument has the same shape as a subexpression of its first argument. To solve the problem at hand, the second expression has to be a flat list of atomic expressions. Try Pick on the first example. You'll get a semi-explanatory error message. – Michael E2 May 04 '16 at 15:46
  • Ah, okay. Thanks for the tip! – dvitek May 04 '16 at 16:29
0

This is unlikely to be efficient, due to the Complement, and it's quite possible that the very first operation we do on the array "unpacks" it, but here's one possibility.

Let's suppose we have a sample array,

array = SparseArray@RandomInteger[{0, 5}, {10, 10}]; array // MatrixForm

enter image description here

We set a new background for this SparseArray by doing

array1 = SparseArray[array, Automatic, 1];

You can change 1 to whatever quantity you're looking for. This can include expressions like x^2 if those are the positions you are looking for.

Then, array1["NonzeroPositions"] returns a List of the non-background positions, which is everything but the positions of the 1s. We can then take the Complement of this list in the set of all possible position indices:

pos = Complement[Table[{i, j}, {i, 10}, {j, 10}]~Flatten~1, array1["NonzeroPositions"]]
(* {{1, 1}, {1, 2}, {1, 4}, {2, 3}, {2, 4}, {4, 10}, {5, 5}, {5, 9},
    {6, 2}, {6, 9}, {7, 3}, {7, 9}, {8, 6}, {8, 10}, {9, 7}, {10, 4}} *)

To verify:

Extract[array, array1["NonzeroPositions"]]
Extract[array, pos]
(* {4, 4, 3, 0, 5, 2, 0, 4, 2, 4, 0, 0, 5, 0, 3, 3, 5, 3,
    4, 2, 5, 0, 3, 5, 2, 3, 3, 4, 5, 4, 5, 4, 5, 2, 5, 4,
    5, 5, 5, 2, 3, 4, 5, 5, 0, 5, 0, 5, 3, 5, 2, 3, 5, 4,
    3, 4, 0, 2, 3, 4, 0, 5, 3, 4, 4, 4, 2, 3, 4, 4, 0, 3,
    0, 3, 5, 3, 4, 4, 2, 2, 3, 5, 4, 5} *)
(* {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1} *)
march
  • 23,399
  • 2
  • 44
  • 100
  • Hi March, this is roughly the same approach that ciao used in his comment. There the Unitize has a hope of being SparseArray-fast. Unfortunately for my purposes the array is far too big for Complement to be wieldy - this is the entire point of SparseArrays! – dvitek May 03 '16 at 16:31
  • @dvitek. Right, agreed, although in defense, I was working on this solution when ciao posted his comment. As far as the Complement goes, I feel like there's a way to quickly "invert" the list in this case, since it's so highly structured, but off the top of my head, I can't come up with a solution. – march May 03 '16 at 19:53
  • 1
    The changing the background is a cool trick, I'll admit :). It turns out that you can do lots of stuff that's SparseArray-fast - a simple Boole[#===1]& works for my purposes, and can also be used on expressions like x^2. – dvitek May 03 '16 at 20:04
  • For what it's worth, you can do the Complement far more quickly by walking through the list of indices and throwing anything in that you don't see - it's all highly structured so you're golden. But for my purposes (truly sparse arrays) it's quite slow, unfortunately. Thanks for your help! – dvitek May 03 '16 at 20:06