Skip to main content

sparse arrays - SparseArray row operations



Given a large (very) sparse matrix, A, how can I efficiently "operate" on only the nonzeros in a given row? For example: For each row in A, generate a list of column indices that have a magnitude (absolute value) greater than a threshold, r.


My current approach is similar in spirit to a previous answer by Leonid (Efficient way to combine SparseArray objects?). I shamelessly steal his implementation of spart.


My current approach


First, let's build a random SparseArray for testing purposes.


randomSparseArray[n_Integer, r_Integer] := SparseArray[
Table[{Random[Integer, {1, n}], Random[Integer, {1, n}]} ->RandomReal[],{r}]
,{n, n}];

A = randomSparseArray[10, 20]; (*10x10 matrix with 20 nonzeros*)


When Part is used to grab a row of A, a SparseArray representation is returned. Okay so let's rip this representation apart into its raw form using spart. For example the 3rd row of A


HoldPattern[spart[SparseArray[s___], p_]] := {s}[[p]];    
raw = spart[A[[3]], 4];
(*{1, {{0, 2}, {{10}, {5}}}, {0.534313, 0.981372}}*)

And then use a combination of Part and Position to extract the column indices from raw . (here I use a threshold of 0.5)


ind = Flatten@Position[Flatten@raw[[3]],x_/;Abs[x]>0.5];
Flatten[raw[[2, 2]]][[ind]]

Wrapping this all up into a function:



thresholdIndex[A_SparseArray, r_] := Module[{raw, ind},
raw = spart[A, 4];
ind = Flatten@Position[Flatten@raw[[3]], x_ /; Abs[x] > r];
Flatten[raw[[2, 2]]][[ind]]
]

We can now use Map to hit each row of A and search for column indices for elements with magnitudes greater than say, 0.5


Map[thresholdIndex[#, .5] &, A]

This whole process feels a little roundabout. First Mathematica has to extract a new SparseArray representation of a row (is this expensive?), we then chop it into pieces in order to work on it. Is there an easier way to do this while still maintaining performance?



My other idea is to apply spart to the original matrix and work with the CSR (compressed-sparse-row) representation. But this seems to defeat the purpose of even using a SparseArray from the get go.


Sidebar


I've been starting to implement an efficient algebraic multigrid (AMG) package in Mathematica. AMG is a fast iterative method traditionally used to solve large sparse matrix equations produced from the discretization of elliptic PDEs. One of the steps in the algorithm is very similar in flavor to thresholdIndex. As a sidebar, does anyone know if a nice AMG has been implemented in Mathematica? I would like to build a tool similar to PyAMG (http://code.google.com/p/pyamg/).




Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

list manipulation - Selecting multiple columns from a matrix?

Sample data: data = { {{2013, 1, 1}, 24.13, 167.67, 231.82}, {{2013, 1, 2}, 32.15, 170.92, 225.99}, {{2013, 1, 3}, 35.43, 172.68, 221.67}, {{2013, 1, 4}, 36.73, 173.05, 218.32}, {{2013, 1, 5}, 58.19, 165.96, 197.05}, {{2013, 1, 6}, 69.99, 163.50, 187.52}, {{2013, 1, 7}, 71.37, 154.21, 175.58}, {{2013, 1, 8}, 72.51, 149.66, 163.25}}; I want a DateListPlot with three graphs, so for a matrix formed by columns 1 and 2, one for columns 1 and 3, and 1 for columns 1 and 4. At the moment I'm using this code: data2 = Transpose[{data[[All, 1]], data[[All, 2]]}]; data3 = Transpose[{data[[All, 1]], data[[All, 3]]}]; data4 = Transpose[{data[[All, 1]], data[[All, 4]]}]; DateListPlot[{data2, data3, data4}, Joined -> True, Filling -> {3 -> {1}}] but I have a hunch that this can be done more efficiently. I don't like the Transpose s in particular. Any ideas? edit (for extra credit) What if I need to multiply the second column by 2, which in my solution is simp...

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...