Is it possible to accelerate cuSpatial with NVIDIA RT Cores

[pwrliang]

Hi cuSpatial community,

Our recent work, RayJoin(paper, code), leverages RT Cores to accelerate spatial join queries, specifically Line Segment Intersection (LSI) and Point in Polygon (PIP), by transforming spatial queries into ray tracing problems. This approach enables the queries to run on RT Cores, taking advantage of hardware acceleration. RayJoin demonstrates real-time performance, processing in milliseconds for datasets containing millions of polygons. We evaluated RayJoin against cuSpatial and other CPU/GPU spatial join solutions, noting significant performance gains. For instance, our PIP query completed in 227 milliseconds when joining two OpenStreetMap datasets with 1.6M and 303K polygons respectively on an RTX 3090, while cuSpatial took about 2120 seconds. Additionally, we implemented the LSI query on RT Cores, similar to cuSpatial’s Pairwise Linestring Intersections, but without the restriction to paired inputs, allowing it to identify arbitrary intersections between two sets of line segments.

Here are my thoughts. I believe that GIS practitioners, who are typical users of cuSpatial and often use Jupyter Notebook to interact with spatial data. A high-performance spatial join processing would enable quicker exploration of spatial data and expedite users' workflow. Furthermore, as both cuSpatial and RT Cores are NVIDIA products, integrating RT Cores-accelerated spatial queries into cuSpatial could be cool.

[jarmak-nv]

Hi @pwrliang - thanks for checking us out!

Your approach looks and seems very interesting/fast - really love seeing RT cores getting used for more things out in the wild.

cuSpatial's perf looks pretty rough in the paper, but I'm also not too shocked. It doesn't run at top performance on the GeForce cards since we require FP64 precision. The RTX3090 has ~35.6 FP32 TFLOPS vs ~0.55 FP64 TFLOPS.

This could be an interesting way to increase our perf on but we'd need to do exploration to see how we could integrate it with the existing GeoArrow data structures and FP64 restrictions. We do have two issues (1, 2) open on looking into how to support FP32 but currently it's not on the horizon.

We're always happy to get feature reqs or PRs if you'd like to formalize the request!

[pwrliang]

Hi Ben, thanks for your reply. We noticed that the consumer-grade GPUs have limited FP64 ALUs. Our method works around this limitation by scaling FP64 to INT128. The intersection tests are conducted with INT128 operations, so we eliminate the expensive FP64 ops and achieve high precision. We will evaluate cuSpatial on A100 and see how the performance improves. We will also investigate the solution for integrating our RT-solution under cuSpatial's codebase.

[trxcllnt]

@pwrliang Do you have the code for the LKEU ⊲⊳ PKEU benchmark that OOM'd in cuSpatial? I want to check whether #1381 fixes the issue.

I tried reconstructing the datasets from the lakes/parks datasets link in the RayJoin readme and running misc/cdb2shp.py to construct shapefiles, but I'm not sure how to transform the resulting GeoDataFrames into the point and polygon inputs required by the quadtree.

Thanks!

[trxcllnt]

Yes, there were three bugs PR #1381 addresses:

1. A counting/overflow error due to accumulating into uint32_t, fixed by accumulating into uint64_t (fixed here)
2. That counting bug was masking a bug in how we were calling Thrust. thrust::count_if has a limit on the number of input elements it can process (INT32_MAX), and thrust::copy_if's limit is I think UINT32_MAX. Since we always used thrust::copy_if and never processed more than UINT32_MAX elements, we never hit this limitation. The solution is to count/copy in chunks of no larger than INT32_MAX (noted here).
3. The routine could OOM due to attempting to allocate output buffers larger than available GPU memory. We have to over-allocate (assuming every poly will intersect with every point in the worst case), even though in most cases the intersection count will not need all that space. This over-allocation can easily OOM. The solution is to first try the over-allocation approach (because it's fastest), then fall back to computing an exact count for the output buffers (fixed here).

[trxcllnt]

@pwrliang is there any reason you set max_depth = 7, max_size = 50 in your test?

[pwrliang]

These values are just trying to make sure the code I provided works. In the evaluation of the paper, I set max_depth = 7 and max_size = # of points / 4^7 / 4 as suggested by https://docs.rapids.ai/api/cuspatial/legacy/user_guide/cuspatial_api_examples/

[trxcllnt]

I can confirm PR #1381 fixes the OOM issue.

Starting from your script, I added the ability to serialize the transformed structures to Arrow tables (so rerunning is fast), compute an optimum max_depth and max_size over a user-provided threshold, and compare both single and double-precision performance.

These are the results on my 4090:

$ python rayjoin_pip.py 10
dtype: float32
  params: {'max_depth': 10, 'max_size': 18.610915184020996, 'scale': 0.0977388471364975}
  tree build (83.222ms)
  compute bboxes (34.97ms)
  num poly/quad pairs: 611328 (5.442ms)
  num intersections: 461937 (2928.46ms)
  point in polygon (2968.92ms)
  total time: 3052.172ms

dtype: float64
  params: {'max_depth': 10, 'max_size': 18.610915184020996, 'scale': 0.0977388490235464}
  tree build (73.656ms)
  compute bboxes (36.878ms)
  num poly/quad pairs: 611294 (5.507ms)
  num intersections: 464331 (17671.38ms)
  point in polygon (17713.809ms)
  total time: 17787.48ms

Passing a higher threshold yields a lower max_depth/higher max_size, and performance is impacted due to testing more points per quadrant against each poly. So using the largest max_depth and smallest max_size is best for performance.

[pwrliang]

Thank you for the bug fix and your benchmark scripts. I will reevaluate cuSpatial after your fix is merged into the upcoming release. This fix enables me to run PIP queries in one shot. Without this fix, I have to execute the queries in batch mode as a workaround, resulting in significantly poorer performance. Unfortunately, our paper is already camera-ready and we cannot update the performance numbers anymore, but I will post the updated performance numbers on the RayJoin GitHub page.