今年も早いもので気がついたら Advent Calendar の季節ですが、今回のこちらの記事は RDBMS-GIS(MySQL,PostgreSQLなど) Advent Calendar 2020 - Qiita の 12/5(土) のものです。
2020年はこの辺からスタートして、かなり性能的に面白いところまで持ってくることができたので、年末の締めにGPU版PostGISとGiSTインデックス対応についてご紹介させていただこうかと。
【緩募】PostGISのGPU対応に興味がある方
— 海外 浩平｜KaiGai Kohei🌻 (@kkaigai) 2020年1月29日
来年度のR&Dのターゲットとして地理情報処理の高速化を考えているのですが、当方だけではGIS系のデータやワークロードを保有していないため、これらを提供できる方、うまく問題を改善できれば実システムに適用できるぜという方を探しています。
何をやりたいか？
元々、PG-StromではIoT/M2M領域でありがちな『時間経過とともに溜まっていくログデータの処理、分析』にフォーカスを当てていたが、携帯電話や自動車といった移動体デバイスの生成するログを考えると、その中にGPSの生成する『緯度、経度』情報が含まれる事もしばしばである。
ホンマに商売になるのか分からないところではあったが、こういったデータをPostGISを使って分析したい、という問合せはちょくちょくもらっていたので、相応に需要はあるのだろうという事で、春先からGPU版PostGIS関数の実装を進めていた。
ターゲットとするワークロードは以下のようなもので、携帯電話や自動車から収集したログデータ（これには、典型的にはデバイスID、タイムスタンプ、緯度・経度、その他の属性が含まれる）と、例えばプッシュ広告を配信する地域や、事故・渋滞情報を配信するエリアを定義するエリア情報を突合するという処理を想定している。
つまり、単純化すると多角形領域の中に含まれる点を抽出するという処理になる。
これをナイーブに実装するだけであれば、st_containsやst_dwithinといった関数をGPUで実装するだけで*1よいのだが、例えば10万個の多角形領域 × 1000万個の座標データの組み合わせを検査するとなればその組み合わせは1兆通りにのぼり、１枚あたり数千コアを搭載するGPUといえども、なかなかにしんどい処理となる。
こういったケースで絞り込みを効率的に行うため、PostgreSQL/PostGISにはGiSTインデックスという仕組みが存在する。
GiSTインデックス（R木）
位置データの場合、緯度・経度という二次元的な広がりがあるため、ある値が別の値より大きいのか、小さいのかという関係を定義する事ができない。
そのため、値の大小関係を利用して目的の要素へと直線的に木構造を降下していくB木のようにはいかないが、それでも効率的に絞り込みを行うためR木と呼ばれるインデックス構造を作る事ができる。
R木インデックスでは、多角形要素（Polygon）はどれだけ複雑な形状を持っていようとも、その形状を完全に包含する長方形（Bounding-Box）であるとして扱われる。ある点要素がその長方形の中に含まれていれば多角形要素の中に含まれる可能性があるが、その長方形の外側であれば、複雑な多角形要素との当たり判定を行うまでもなく『包含されていない』と結果を返す事ができる。
以下の図のケースでは、R木インデックスのRoot要素であるR1はその子要素であるR3、R4、R5を完全に包含し、R2は同じくその子要素であるR6、R7を完全に包含する長方形として定義されてる。
もし以下の図の位置が検索キーとして与えられた場合、R1およびR2は共に検索キーを包含し、その子要素をチェックする事となる。
次に、R1の子要素R3、R4、R5と検索キーの重なりを検査する。すると、R4は検索キーを包含するが、R3とR5は検索キーを包含しない。
したがって、R3の子要素R8、R9、R10およびR5の子要素R13、R14はそもそもチェックするまでもなく、検索キーを包含しないという事になる。
次に、R4の子要素R11とR12と検索キーの重なりを検査する。すると、R12は検索キーを包含するが、R11は包含しない。
この階層はR木のLeaf要素であるのでこれ以上探索が深くなることはなく、R12がインデックスする多角形要素が本当に検索キーを包含するのかどうか、st_contains関数などを用いてチェックする。
一方、Root要素で検索キーを包含していたR2の子要素R6、R7は共に検索キーを包含しないため、ここでインデックスの探索は打ち切りとなる。
この例では単純化されているが、実際にはR木の一階層を降下する毎に100要素程度の『当たり判定』をシーケンシャルに行っていくため、実はインデックス探索とはいえ、そこそこ計算ヘビーな処理にはなってしまう。
一方でGPU向きの性質としては、突合処理中のR木インデックスはRead-Onlyなデータ構造であるため並列度を上げやすく、その気になれば数千コアを同時に稼働してのインデックス探索を行う事ができる。
簡単なベンチマーク
PG-Stromに実装したGPU版PostGISの機能には、多角形や点などジオメトリ要素間の包含関係を判定するst_contains関数も含まれており、これを用いてのエリア定義情報×座標情報の突合処理のベンチマークを行ってみる事にする。
サンプルとして使用するのは、国土地理院の公開している市区町村の形状データと、日本列島を概ね包含する領域にランダムに打った点との突合処理。
この中で、東京都に属する市区町村に含まれるエリアに打たれた点の数を、市区町村ごとにカウントするものとする。
実行するクエリは以下の通り。geo_japanテーブルには市町村の形状データを、fgeopointテーブルにはランダムに生成した点のデータを挿入する。
SELECT n03_001,n03_004,count(*) FROM geo_japan j, fgeopoint p WHERE st_contains(j.geom, st_makepoint(x,y)) AND j.n03_001 like '東京都’ GROUP BY n03_001,n03_004;
ランダムなデータの生成は以下の通り。
=# CREATE TABLE
geopoint (
gid int primary key,
x float8,
y float8);
CREATE TABLE
=# INSERT INTO geopoint (SELECT x, pgstrom.random_float(0, 123.0, 154.2),
pgstrom.random_float(0, 20.0, 46.2)
FROM generate_series(1,10000000) x);
INSERT 0 10000000
postgres=# SELECT * FROM geopoint LIMIT 4;
gid | x | y
-----+--------------------+--------------------
1 | 133.5737876430963 | 23.438477765972948
2 | 133.47253950874904 | 22.512966607004856
3 | 136.9879882356096 | 22.22637613640464
4 | 126.3652188637132 | 27.186177021165463
(4 rows)GPU側はこれと同じ内容をGPUメモリストアに挿入する。
GPUメモリストアに関しては、後日、PostgreSQL Advent Calendar 2020の16日目でも記載する予定ですので、そちらをご覧あれ。
=# CREATE FOREIGN TABLE
fgeopoint (
gid int,
x float,
y float)
SERVER gstore_fdw
OPTIONS (gpu_device_id '0',
base_file '/opt/nvme/fgeopoint.base',
redo_log_file '/opt/pmem/fgeopoint.redo',
max_num_rows '12000000',
primary_key 'gid');
CREATE FOREIGN TABLE
=# INSERT INTO fgeopoint (SELECT * FROM geopoint);
INSERT 0 10000000これで準備が整った。さぁ、実行してみる事にしよう。
先ずはCPUのみのケース。予め並列度を上げておく。
=# set pg_strom.enabled = off;
SET
=# set max_parallel_workers_per_gather = 99;
SET
=# SELECT n03_001,n03_004,count(*)
FROM geo_japan j, geopoint p
WHERE st_contains(j.geom, st_makepoint(x,y))
AND j.n03_001 like '東京都'
GROUP BY n03_001,n03_004;
n03_001 | n03_004 | count
---------+------------+-------
東京都 | あきる野市 | 90
東京都 | 三宅村 | 53
東京都 | 三鷹市 | 14
： ：
東京都 | 青ヶ島村 | 4
東京都 | 青梅市 | 109
(63 rows)
Time: 30539.097 ms (00:30.539)約30秒と言ったところ。
なおEXPLAIN ANALYZEを見てみると、応答時間30.7秒のうち30.68秒をNested Loop+Index Scanの箇所で消費しており、ここがワークロードの肝である事がわかる。
postgres=# EXPLAIN ANALYZE
SELECT n03_001,n03_004,count(*)
FROM geo_japan j, geopoint p
WHERE st_contains(j.geom, st_makepoint(x,y))
AND j.n03_001 like '東京都'
GROUP BY n03_001,n03_004;
QUERY PLAN
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Finalize GroupAggregate (cost=7483675086.12..7483829745.31 rows=4858 width=29) (actual time=30709.855..30710.080 rows=63 loops=1)
Group Key: j.n03_001, j.n03_004
-> Gather Merge (cost=7483675086.12..7483829550.99 rows=19432 width=29) (actual time=30709.838..30732.270 rows=244 loops=1)
Workers Planned: 4
Workers Launched: 3
-> Partial GroupAggregate (cost=7483674086.06..7483826236.39 rows=4858 width=29) (actual time=30687.466..30687.572 rows=61 loops=4)
Group Key: j.n03_001, j.n03_004
-> Sort (cost=7483674086.06..7483712111.50 rows=15210175 width=21) (actual time=30687.452..30687.475 rows=638 loops=4)
Sort Key: j.n03_001, j.n03_004
Sort Method: quicksort Memory: 73kB
Worker 0: Sort Method: quicksort Memory: 74kB
Worker 1: Sort Method: quicksort Memory: 74kB
Worker 2: Sort Method: quicksort Memory: 76kB
-> Nested Loop (cost=0.41..7481859623.72 rows=15210175 width=21) (actual time=71.496..30686.278 rows=638 loops=4)
-> Parallel Seq Scan on geopoint p (cost=0.00..88695.29 rows=2500029 width=16) (actual time=0.012..207.553 rows=2500000 loops=4)
-> Index Scan using geo_japan_geom_idx on geo_japan j (cost=0.41..2992.66 rows=1 width=1868) (actual time=0.012..0.012 rows=0 loops=10000000)
Index Cond: (geom ~ st_makepoint(p.x, p.y))
Filter: (((n03_001)::text ~~ '東京都'::text) AND st_contains(geom, st_makepoint(p.x, p.y)))
Rows Removed by Filter: 0
Planning Time: 0.156 ms
Execution Time: 30732.422 ms
(21 rows)次に、GPUによるGiSTインデックスの探索を試してみる。
使用したのは、先週届いたばかりの最新鋭モデル、NVIDIA A100。
こやつを利用して実行した結果がコレ。
=# SELECT n03_001,n03_004,count(*)
FROM geo_japan j, fgeopoint p
WHERE st_contains(j.geom, st_makepoint(x,y))
AND j.n03_001 like '東京都'
GROUP BY n03_001,n03_004;
n03_001 | n03_004 | count
---------+------------+-------
東京都 | あきる野市 | 90
東京都 | 三宅村 | 53
： ：
東京都 | 青ヶ島村 | 4
東京都 | 青梅市 | 109
(63 rows)
Time: 316.673 msどやっ！！316msで100倍近い高速化！
EXPLAIN ANALYZEを覗いてみると、GpuJoinのdepth=1でGpuGiSTJoinが選択され、13.28MBのGiSTインデックスをGPUへロードして結合処理を行っている事がわかる。
=# EXPLAIN ANALYZE
SELECT n03_001,n03_004,count(*)
FROM geo_japan j, fgeopoint p
WHERE st_contains(j.geom, st_makepoint(x,y))
AND j.n03_001 like '東京都'
GROUP BY n03_001,n03_004;
QUERY PLAN
------------------------------------------------------------------------------------------------------------------------------------------
GroupAggregate (cost=6141933.29..6142042.59 rows=4858 width=29) (actual time=329.118..329.139 rows=63 loops=1)
Group Key: j.n03_001, j.n03_004
-> Sort (cost=6141933.29..6141945.43 rows=4858 width=29) (actual time=329.107..329.110 rows=63 loops=1)
Sort Key: j.n03_001, j.n03_004
Sort Method: quicksort Memory: 29kB
-> Custom Scan (GpuPreAgg) (cost=6141575.10..6141635.83 rows=4858 width=29) (actual time=328.902..328.911 rows=63 loops=1)
Reduction: Local
Combined GpuJoin: enabled
-> Custom Scan (GpuJoin) on fgeopoint p (cost=3759781.06..6712457.81 rows=60840000 width=21) (never executed)
Outer Scan: fgeopoint p (cost=0.00..100000.00 rows=10000000 width=16) (never executed)
Depth 1: GpuGiSTJoin(plan nrows: 10000000...60840000, actual nrows: 10000000...2553)
HeapSize: 7841.91KB (estimated: 3113.70KB), IndexSize: 13.28MB
IndexFilter: (j.geom ~ st_makepoint(p.x, p.y)) on geo_japan_geom_idx
Rows Fetched by Index: 4952
JoinQuals: st_contains(j.geom, st_makepoint(p.x, p.y))
-> Seq Scan on geo_japan j (cost=0.00..8928.24 rows=6084 width=1868) (actual time=0.164..17.723 rows=6173 loops=1)
Filter: ((n03_001)::text ~~ '東京都'::text)
Rows Removed by Filter: 112726
Planning Time: 0.344 ms
Execution Time: 340.415 ms
(20 rows)
*1:いや、そこそこ大変ではあるが…。
