Question: Goal I'm working on a project that requires me visualize how far a given budget can reach on a map with different price zones. Example with 3 price zones, water, countryside, and city, each having different costs per meter traveled:  Red point is the origin, green point is just a sample point, it has no bearing on the result. Gradient goes from yellow (cheap) to purple (max budget) Issue Time complexity is currently really bad, causing time to grow very fast as I increase the resolution (decreasing size in meters of each cell). This is most likely due to the iterative creation of the scipy.sparse.lil_array when creating a array adjacency graph. Details My current processes to achieve this is by following 6 major steps (time added to show bottlenecks):  Create a Python list of shapely.geometry.Point which represents all points in a 2d grid which we will sample the prices from the geopandas geometry. The list is initialized with empty values before modifications. I iteratively set each index (one at a time) to a point created depending on the grid cell size. And lastly I convert the python list into a pandas dataFrame to later intersect with the geopandas regions. (13.6s)  I now use the dataFrame with points and convert it into a GeoDataFrame and the use gpd.sjoin against my diffrent price areas to receive in this case 3 geoDataFrame's each containing only the points which intersected with the given price area. (0.7s)  After mapping each point with their price, I want to combine it back into a 2d array with the value being the price to do this I use a Python 2d List and index and replace values using a conversion from point coordinates to index coordinates. (14.9s)  We now have a 2d array that contains the cost/weight for each cell. Now I want to convert it into a graph (matrix adjacency graph). To do this I use an iterative procces and inserts each edge into a scipy.sparse.lil_array and then after each element has been connected convert the array into a scipy.sparse.csr_array (136.7s)  I now call scipy.sparse.csgraph.dijkstra and use the distance array created. (0.4s)  Due to a fixed price needed for each node I once again iteratively loop through the array adding the fixed cost. (0.2s)   In addition to step 4 I should mention that I use more than 4 edges per node. Namely using a 5x5 centered on the node (but ignoring duplicate cardinal and diagonal paths), but this introduces a bit of a error margin as it is possible to avoid cost of a cell by jumping diagonal over it, but the error margin should be small enough that it can be ignored. This is due to the otherwise square result:  More info regarding step 4: To create the edges I iterate trough a list for each node adjecency_offsets= [ (-3, -2, 3.60), (-3, -1, 3.16), (-3, 1, 3.16), (-3, 2, 3.60), (-2, -3, 3.60), (-2, -1, 2.24), (-2, 1, 2.24), (-2, 3, 3.60), (-1, -3, 3.16), (-1, -2, 2.24), (-1, -1, 1.42), (-1, 0, 1 ), (-1, 1, 1.42), (-1, 2, 2.24), (-1, 3, 3.16), ( 0, -1, 1 ), ( 0, 1, 1 ), ( 1, -3, 3.16), ( 1, -2, 2.24), ( 1, -1, 1.42), ( 1, 0, 1 ), ( 1, 1, 1.42), ( 1, 2, 2.24), ( 1, 3, 3.16), ( 2, -3, 3.60), ( 2, -1, 2.24), ( 2, 1, 2.24), ( 2, 3, 3.60), (-3, -2, 3.60), (-3, -1, 3.16), (-3, 1, 3.16), (-3, 2, 3.60) ]  first integer is the row offset, 2nd is column offset, and 3rd value is the distance multiplier to reach the cell from the node (using approximations instead of sq(2) etc) Then by using following code we create a graph and populate the adjacency array def Step_4(arr: list[list[int]]): # arr is a 2d-array of cost for each point (the index) height = len(arr) width = len(arr[0]) num_nodes = width*height tmp_array = scipy.sparse.lil_array((num_nodes, num_nodes), dtype=np.float32) # we dont need more precision than this # For every row and column in 2d array for row in range(height): for col in range(width): # Create a src vertex src = row * width + col # and find all edges for (r_off, c_off, w_mul) in adjecency_offsets: # ignore invalid destinations if (row + r_off < 0 or col + c_off < 0): continue if (row + r_off >= height or col + c_off >= width): continue # valid destination dst = (row + r_off) * width + col + c_off # Multiply cost of cell with distance to cell float_w = w_mul * arr[row + r_off][col + c_off] # Set weight in adjacency array tmp_array[src, dst] = float_w # return crs version due to faster pathfinding return tmp_array.tocsr()  Question Is this a good approach to this problem? (even if it has blocky unnatural edges due to it using graph with limited amount of edges for each node) If so how can I more efficiently create my adjacency graph (Reduce time spent in step 4)? What I have tried  I have tried drawing lines from origin in a circle and saving how long the line can get before the budget runs out. this gave me poor runtime (probably a coding issue) but also can not handle finding a cheaper route by going further thru cheaper areas.  Using a list adjacency graph. whilst this reduced creation time by a lot, it also increased path finding time by a larger margin (needed to use a self implemented version of Dijkstra based on Wikipedia using a heap and early exit using budget) in total the total time of creating and path finding is around 2x slowdown compared to the adjacency matrix method mentioned above.  I have tried using DOK format instead but that resulted in a time of 674s instead of 120s. So not sure whats happening there.  

Answer: The complexity of this code seems fine. The main issue is that your code is a pure-Python code: it is not vectorized (i.e. it does not spent most of its time in fast native functions). As a result, it is very slow because Python codes are generally executed using the CPython interpreter. The key to fix this issue is to make this function native (e.g. using Cython, Numba, by rewriting this in languages like C/C++/Rust), or using modules like Numpy (e.g. using broadcasting operation). In this case, I choose to use Numba (so to keep a similar code and show how inefficient interpreted CPython code is). Numba is Python module able to compile specific Python functions (mainly operating on Numpy arrays) to native ones at runtime (Just-in-Time compiler a.k.a. JIT). Here is the optimized Numba code: import scipy import numpy as np @nb.njit('(int32[:,::1], int32[::1], int32[::1], float32[::1])') def compute(arr, row_offsets, col_offsets, weight_offsets): # arr is a 2d-array of cost for each point (the index) height = len(arr) width = len(arr[0]) num_nodes = width * height out_rows = np.zeros(num_nodes * weight_offsets.size, dtype=np.int32) out_cols = np.zeros(num_nodes * weight_offsets.size, dtype=np.int32) out_vals = np.zeros(num_nodes * weight_offsets.size, dtype=np.float32) nnz = 0 for row in range(height): for col in range(width): src = row * width + col for (r_off, c_off, w_mul) in zip(row_offsets, col_offsets, weight_offsets): if (row + r_off < 0 or col + c_off < 0) or (row + r_off >= height or col + c_off >= width): continue dst = (row + r_off) * width + col + c_off float_w = w_mul * arr[row + r_off][col + c_off] if float_w != 0: out_rows[nnz] = src out_cols[nnz] = dst out_vals[nnz] = float_w nnz += 1 return (out_rows, out_cols, out_vals, nnz) def Step_4_fast(list_of_list): arr = np.array(list_of_list, dtype=np.int32) row_offsets = np.array([item[0] for item in adjecency_offsets], dtype=np.int32) col_offsets = np.array([item[1] for item in adjecency_offsets], dtype=np.int32) weight_offsets = np.array([item[2] for item in adjecency_offsets], dtype=np.float32) out_rows, out_cols, out_vals, nnz = compute(arr, row_offsets, col_offsets, weight_offsets) return scipy.sparse.csr_array((out_vals[:nnz], (out_rows[:nnz], out_cols[:nnz])), shape=(arr.size, arr.size))  The main idea is to build a COO matrix manually (i.e. "ijv triplet") and then build a Scipy's CSR matrix from it. The COO matrix takes a significant space but be aware that CPython list of CPython list was already pretty inefficient in term of memory usage: each float item typically takes 24 bytes on mainstream machines and each float reference also takes 8 bytes which means 32 bytes per item instead of just 4 (so 8 times more than needed). I used 32-bit integer so to reduce the memory usage and improve performance. You should check if this is fine in your real-wold use-case though (i.e. check there are no overflow). Otherwise, you should just use 64-bit ones.  Benchmark Here are input for the benchmark: np.random.seed(42) # adjecency_offsets is like in the OP question arr = np.random.randint(1, 100_000, (500, 500)).tolist()  Here are performance results on my i5-9600KF CPU: Step_4(arr): 16.43 seconds Step_4_fast(arr): 0.16 seconds  Thus, the Numba code is about 100 times faster!  Notes One downside of the Numba code is that the function is compiled at runtime. It takes about 0.4 second on my machine. If this is a problem, then you can use Cython instead. Alternatively, you can cache the compilation so only the first compilation is slower (using the cache=True compilation flag). Note that only ~30% of the time is spent in the compute function showing how fast it is. ~60% of the time is spent in the COO to CSR matrix conversion. ~10% of the time is spent in the list-of-list to Numpy array conversion. Thus, for better performance, it would be a good idea to use a better data structure than a COO matrix (a kind of CSR one but with growing lists) and avoid list-of-list to Numpy array conversions (by not using lists at all but Numpy arrays instead). That being said, this is certainly not easy to do that efficiently in Numba (simpler in native languages). Another way to make the code faster is to parallelize the operations once the code is modified to directly operate on CSR-like matrices (array of growing list). However this is far from being easy in Python (even in Cython/Numba). In a native language, each thread can operate on its block of arr and modify a thread-local CSR-like matrix. The CSR-like matrices should then be merged if possible also in parallel (pretty difficult part). This last implementation should be an order of magnitude faster than the already much faster above Numba code on a mainstream CPU (e.g. with 6~8 cores).