Skip to main content

functions - Why is the new PositionIndex horribly slow?


This issue has largely been mitigated in 10.0.1. New timings for the final test below are:


Needs["GeneralUtilities`"]
a = RandomInteger[9, 5*^5];
myPosIdx[a] // AccurateTiming
cleanPosIdx[a] // AccurateTiming (* see self-answer below *)

PositionIndex[a] // AccurateTiming


0.0149384

0.0149554

0.0545865

Still several times slower here than the readily available alternatives but no longer devastating.





Disconcertingly I have discovered that the new (v10) PositionIndex is horribly slow.


Using Szabolcs's clever GatherBy inversion we can implement our own function for comparison:


myPosIdx[x_] :=
<|Thread[x[[ #[[All, 1]] ]] -> #]|> & @ GatherBy[Range @ Length @ x, x[[#]] &]

Check that its output matches:


RandomChoice[{"a", "b", "c"}, 50];

myPosIdx[%] === PositionIndex[%]



True

Check performance in version 10.0.0 under Windows:


a = RandomInteger[99999, 5*^5];
myPosIdx[a] // Timing // First
PositionIndex[a] // Timing // First



0.140401

0.920406

Not a good start for the System` function, is it? It gets worse:


a = RandomInteger[999, 5*^5];
myPosIdx[a] // Timing // First
PositionIndex[a] // Timing // First



0.031200

2.230814

With fewer unique elements PositionIndex actually gets slower! Does the trend continue?


a = RandomInteger[99, 5*^5];
myPosIdx[a] // Timing // First
PositionIndex[a] // Timing // First



0.015600

15.958902

Somewhere someone should be doing a face-palm right about now. Just how bad does it get?


a = RandomInteger[9, 5*^5];
myPosIdx[a] // Timing // First
PositionIndex[a] // Timing // First



0.015600

157.295808

Ouch. This has to be a new record for poor computational complexity in a System function. :o



Answer



First let me note that I didn't write PositionIndex, so I can't speak to its internals without doing a bit of digging (which at the moment I do not have time to do).


I agree performance could be improved in the case where there are many collisions. Let's quantify how bad the situation is, especially since complexity was mentioned!


We'll use the benchmarking tool in GeneralUtilities to plot time as a function of the size of the list:


Needs["GeneralUtilities`"]

myPosIdx[x_] := <|Thread[x[[#[[All, 1]]]] -> #]|> &@
GatherBy[Range@Length@x, x[[#]] &];
BenchmarkPlot[{PositionIndex, myPosIdx}, RandomInteger[100, #] &, 16, "IncludeFits" -> True]

which gives:


PositionIndex benchmark


While PositionIndex wins for small lists (< 100 elements), it is substantially slower for large lists. It does still appear to be O(nlogn), at least.


Let's choose a much larger random integer (1000000), so that we don't have any collisions:


enter image description here


Things are much better here. We can see that collisions are the main culprit.



Now lets see how the speed for a fixed-size list depends on the number of unique elements:


BenchmarkPlot[{PositionIndex, myPosIdx}, RandomInteger[#, 10^4] &, 
2^{3, 4, 5, 6, 7, 8, 9, 10, 11, 12}]

enter image description here


Indeed, we can see that PositionIndex (roughly) gets faster as there are more and more unique elements, whereas myPosIdx gets slower. That makes sense, because PositionIndex is probably appending elements to each value in the association, and the fewer collisions the fewer (slow) appends will happen. Whereas myPosIdx is being bottlenecked by the cost of creating each equivalence class (which PositionIndex would no doubt be too, if it were faster). But this is all academic: PositionIndex should be strictly faster than myPosIdx, it is written in C.


We will fix this.


Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]