読者です 読者をやめる 読者になる 読者になる

進捗)SSD-to-GPU ダイレクトSQL実行機能

ここ暫くブログでまとめていなかった、SSD-to-GPUダイレクトSQL実行機能の進捗について。

この機能をかいつまんで言うと、NVMe-SSDに格納されているPostgreSQLのデータブロックをGPU RAMに直接転送し、そこでSQLのWHERE句/JOIN/GROUP BYを実行することで見かけ上のI/O量を削減するという代物である。
NVIDIAのTesla/Quadro GPUが対応するGPUDirect RDMA機能を使い、SSD<=>GPU間のデータ転送を仲介するLinux kernel moduleを使えば、CPU/RAMにデータをロードする前にGPU上での処理を行うことができる。

しばらく前からScan系の処理には対応していたが、JOIN/GROUP BYへの対応を加え、さらにPostgreSQL v9.6のCPU並列にも追従したということで、簡単なベンチマークなら取れるくらいまで開発は進んできている。

という事で、現時点での実力がどの程度なのか、手元の環境を使って測定してみることにした。
使用したのはDELL R730にTesla K80*1Intel SSD 750(400GB)を2枚搭載したモデル。
GPUDirect RDMAの制約から、SSDGPUが同一のCPU又はPCIeスイッチに接続されている必要があるため、共にCPU2配下のPCIeスロットに接続している。

ベンチマークに使用したのはTPC-Hを簡略化したStar Schema Benchmark(SSBM)と呼ばれるテストで、中核となるlineorderテーブルのサイズが明らかにRAMに収まり切らないようscaling factorを調整している。

SSBMでは何種類かクエリが定義されているが、基本的にはWHERE句/JOIN/GROUP BYという、集計クエリの典型的なものである。例えば以下のようなものである。

SELECT sum(lo_revenue), d_year, p_brand1
  FROM lineorder, date1, part, supplier
 WHERE lo_orderdate = d_datekey
   AND lo_partkey = p_partkey
   AND lo_suppkey = s_suppkey
   AND p_category = 'MFGR#12‘
   AND s_region = 'AMERICA‘
 GROUP BY d_year, p_brand1
 ORDER BY d_year, p_brand1;

それを実行してみた結果が以下の通り。じゃん!

青色がPostgreSQL(ファイルシステム経由I/O)での実行結果、橙色がPG-Stromでの実行結果で、それぞれ色の濃い方がSSD x1枚での実行結果、色の薄い方がSSD x2枚(md-raid0)での実行結果。

これを見る限り、SSD x1枚の場合、PG-Stromは概ね理論限界*2に達しているが、SSD x2枚のスループットである4.4GB/sにはまだ届いていない。ただ、これはGPUがTesla K80というやや古いモデルであった事や、GROUP BYでの集約演算がNumeric型の総和や平均など、Kepler世代のGPUには比較的辛い処理だという事は考慮しなければいけないだろう。(例えばQ4-1やQ4-2ではI/Oよりも集計処理にボトルネックがあるように見える)
一方PostgreSQLの場合はSSD x1枚で600-750MB/s程度、SSD x2枚で1.6GB/s程度のスループットなので、PG-Stromの場合はI/Oネックな集計処理において2.3~3.0倍程度の優位性があるという事になるだろう。

この辺、SSDGPU、あるいはCPUのどこに律速要因があるのか追及するのはなかなか骨の折れる作業ではあるが、近々、最新モデルであるTesla P40で実行できるようになるので、最新世代のGPUではどの程度のスループットまで耐えられるのか、試してみたいところである。

本日の記事のより詳しい内容は、先日、BigData基盤研究会#7の方で喋らせて頂いた資料を公開しているので、そちらも併せて参照していただければ。

www.slideshare.net

*1:ちょっと古い

*2:Intel SSD 750(400GB)のSeqReadスペックは2200MB/s

PCIeスロット接続型NVMe-SSDまとめ(2017年4月時点)

PCIeスロット接続型のNVMe-SSDのスペック等、各社どうだったけな~と探すものの、案外まとまった情報がないので自分でまとめてみた。ブックマーク代わりです。
基本的に各社のWebに掲載されているカタログスペックを転記。手作業なので内容の正確さに関しては保証しかねます。(指摘頂ければ直します)

Intel

コンシューマ向け

モデル Intel SSD 750
400GB(MLC)
Intel SSD 750
800GB(MLC)
Intel SSD 750
1.2TB(MLC)
形状 HHHL HHHL HHHL
PCIe PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4
容量 400GB 800GB 1.2TB
SeqRead 2100MB/s 2100MB/s 2400MB/s
SeqWrite 800MB/s 800MB/s 1200MB/s
RandRead 420kIOPS 420kIOPS 440kIOPS
RandWrite 210kIOPS 210IOPS 290kIOPS
リリース Q3'15 Q3'15 Q2'15

エンタープライズ向け

モデルl DC P3608
1.6TB(MLC)
DC P3608
3.2TB(MLC)
DC P3608
4.0TB(MLC)
DC P3700
400GB(MLC)
DC P3700
800GB(MLC)
DC P3700
1.6TB(MLC)
DC P3700
2.0TB(MLC)
形状 HHHL HHHL HHHL HHHL HHHL HHHL HHHL
PCIe PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4
容量 1.6TB 3.2TB 4.0TB 400GB 800GB 1.6TB 2.0TB
SeqRead 5000MB/s 4500MB/s 5000MB.s 2700MB/s 2800MB/s 2800MB/s 2800MB/s
SeqWrite 2000MB/s 2600MB/s 3000MB/s 1080MB/s 1900MB/s 1900MB/s 1900MB/s
RandRead 850kIOPS 850kIOPS 850kIOPS 450kIOPS 460kIOPS 450kIOPS 450kIOPS
RandWrite 150kIOPS 80kIOPS 50kIOPS 75kIOPS 90kIOPS 150kIOPS 175kIOPS
リリース Q3'15 Q3'15 Q3'15 Q2'14 Q2'14 Q2'14 Q2'14

※ Optane DCモデルに関してはARKにデータが上がってから更新予定。

Samsung

エンタープライズ向け

モデル PM1725a
1.6TB(V-NAND)
PM1725a
3.2TB(V-NAND)
PM1725a
6.4TB(V-NAND)
PM1725
3.2TB(V-NAND)
PM1725
6.4TB(V-NAND)
形状 HHHL HHHL HHHL HHHL HHHL
PCIe PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x8
容量 1.6TB 3.2TB 6.4TB 3.2TB 6.4TB
SeqRead 5840MB/s 6200MB/s 6300MB/s 5600MB/s 5600MB/s
SeqWrite 2100MB/s 2600MB/s 2600MB/s 1800MB/s 1800MB/s
RandRead 1000kIOPS 1000kIOPS 1000kIOPS 1000kIOPS 1000lIOPS
RandWrite 140kIOPS 180kIOPS 180kIOPS 130kIOPS 130kIOPS
リリース Q3'16 Q3'16 Q3'16 Q3'15 Q3'15

HGST

エンタープライズ向け

モデル Ultrastar
SN260(1.6TB)
Ultrastar
SN260(3.2TB)
Ultrastar
SN260(6.4TB)
Ultrastar
SN150(1.6TB)
Ultrastar
SN150(3.2TB)
形状 HHHL HHHL HHHL HHHL HHHL
PCIe PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x4 PCIe 3.0 x4
容量 1.6TB 3.2TB 6.4TB 1.6TB 3.2TB
SeqRead 6100MB/s 6100MB/s 6100MB/s 3000MB/s 3000MB/s
SeqWrite 2200MB/s 2200MB/s 2200MB/s 1600MB/s 1600MB/s
RandRead 1200kIOPS 1200kIOPS 1200kIOPS 743kIPOS 743kIPOS
RandWrite 200kIOPS 200kIOPS 200kIOPS 140kIPS 140kIOPS
リリース Q4'16 Q4'16 Q4'16 Q2'15 Q2'15

Seagate

エンタープライズ向け

モデル Nytro XP7200
7.7TB(MLC)
Nytro XP6500
1.5TB
Nytro XP6500
4.0TB
Nytro XP7102
800GB
Nytro XP7102
1.6TB
形状 FHHL FHHL FHHL HHHL HHHL
PCIe PCIe 3.0 x16 PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x4 PCIe 3.0 x4
容量 7.7TB 1.3TB 3.4TB 800GB 1.6TB
SeqRead 10000MB/s 4.0GB/s 4.0GB/s 2500MB/s 2500MB/s
SeqWrite 2300MB/s 1.5GB/s 2.0GB/s 850MB/s 900MB/s
RandRead 950kIOPS 300kIPOS 275kIOPS 245kIOPS 245kIOPS
RandWrite 64kIPS 100kIOPS 75kIOPS 35kIOPS 40kIOPS
リリース Q3'16 ??? ??? Q3'16 Q3'16

ショックだ・・・Nytro XP7200の3.8TB版がサイトから消えている・・・。 orz

OCZ社

コンシューマ向け

モデル OCZ RD400a
128GB(MLC)
OCZ RD400a
256GB(MLC)
OCZ RD400a
512GB(MLC)
OCZ RD400a
1.0TB(MLC)
形状 HHHL HHHL HHHL HHHL
PCIe PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4
容量 128GB 256GB 512GB 1.0TB
SeqRead 2200MB/s 2600MB/s 2600MB/s 2600MB/s
SeqWrite 620MB/s 1150MB/s 1600MB/s 1550MB/s
RandRead 170kIPS 210kIOPS 190kIOPS 210kIOPS
RandWrite 110kIOPS 140kIPS 120kIOPS 130kIOPS
リリース Q2'16 Q2'16 Q2'16 Q2'16

Plextor

コンシューマ向け

モデル M8Pe(Y)
128GB
M8Pe(Y)
256GB
M8Pe(Y)
512GB
M8Pe(Y)
1.0TB
形状 HHHL HHHL HHHL HHHL
PCIe PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4 PCIe 3.0 x4
容量 128GB 256GB 512GB 1.0TB
SeqRead 1600MB/s 2000MB/s 2300MB/s 2500MB/s
SeqWrite 500MB/s 900MB/s 1300MB/s 1400MB/s
RandRead 120kIPS 210kIOPS 260kIPS 280kIOPS
RandWrite 130kIOPS 230kIPS 250kIOPS 240kIPS
リリース Q4'16 Q4'16 Q4'16 Q4'16

Kingstone社

エンタープライズ向け

モデル DCP1000
800GB
DCP1000
1.6TB
DCP1000
3.2TB
形状 HHHL HHHL HHHL
PCIe PCIe 3.0 x8 PCIe 3.0 x8 PCIe 3.0 x8
容量 800GB 1.6TB 3.2TB
SeqRead 6800MB/s 6800MB/s 6800MB/s
SeqWrite 5000MB/s 6000MB/s 6000MB/s
RandRead 900kIOPS 1100kIPS 1000kIOPS
RandWrite 145kIOPS 200kIOPS 180kIPS
リリース Q1'17 Q1'17 Q1'17

※なぜかメーカーのサイトから消えている・・・。

AWSのP2.*インスタンスで PG-Strom を試す

従前、AWSの提供するGPUインスタンス g2.* に搭載されているGPUはGRID K520というちょっと古いモデルで、PG-Stromは非対応だった。
理由は、一年ほど前にComputing Capability 3.5以降で対応のDynamic Parallelism機能を使うように全面的に作り直したからで、詳細は以下のエントリを参照。

kaigai.hatenablog.com

その後、昨年の10月にAWSは新世代*1GPUインスタンスを新たにリリースした。

japan.zdnet.com

これでPG-Stromの動作要件を満たすようになった上に、特にメモリ搭載量で相応の強化が行われたため、例えばPGconf.ASIAで発表を行った創薬領域の類似度サーチのような、I/Oが支配的でないようなワークロードであれば相応の効果が見込める、ハズである。

発表から少し間が空いてしまったが、p2.*インスタンスで動作検証を行い、当該環境をAMIイメージとして公開してみた。

環境構築

基本的にはAMIを探してポチポチするだけである。


最初に、リージョンを p2.* インスタンス提供済みのリージョンに切り替える。
現状、東海岸の「バージニア北部」、西海岸の「オレゴン」、欧州の「アイルランド」での提供が開始されている模様。東京はまだである*2
とりあえず、今回は「オレゴン」リージョンを使用。

また、AWSアカウントを作成直後は p2.* インスタンスの起動制限数が 0 に設定されているため、p2.xlargeインスタンスを起動できるよう、EC2ダッシュボードの『制限』から『制限緩和のリクエスト』を行う必要がある。豆知識。


次に、コミュニティAMIから『PGStrom』で検索すれば、夜なべして作ったAMIイメージがヒットする。
これは PG-Strom v1.0 + PostgreSQL v9.5.5 + CUDA 7.5環境を CentOS7(x86_64) 上に構築したもので、テスト用のデータも既にセットアップ済みのものである。


次に、起動すべきインスタンスを選択する。
p2.*インスタンスは以下の3種類提供されており、今回は『p2.xlarge』タイプを使用した。
なお、Amazon様も『高性能データベース』用途ときちんと書いていてくれておりますw

一応、PG-StromはマルチGPUにも対応しているが、CPU並列を使えないPostgreSQL v9.5系列ではGPUへのデータ供給の方が先にボトルネックになってしまうので、あまり意味はない。
この辺は、現在開発中の PostgreSQL v9.6 ベースの実装でうまくハンドリングできるようになる。

インスタンスの構成を確認して起動。ストレージはroot区画に24GBをアタッチするだけのシンプルなものがデフォルト構成。

しばらく待っているとインスタンスが立ち上がってくる。

『ec2-user』でログインする事ができる。

動作確認

早速動作確認を行う。ec2-userはlocalhost経由でPostgreSQLサーバに接続できるよう設定されている。

[ec2-user@ip-172-31-24-196 ~]$ psql postgres
psql (9.5.5)
Type "help" for help.

postgres=#

Tesla K80 GPUが搭載されている事が分かる。
物理的には一枚のカード上にGPU 2個を搭載するデバイスなので、一枚のカードを他の誰かとシェアしている事になる。

postgres=# SELECT * FROM pgstrom_device_info();
 id |                       property                        |    value
----+-------------------------------------------------------+--------------
  0 | Device name                                           | Tesla K80
  0 | Total global memory size                              | 11519 MBytes
  0 | max threads per block                                 | 1024
  0 | Maximum block dimension X                             | 1024
  0 | Maximum block dimension Y                             | 1024
  0 | Maximum block dimension Z                             | 64
  0 | Maximum grid dimension X                              | 2147483647
  0 | Maximum grid dimension Y                              | 65535
  0 | Maximum grid dimension Z                              | 65535
  0 | Maximum shared memory available per block             | 48 KBytes
  0 | Memory available on device for __constant__           | 64 KBytes
  0 | Warp size in threads                                  | 32
  0 | Maximum pitch in bytes allowed by memory copies       | 2147483647
  0 | Maximum number of 32bit registers available per block | 65536
  0 | Typical clock frequency in kilohertz                  | 823 MHz
  0 | Alignment requirement for textures                    | 512
  0 | Number of multiprocessors on device                   | 13
  0 | Has kernel execution timeout                          | False
  0 | Integrated with host memory                           | False
  0 | Host memory can be mapped to CUDA address space       | True
  0 | Compute mode                                          | Default
  0 | Alignment requirement for surfaces                    | 512
  0 | Multiple concurrent kernel support                    | True
  0 | Device has ECC support enabled                        | True
  0 | PCI bus ID of the device                              | 0
  0 | PCI device ID of the device                           | 30
  0 | Device is using TCC driver model                      | False
  0 | Peak memory clock frequency                           | 2505 MHz
  0 | Global memory bus width                               | 384 bits
  0 | Size of L2 cache in bytes                             | 1536 KBytes
  0 | Maximum threads per multiprocessor                    | 2048
  0 | Number of asynchronous engines                        | 2
  0 | Device shares unified address space                   | True
  0 | PCI domain ID of the device                           | 0
  0 | Major compute capability version number               | 3
  0 | Minor compute capability version number               | 7
  0 | Device supports stream priorities                     | True
  0 | Device supports caching globals in L1                 | True
  0 | Device supports caching locals in L1                  | True
  0 | Maximum shared memory per multiprocessor              | 112 KBytes
  0 | Maximum number of 32bit registers per multiprocessor  | 131072
  0 | Device can allocate managed memory on this system     | True
  0 | Device is on a multi-GPU board                        | False
  0 | Unique id if device is on a multi-GPU board           | 0
(44 rows)

早速、シンプルな集計のクエリを叩いてみるが、その前にメインのテーブルに対してpg_prewarmを実行しておくが吉。
gp2ストレージタイプの場合、バッファに載っていないデータをロードするのに割と時間がかかるようなので、I/Oネックになってしまうと、CPUとかGPUとかいう次元の遥か手前で引っかかってしまう。

postgres=# SELECT pg_prewarm('t0'::regclass);
 pg_prewarm
------------
     833334
(1 row)

実行計画を確認する。JOIN+GROUP BYをGPUで実行するようにプランが選択された事がわかる。

postgres=# EXPLAIN SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
                                         QUERY PLAN
---------------------------------------------------------------------------------------------
 HashAggregate  (cost=3578311.16..3578311.48 rows=26 width=12)
   Group Key: t0.cat
   ->  Custom Scan (GpuPreAgg)  (cost=14139.24..2873626.84 rows=234 width=48)
         Reduction: Local + Global
         GPU Projection: cat, ax
         ->  Custom Scan (GpuJoin) on t0  (cost=10139.24..2838812.04 rows=98599882 width=12)
               GPU Projection: t0.cat, t1.ax
               Depth 1: GpuHashJoin, HashKeys: (t0.aid)
                        JoinQuals: (t0.aid = t1.aid)
                        Nrows (in/out: 98.60%), KDS-Hash (size: 13.47MB, nbatches: 1)
               ->  Seq Scan on t1  (cost=0.00..1935.00 rows=100000 width=12)
(11 rows)

1億行のテーブルのJOINとGROUP BYで11.3秒要している。
セッションに接続後の初回実行だと、これに加えて、GPUコンテキストの初期化やGPUコードのJITコンパイルの時間も余分にかかるが、オンプレに比べるとAWSの環境ではこの辺が多少もっさりする印象がある。

postgres=# SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
 cat |  count  |       avg
-----+---------+------------------
 nnn | 3845736 | 49.9122204644341
 ccc | 3842975 | 49.9253161146813
 ddd | 3841209 | 49.9346989307005
 aaa | 3848221 | 49.9265219996308
 kkk | 3843481 | 49.9232348058601
 fff | 3845484 | 49.9434949581969
 iii | 3846743 | 49.9227199719054
 jjj | 3846076 | 49.9292863471083
 qqq | 3845646 |  49.945213365697
 hhh | 3842519 | 49.9266693322143
 ttt | 3846725 | 49.9232478784252
 ooo | 3847927 | 49.9346102129999
 zzz | 3847116 | 49.9320751724648
 lll | 3848447 |  49.927865906661
 www | 3846691 | 49.9537192102014
 bbb | 3845315 | 49.9281298849625
 ppp | 3850842 | 49.9149069792416
 eee | 3846285 |  49.931458202446
 xxx | 3845570 | 49.9577920754281
 ggg | 3845044 | 49.9409383291351
 rrr | 3847816 | 49.9341189910578
 uuu | 3845813 | 49.9295202543591
 vvv | 3849157 | 49.9253944163053
 yyy | 3843414 | 49.9364087656463
 mmm | 3848758 | 49.9033622681507
 sss | 3846990 | 49.9213589517191
(26 rows)

Time: 11326.834 ms

一方、PG-Stromを無効にした場合は以下のように、同じクエリの実行に51.0秒を要している。トゥットゥルー。

postgres=# EXPLAIN SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
                                 QUERY PLAN
-----------------------------------------------------------------------------
 HashAggregate  (cost=3937016.94..3937017.26 rows=26 width=12)
   Group Key: t0.cat
   ->  Hash Join  (cost=3185.00..3197517.82 rows=98599882 width=12)
         Hash Cond: (t0.aid = t1.aid)
         ->  Seq Scan on t0  (cost=0.00..1833334.00 rows=100000000 width=8)
         ->  Hash  (cost=1935.00..1935.00 rows=100000 width=12)
               ->  Seq Scan on t1  (cost=0.00..1935.00 rows=100000 width=12)
(7 rows)

postgres=# SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
 cat |  count  |       avg
-----+---------+------------------
 nnn | 3845736 | 49.9122204644368
 ccc | 3842975 |  49.925316114681
 ddd | 3841209 |  49.934698930695
 aaa | 3848221 | 49.9265219996321
 kkk | 3843481 | 49.9232348058579
 fff | 3845484 | 49.9434949581975
 iii | 3846743 |    49.9227199719
 jjj | 3846076 | 49.9292863471066
 qqq | 3845646 | 49.9452133657041
 ttt | 3846725 | 49.9232478784215
 hhh | 3842519 |  49.926669332217
 zzz | 3847116 |  49.932075172461
 ooo | 3847927 | 49.9346102129966
 lll | 3848447 | 49.9278659066662
 www | 3846691 | 49.9537192102029
 bbb | 3845315 | 49.9281298849615
 ppp | 3850842 | 49.9149069792369
 eee | 3846285 |  49.931458202444
 xxx | 3845570 | 49.9577920754269
 ggg | 3845044 | 49.9409383291352
 uuu | 3845813 | 49.9295202543602
 rrr | 3847816 | 49.9341189910555
 vvv | 3849157 | 49.9253944163048
 yyy | 3843414 | 49.9364087656454
 mmm | 3848758 | 49.9033622681445
 sss | 3846990 | 49.9213589517182
(26 rows)

Time: 51031.084 ms

p2.* インスタンスの登場でPG-Stromを使用するための環境の準備はずいぶんお手軽になった。デプロイ即利用できるようAMIイメージも用意してみたので、まだPG-Stromを触った事がないという方は、ぜひ試してみて頂ければと思う。

*1:と言っても、Tesla K80なんですが…。

*2:予定なしとか悲しい事は言わないでほしいけど

2017年の開発ロードマップについて考える

あけましておめでとうございました。(やや出遅れ感)

新年という事で、この一年、どういった技術開発に取り組んでいきたいかをざーっと書き出してみる事にする。
これらのうち、いくつかはPostgreSQL本体機能の強化を伴うものであったりするので、ある程度計画的にモノゴトを進めないといけないワケで…。

PG-Strom v2.0

先ず最優先で取り組むのが、PostgreSQL v9.6への対応。
CPUパラレル実行と、新しいオプティマイザへの対応でかなり大きなアーキテクチャ上の変更を伴ったものの、全体としてはよりシンプルな設計に落とし込む事ができている。

ちなみに、現状だとこの程度までは動くようになっている。
集約演算がGroupAggregateGpuPreAggの二段階に分解されており、GpuPreAggGatherの配下で並列に動作している事に注目。

postgres=# EXPLAIN (ANALYZE, VERBOSE)
           SELECT cat, count(*), avg(aid), max(bid)
             FROM t0
            WHERE aid < 50000 and cid > 50000
            GROUP BY cat;

                                           QUERY PLAN
------------------------------------------------------------------------------------------------
 GroupAggregate  (cost=91050.50..91070.84 rows=26 width=48)
                 (actual time=1754.432..1755.056 rows=26 loops=1)
   Output: cat, pgstrom.sum((pgstrom.nrows())),
                pgstrom.favg((pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                             (pgstrom.psum((aid)::bigint))))),
                max((pgstrom.pmax(bid)))
   Group Key: t0.cat
   ->  Sort  (cost=91050.50..91052.32 rows=728 width=48)
             (actual time=1754.381..1754.524 rows=910 loops=1)
         Output: cat, (pgstrom.nrows()), (pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                                                       (pgstrom.psum((aid)::bigint)))),
                 (pgstrom.pmax(bid))
         Sort Key: t0.cat
         Sort Method: quicksort  Memory: 160kB
         ->  Gather  (cost=90938.54..91015.89 rows=728 width=48)
                     (actual time=1749.313..1753.732 rows=910 loops=1)
               Output: cat, (pgstrom.nrows()), (pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                                                             (pgstrom.psum((aid)::bigint)))),
                       (pgstrom.pmax(bid))
               Workers Planned: 4
               Workers Launched: 4
               ->  Parallel Custom Scan (GpuPreAgg) on public.t0  (cost=89938.54..89943.09 rows=182 width=48)
                                                            (actual time=1670.139..1673.636 rows=182 loops=5)
                     Output: cat, (pgstrom.nrows()),
                             pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                                          (pgstrom.psum((aid)::bigint))),
                             (pgstrom.pmax(bid))
                     Reduction: Local
                     GPU Projection: t0.cat, pgstrom.nrows(), pgstrom.nrows((t0.aid IS NOT NULL)),
                                     pgstrom.psum((t0.aid)::bigint), pgstrom.pmax(t0.bid), t0.aid, t0.cid
                     Outer Scan: public.t0  (cost=4000.00..87793.21 rows=2496387 width=12)
                                            (actualtime=14.663..274.534 rows=500187 loops=5)
                     Outer Scan Filter: ((t0.aid < 50000) AND (t0.cid > 50000))
                     Rows Removed by Outer Scan Filter: 1499813
                     Extra: slot-format
                     Worker 0: actual time=1724.195..1728.058 rows=182 loops=1
                     Worker 1: actual time=1570.952..1573.837 rows=182 loops=1
                     Worker 2: actual time=1738.205..1742.053 rows=182 loops=1
                     Worker 3: actual time=1569.055..1571.961 rows=182 loops=1
 Planning time: 0.907 ms
 Execution time: 1759.557 ms
(25 rows)

また、PG-Strom v2.0では、PostgreSQL v9.6へのキャッチアップだけではなく、いくつか目玉となる機能を準備中である。
一つは、これまで何度か紹介している SSD-to-GPU P2P DMA 機構。そしてもう一つは、BRINインデックスへの対応である。

SSD-to-GPU P2P DMA

SSD-to-GPU P2P DMA (NVMe-Strom) は、NVIDIA社製GPUのGPUDirect RDMA機構を利用したもので、PostgreSQLのデータブロックが格納されているNVMe-SSDのデータブロックからGPUへとダイレクトにデータ転送を行う。ファイルシステムを介する事によるオーバーヘッドや、RAMへの無駄なコピーが発生しないため、スループットを稼げるという特長がある。
現状では、GpuScanワークロード下においてNVMe-SSD 1個から成る区画からのデータ転送に対応しており、シングルプロセス性能で1.4GB/sのスキャン性能を出している。
PostgreSQL v9.6対応の過程で、GpuJoinやGpuPreAggの直下にテーブルスキャンが入る場合、これらのロジックはGpuScanがテーブルをスキャンするための関数を直接呼ぶように改良されているので、特別な事は何もしなくても『ストレージからデータを読んだ時点で既にJOIN/GROUP BYが完了している』という状態を作り出す事はできるはず。

PG-Strom v2.0に向けた課題はSoft-RAID0/1への対応。Linuxの場合、基本的には128KB単位で順番にストライピングがかかっているだけなので、技術的にはそう難しい話ではないと考えている。
DC用途向けに、PCIe x8スロット接続で5~6GB/s程度のSeqRead性能を持つNVMe-SSD製品が各社から出てきているので、計算上は、SSD二枚から全力でGPUにデータを流し込む事ができれば、GPUの持つPCIe x16スロットの帯域を飽和させられる事になる。

BRINインデックス対応

BRINインデックス自体はPostgreSQL v9.5から搭載されている機能で、特に時系列データのように

  1. ある一定範囲の値を持つデータが物理的な近傍に集まっている
  2. データの更新頻度が小さい
  3. データサイズが大きい

といった特徴を備えたデータセットに向いており、例えば、センサデータをPostgreSQLに収集して解析するといったワークロードに有効な機能。

BRINインデックス自体は、永安さんのこちらの記事が詳しいです。
pgsqldeepdive.blogspot.jp

PG-Stromとしては、搭載RAMが比較的小さなGPUを使うという事もあり、B-treeのようなランダムメモリアクセスを前提としたインデックスへの対応は厳しい。
ただ、条件句の評価はGPUの数千コアを使って並列処理が可能であるものの、インデックスの選択率が高くなると分の悪い勝負なので『このブロックは明らかに該当行なし』という事が分かっているなら、それを読み飛ばしたい。

GpuScanがBRINインデックスを理解し、必要のないブロックを読み飛ばす事ができるようになれば、例えばIoTのキーワードに絡めてセンサデータの集積・解析用途に、という使い方もできるハズである。
特に、PG-Stromはカラムナーを前提としたDWHではないので、生データをそのまま処理させても高速化できるという点は強みになるだろう。

PL/CUDA

PL/CUDAに関しては、言語バインディングに関してする事は(できる事は)多くないので、その周辺領域を拡充していきたい。

一つは、PostgreSQLにおける可変長データの1GB越え。
現状、全ての可変長データの基盤となっている varlena 構造は、最大でも1GBまでのデータしか持てないため、PL/CUDA関数の引数としてarray-matrixを渡す時には、例えば問題領域をうまく分割してデータサイズを1GB未満に抑えてやらないといけない。

昨年10月のCBI学会で発表した研究でも、1000万件の化合物データ(1.3GB)をロードするにはサイズが大きすぎたので、安全マージンも見て4分割した上でGPU側へロードしている。

しかし、昨今のGPUでは10GB近く、あるいはそれ以上のメモリを搭載するのが常識的になりつつあり、問題領域を1GB以内に抑えねばならない、、、というのはユーザにいかにも不都合である。

先にpgsql-hackersにデザインプロポーザルを投げたところ、Robert Haasから『varlenaとは別の体系で可変長データを保持するフォーマットを作成すべき』とサジェスチョンがあり、自分もこの方針には同意。3月のcommit-festまでにはパッチを投稿し、2018年リリース予定のPostgreSQL v11でのマージを目指したい。

もう一つは、現状、複雑な計算ロジックを個々のユーザ毎に書かねばならないという点である。
2017年1月時点でPL/CUDAを実証できたワークロードとしては、k-NN法類似度サーチや、k-meansクラスタリングがあるが、例えば MADLib のような統計解析パッケージで提供されているアルゴリズムの、全部とは言わないまでも、使用頻度が高く計算負荷の高いものをGPUで計算するようパッケージ化できれば、ユーザの裾野はより広がるだろうし、仮にカスタマイズが必要となっても骨格となるアルゴリズムGPU実装が既に存在する事で省力化が可能となるハズである。

PG-Strom v3.0へ向けた種まき

In-memory Columnar Cache

NVMe-SSDとの密連携の他にI/O系処理を高速化する方策として、通常のPostgreSQLテーブル(行形式)の脇に、予めDMAバッファ上に列形式にデータを再編したキャッシュを持たせる機構を考えている。
このキャッシュはBackground-workerにより非同期で作成され、そのため、スキャンする区間のうち一部領域しか列形式のキャッシュが構築されていないかもしれない。しかし、その場合でもPG-StromはGPUで行形式データを扱えるので(多少のパフォーマンス差に目をつぶれば)大きな問題とはならない。

データ本体とキャッシュを別に持つ場合、必ず一貫性制御が複雑で頭の痛い問題として立ちはだかる。
イデアとしては、これを避けるために ALL_VISIBLE フラグが 1 であるブロックのみを列形式キャッシュに持たせる。
ALL_VISIBLE=1であるブロックは、MVCC制御に関わらず、全てのレコードが全てのトランザクションから可視である事が保証されている。そのため、複雑な同時実行制御に頭を悩ませる必要はなく、全てのレコードの内容を単純に列データとして展開すればよい。

問題は、PostgreSQL側でテーブルが更新され、ALL_VISIBLE フラグが 0 にクリアされた時のinvalidation処理である。
現状、ここにフックを挟む事はできないので、PostgreSQL側の機能強化を行う必要がある。
デザインプロポーザルを出し、この処理を行うにふさわしい場所とフックの仕様を固めていきたい。

GpuSort(+LIMIT)

PG-Strom v2.0では、実はGpuSort機能の廃止を予定している。これは、ソートという問題の性質上、ある程度問題規模が大きくならないとGPUによる処理時間メリットが出てこない一方で、一度にGPUでソートを行う件数が多くなればなるほど、初期データのローディングに時間がかかるようになり、非同期・多重処理のメリットを得にくいからである。
そのため、GPUでソートを行うセグメントサイズを小さくして単位ローディング時間を短くする一方で、CPUでのMergeSort処理の割合が大きくなるか、それとも、セグメントサイズを大きくするかというジレンマに悩まされてきた。
(で、最終的にはそれほど速くならない事が多かったり・・・。)

根本的な原因は、GPUで処理を行ってもソートではデータ件数を減らす事が不可能な点にある。なので、GpuSortを有効に活かすには、何がしか『データ件数を減らせる』パターンでのソートに限った方が利口だ。
例えば、LIMIT句で『上位xx件を取り出す』という事が明らかな場合に限り、GpuSortを使用するというパターンであれば十分に効果を発揮する事ができるだろう。
これにはPostgreSQL本体側で、『LIMIT句でxx件のデータが要求されますよ』という事を下位のノードに伝えてやる仕組みが必要だが、PostgreSQL v10 向けにパッチを出しており、ある程度レビューも進んでいるので間に合うだろうとは踏んでいる。

こんな感じで2017年の開発ロードマップについて考えてみたが、さてさて、大晦日に振り返ってどの程度きちんとやれているでしょうか。

Beyond the 1GB limitation of varlena

This article is a part of the PostgreSQL Advent Calendar 2016.

According to the request by Joe Conway (@josepheconway), I wrote this article in English.

I like to share the discussion we had at the PostgreSQL developer unconference on the day before PGconf.ASIA 2016, and the related idea of mine.

PostgreSQL supports variable length data types; like text, bytea, json, numeric, array and so on.
These data types individually have their own characteristics and own internal format, however, all of them are built on a common structure to represent variable length field; that is called varlena.

It has very simple internal format. Contents of the variable length fields are followed by either 1-byte or 4-bytes header.

We can identify the header type by the least bit of the first 1 byte. If least bit of the first byte, it means 1-byte header, elsewhere, 4-byte header.
In case of 4-byte header, the second bit is also used to show whether it is compressed or not. Therefore, rest of 30-bits can be available to represent the contents length, so, maximum length of the variable length field is 1GB.
A special case exists if the first byte is 00000001. It is an external TOAST reference which consists of OID of TOAST table and unique ID within TOAST table.

Below is the source code comment at include/postgres.h.

    :
 * Bit layouts for varlena headers on little-endian machines:
 *
 * xxxxxx00 4-byte length word, aligned, uncompressed data (up to 1G)
 * xxxxxx10 4-byte length word, aligned, *compressed* data (up to 1G)
 * 00000001 1-byte length word, unaligned, TOAST pointer
 * xxxxxxx1 1-byte length word, unaligned, uncompressed data (up to 126b)
   :

My concern is we have no way to represent a larger variable length datum more than 1GB.

It is quite natural for users who want to process the maximum available data size as long as system capability allows. It is not uncommon for the recent GPU models to have 10GB class device memory or more.
On the other hands, due to the restriction of variable length datum of PostgreSQL, we cannot have a datum larger than 1GB. It also means we cannot provide a large matrix (= 2D-array) onto PL/CUDA function at once. It is an unignorable restriction if a problem user want to solve by PL/CUDA is unavailable or expensive to split into multiple portions.

According to the background, we discussed a few options to support variable length datum larger than 1GB.

64bits varlena header

The first idea was straightforward but got less interest because 99% of existing variable length data types are satisfied with 1GB limitation.
If we have one more varlena header, VARDATA() and VARSIZE() which are widely used in the PostgreSQL core and extensions need to have branch operation inside the macro, and it is not easy to justify the penalty for this niche usage.

Use of large-object and its ID instead

The second idea was a suggestion from audience. Now PostgreSQL has a feature of large object which allows to store up to 4TB data chunk with a unique identifier. If PL/CUDA function supports ID of large object, instead of varlena datum, PL handler can extract larger data chunk by itself.
However, here is a problematic scenario. It requires users to construct large objects preliminary. It is inconvenient when user wants to deal with a large matrix constructed on the fly. For example, the query below constructs a matrix based on the underlying table scan with qualifier. The qualifier can be changed for each execution.

SELECT * FROM matrix_unnest(
  (SELECT my_plcuda_function(Q.matrix, D.matrix)
     FROM (SELECT array_matrix(a,b,c) matrix
             FROM table_q
            WHERE tag LIKE '%abc%') Q,
          (SELECT array_matrix(x,y,z) matrix
             FROM table_d
            WHERE tag LIKE '%xyz%') D
  )
)

In addition, creation of a large object makes unnecessary i/o; which leads performance slow-down.

Special type that includes indirect reference

A overall consensus was to define a special data type which support indirect references to data chunks less than 1GB. If only narrow use-case wants to have a datum larger than 1GB, it is quite natural to define a special data type for the purpose.

My interest is representation of matrix in database.

In case of matrix, it has less necessity to have all the items in a flat memory chunk.

If we have 8GB of matrix, we can split it into 9 portions to keep individual chunk size less than 1GB.

Then, once we define a special matrix type that consists of small metadata and indirect references to the chunks less than 1GB, it is available to represent a big matrix larger than 1GB as a usual SQL data type.

A remaining problem is serialization/deserialization when the matrix data type is saved/loaded.
Right now, PostgreSQL saves contents portion of the variable length datum onto the storage as is, thus, pointers to reference sub-matrix has to be serialized appropriately, however, we have no infrastructure to manipulate type specific data structure on toast_insert_or_update() which set up physical tuple for INSERT/UPDATE.
Likely, pg_type needs to have an optional callback functions for serialization/deserialization. It shall be a mandatory requirement if data structure has indirect reference to other memory chunks.

I expect the serialization callback will use TOAST relation to save a large but less than 1GB chunks, then put its unique ID instead of the pointers. We will be able to have another advantage in this approach, because all the sub-matrix we have to update are the portion actually updated. If some of sub-matrix were not updated, we don't need to delete old toast pages and insert new ones. It will make a performance benefit than existing flat varlena structure.

The timing for deserialization needs a bit more consideration. Because heap_tuple_fetch_attr handles deserialization of the existing flat varlena, but no type OID is supplied to the function. It is not a good option to change the function prototype because many existing code already uses this function without type OID.
We have two another options here. The first one packs type OID within the serialized structure. It needs to define a new VARTAG_* label to distinct from the existing flat varlena. The second one is delayed load because indirectly referenced data chunks will not be used without functions/operators which support the data type. It enables not to load unreferenced chunk, however, it is uncertain whether functions/operators can manipulate a value supplied as an argument. *1

Last, even if we can have variable length datum larger then 1GB from the viewpoint of data format, it is never small data chunks. It involves not a small amount of i/o (or memory copy) stuff.
Therefore, it is a significant to have special optimization based on the knowledge for usual use case of the types.

For example, some workloads takes sparse matrix which have small amount of non-zero values, but mostly zero. In this case, type may be able to assume empty sub-matrix are all-zero instead of data size reduction.

Diagonal matrix is also often used, in case when valid values are located around the diagonal axis only.

I hope making a proof of the concept patch near future, then have a discussion at pgsql-hackers.

*1:If not acceptable, we may need to load sub-matrix multiple times when a particular matrix object is referenced by multiple functions/operators.

PGconf.SV 2016 and PL/CUDA

I've participated PGconf Silicon Valley 2016 held at the South San Francisco Conference Center.

I could have a presentation here. Its title is PL/CUDA - Fusion of HPC Grade Power with In-Database Analytics.
Its slides are below:
https://www.slideshare.net/kaigai/pgconfsv2016-plcuda/

What is PL/CUDA?

PL is an abbreviation of Procedural Language. In PostgreSQL, we have some PL options like PL/pyhton, PL/R and so on.
This functionality allows to implement SQL functions with different program languages individually.
For example, this SQL function definition is described with "SQL" language.

CREATE OR REPLACE FUNCTION sample(int, int)
RETURNS int
$$
  SELECT $1 + $2;
$$ LANGUAGE 'sql'

We can describe the code block between the $$ marker in SQL-language.

It means we can use other programming language for SQL function definition if another LANGUAGE is given.
So, it is never strange to have CUDA language here.

PL/CUDA is a procedural language which supports CUDA code block for user defined SQL function.
It allows to embed a heavy computing portion within a SQL query, and invokes GPU kernels transparently.

As literal, DBMS is database management system. It implies heavy computing is not a primary workload for DBMS, and it is almost correct. So, people who want to run heavy analytic algorithm usually exports the data from the database for more computing performance. However, here is two downsides. The first is network traffic to export the data set. If DBMS would have enough computing capability, all we have to export is results of the analytic algorithm. The other is loss of SQL flexibility on the post-process of the analytic algorithm. If we could finish the algorithm core in database, it is quite easy to join the result with other tables, to make a summary using GROUP BY or window functions.

PL/CUDA is a solution to solve these problems. The code block in the function declaration shall be put into GPU kernel, then PG-Strom's infrastructure build the GPU kernel on the fly and execute on GPU devices with thousands threads. Its function arguments and results are automatically move to/from GPU devices.

Actually, PL/CUDA changes the concept of PG-Strom a bit. It originally intended to accelerate SQL by GPU transparently, thus existing SQL shall be accelerated with no query changes.
On the other hands, PL/CUDA assumes the core logic of algorithm is implemented manually for each use case.
Why we changed this concept? It is based on user's feedback.
Even if we don't need to describe CUDA code, right now, nobody implements complicated analytic algorithms in SQL, thus, it means user has to write up a new code to process them either SQL or CUDA. I'm not certain whether SQL is suitable language for this purpose, and most of known algorithm assumes procedural languages.
So, we thought writing a CUDA code may be a valuable option in some use cases, especially, processing of advanced algorithms. We call this type of GPU utilization manual optimization mode, in contradistinction to the existing automatic code generation from SQL; automatic optimization mode.

Array-Matrix

One other thing we need to pay attention is data-format and memory bus utilization rate when GPU kernel handles the data-format.
PostgreSQL adopts row-oriented data format. It is an I/O optimal data format for transactional workloads.
On the other hands, it has a downside when we reference individual columns within a record. Because the location of column is not predictable unless reading the record actually, we have to walk on records from the head, thus, it takes larger number of memory access.
It is unignorable penalty when we run a logic on GPU which has relatively smaller memory cache. In case of column-oriented data format, location of the data can be determined uniquely by pointer of the array and thread-id.

In addition, column-oriented data format has an advantage for GPU devices from the standpoint of memory bus utilization rate.
Unit size of GPU's memory transaction is usually 32bytes/256bits. If multiple processor cores requires coalesced memory location simultaneously, these multiple memory accesses are consolidated to a single memory transaction, then, a single memory transaction will load 8 x 32bits variables at once. This scenario pulls up the maximum utilization rate from the memory bus of GPUs; that is much higher than CPU/RAM bandwidth.
When we run a particular calculation by multiple processor cores, it tends to reference same columns simultaneously. Columns are scattered in row-oriented data structure, so it is hard to pull out the maximum utilization rate of memory bus.

However, columnar-format is not always winner. Transformation from row-format to column-format needs extra cost prior to GPU invocation, and it is hard to justify the data transformation cost for usual SQL operations like JOIN, GROUP BY, etc...
One candidate is highly computing intensive workload like advanced analytic algorithm which we will implement using PL/CUDA.

Fortunately, PostgreSQL already has a data-type which is almost equivalent to the column oriented data format. It is array data type, and we call the array data type with no NULLs and consists of fixed-length data type array-matrix.
An aggregate function array_matrix(variadic datatype[]) consumes multiple rows then produce a 2-dimensional array data. Any elements are re-ordered per column, and data location is predictable by the offset and width of the element data type.

We usually call PL/CUDA function with array-matrix arguments.

Usecase.1 - Similarity Search on Drug Discovery

The first usecase I introduced is a similarity search workloads on drug-discovery *1.

First of all, I like to introduce the background of the problem we tackled.
When researcher tried to find out a candidate of new drug, it has to be "active" to the target protein which is related to the target disease. Some chemical compounds are inactive, some other chemical compounds are active but toxicity. Academic papers, daily published by somewhere, report relationship between a target protein and chemical compounds. So, researcher can gather the public data of chemical compounds relevant to a particular protein.

On the other hands, people said the total number of chemical compounds we can obtain as a test reagent and investigate are about 10,000,000. It is too expensive to purchase and test them in laboratory, with brute-force approach. So, researcher is motivated to pick up chemical compounds which have higher probability of active to the target protein.
It is an approach to search "similar" chemical compounds to the known ones; which are already reported in academic papers and etc...

Similarity is to define a distance between two chemical compounds. Euclid distance is one approach usually used to define geometrical distance between two points.
We used a combination of Jaccard index *2 and k-NN (nearest neighbor) method in our trial.

The k-NN method is used to define distance between a set of items and an item. Once we calculate all the distance between items, than picks up the largest k distances and considers its average the representative distance between them.
Jaccard index can be used when characteristics of two items are described in a bitmap. We consider the ratio of bitmap-AND and bitmap-OR as similarity of two items.

According to the logic, if number of query compounds set is Q, we once need to calculate Q similarities for each database item; that is about 10M. Than, these similarity must be sorted for each database item. Thus, total amount of computing order is almost O(DxQ) + O(DxQlogQ).
It is not small. If Q is 1000, we need to handle 10billion combination to process the logic.

The PL/CUDA function we implemented (knn_gpu_similarity()) did nothing special. It once computes the similarity between two chemical compounds by a processor core for each combination, then run the bitonic-sorting to pick up the largest k-similarities.
This PL/CUDA function takes two array-matrix that are constructed from the relation scan on the fly, or statically pre-built.
Once similarity search gets completed, it unnest the matrix to a generic stream of records by matrix_unnest, then processes the results set using usual SQL features (tables JOIN and window function here).

Its performance is quite impressive. We compared to the same logic based on CPU.
In case of the largest query set (1000), it reduces the response time more than x150 times faster.
It may have a power to change the workstyle of researcher. 3,000sec (=50min) is a timescale for batch job. It is too long to wait in from of the workstation. On the other hands, 20sec is scale for interactive operations; which allows try&error on research process.


Usecase.2 - Data Clustring (k-means)

The second usecase is clustering analysis.
k-means is a representative method for non-hierarchical cluster analysis, used for unsupervized learning.

This algorithm tries to get the final clustering with iteration of simple operations.
First, it assigns a cluster for each item randomly as initial state.

Second, makes a centroid for each cluster according to the initial state.

Next, choose the cluster of every items according to the centroid updated on the last step.

Next, update the centroid according to the latest cluster; updated on the last step.

Then, break the iteration if cluster assignment gets converged or reached to the threshold of the iteration.

The data-set we tried to apply k-means() clustering is a collection of city traffic data at Aarhus, Denmark, by the CityPlus project.
iot.ee.surrey.ac.uk

This dataset contains traffic data (car count, avg speed, ...) in 449 observation point all around the city. The total amount of sensor data is 13.5M records.
In this trial, we applied k-means() algorithm on the raw dataset, then make a summary for each observation point. We assumed each observation points belong to the most frequent cluster.
Then, mapped it on the google map using static map APIs.

The result is below:

It looks to me the blue observation points are bypass highway; many car count and high speed.
On the other hands, green line seems to me the way to the downtown, thus many car but slow speed.
I'm not sure the characteristics of the red line, but appears on the beltway frequently.

Because all the clustering calculation is executed inside of the database, it is quite easy to control the input data.
Below is the same diagram when weekdays and weekend.

  • weekdays


  • weekend

Here is a little difference between them, however, the diagram of weekend categorize the road to the downtown with same color of the highway. So, it has less traffic on the weekend, thus, car speed is higher than the case of weekday.

For the arrangement, all I did is just adding a WHERE clause on the query to kick k-means clustering. Utilization of the SQL flexibility is one of the biggest advantage of PL/CUDA, in addition to the performance benefit.

SELECT report_id, k, c
  FROM (SELECT report_id, k, c,
               row_number() OVER (PARTITION BY report_id
                                  ORDER BY c DESC) rank
          FROM (SELECT report_id, k, count(*) c
                  FROM matrix_unnest(
                        (SELECT gpu_kmeans ( array_matrix(
                                               int4_as_float4(report_id),
                                               avg_measured_time,
                                               avg_speed,
                                               vehicle_count),
                                             5)
                           FROM tr_rawdata
                          WHERE extract('dow' from timestamp) in (0,6)
                        )
                       ) R(report_id int, k int)
                 GROUP BY report_id, k
               ) __summary_1
       ) __summary_2
   WHERE rank = 1;

*1:KaiGai Kohei, Yoshimori Atsushi, Efficient Similarity Search Using Multiple Reference Molecules on PG-Strom architecture, CBI Annual Meeting, 2016

*2:also called Tanimoto index

PL/CUDAでk-means法を実装する

前回のエントリでは、CBI学会で発表を行った、PL/CUDAによる類似化合物の検索について説明した。

今回は、コレとはまた別のワークロードに対する応用という事で、クラスタリング処理のIn-Database実装に挑戦してみた。
トライしてみたのは k-means法 によるクラスタリング。非階層クラスタリングの領域では最も頻繁に使用される(と、言われている)アルゴリズムで、計算量もそこそこ大きい。

k-meansクラスタリングとは

教師なし学習方式の一つで、所与のデータ群を一定数(=k個)のグループに分類するためのアルゴリズムである。
以下のステップを一定回数、またはクラスタに属するデータ群に変化がなくなるまで繰り返す。

1.初期クラスタをランダムに設定する。

2.各クラスタの中心点を計算する。

3.各データ要素の属するクラスタを更新する。各データ要素はクラスタ中心点が最も近傍であるクラスタに割り当てられる。

4.新しいクラスタ定義に基づいて、各クラスタ中心点を再計算する。

5.以下繰り返し・・・。

これをPL/CUDAで実装するとどうなるか?

基本的には上記のアルゴリズムを愚直にCUDAで実装する事になるので、特別な事はしていない。
一回のPL/CUDA関数呼び出しの中で何度も繰り返し処理を行う形になるので、メインのGPU Kernelは1スレッドで起動し、Dynamic Parallelismを使って別のGPU Kernelを起動する形にする方が合理的である。

関数定義の中身に興味がある方は、こちらをどうぞ。
toybox/gpu_kmeans.sql at master · kaigai/toybox · GitHub

In-databaseでクラスタリングする

以下の図はお試し的に10000件のランダムデータを5つのクラスタに分割してみたもの。綺麗に分類されている。

ただ、ランダムデータだけだと面白くないので、来週のPGconf.SVでの発表に向け、そこそこ件数が多くてクラスタリングに適したデータは無いかと探してみたら・・・あった。

iot.ee.surrey.ac.uk

FP7ファンドの支援を受けたCityPulseというプロジェクトが収集しているデータで、デンマークのオーフス(Aarhus)市の自動車通行状況を調査したデータが公開されている。
これは、ある二地点間を通過した自動車の台数や平均速度がタイムスタンプと共に収集されているもので、2014年2月~6月の四ヵ月間のデータが全部で1350万レコード。
これをクラスタリングにかけると、どこの道路とどこの道路の通行状況が似ているのか分かるはず。

と、いう事でやってみた。

まず、PL/CUDAを用いて実装したgpu_kmeans()関数は以下のように定義されている。
第一引数にクラスタ化すべきデータとそのID値をreal[]行列として与え、第二引数にはクラスタ数を与える。

CREATE OR REPLACE FUNCTION
gpu_kmeans(real[],    -- ID + Data Matrix
           int,       -- k-value (number of clusters)
           int = 10,  -- max number of iteration
           int = 1)   -- seed of initial randomness
RETURNS int[]

これをクエリの中で使用すると以下のようになる。
少しGoogle Map Static APIを勉強して、SQLから直接URLを生成するようにしてみた。
ただ、URL長の制限が8KBという事で、線を何百本も引くような地図は作れないため、泣く泣く125件分だけに限定して地図を生成している。

WITH
summary AS (
SELECT report_id, k, c
  FROM (SELECT report_id, k, c,
               row_number() OVER (PARTITION BY report_id
                                  ORDER BY c DESC) rank
          FROM (SELECT report_id, k, count(*) c
                  FROM matrix_unnest(
                          (SELECT gpu_kmeans(array_matrix(
                                               int4_as_float4(report_id),
                                               avg_measured_time,
                                               avg_speed,
                                               vehicle_count),
                                             5)
                             FROM tr_rawdata
                          )
                       ) R(report_id int, k int)
                 GROUP BY report_id, k
               ) __summary_1
       ) __summary_2
   WHERE rank = 1
),
location AS (
SELECT point_1_lat, point_1_lng,
       point_2_lat, point_2_lng,
       CASE k WHEN 1 THEN 'red'
              WHEN 2 THEN 'blue'
              WHEN 3 THEN 'green'
              WHEN 4 THEN 'purple'
              ELSE 'orange'
       END col
  FROM summary s, tr_metadata m
 WHERE s.report_id = m.report_id
),
path_definition AS (
SELECT 'path=color:' || col || '|weight:3|' ||
       point_1_lat::text || ',' || point_1_lng::text || '|' ||
       point_2_lat::text || ',' || point_2_lng::text path_entry
  FROM location
 LIMIT 125 -- becuase of Goole Map API restriction
)
SELECT 'http://maps.google.com/maps/api/staticmap?' ||
       'zoom=11&' ||
       'size=640x480&' ||
       'scale=2&' ||
       string_agg(path_entry, '&') ||
       '&sensor=false'
  FROM path_definition;
$ wget -O map.png "`psql traffic -At -f ~/traffic.sql`"

その結果がこれ。

どうやら、青の区間は高速道路かバイパス道路で、比較的自動車の台数が多くスピードが出ているクラスタに見える。
一方、緑の区間ダウンタウンに向かう道路なので、慢性的な渋滞に悩まされているのだろうか?
赤は少し判断に迷うが、そのどちらでもない、といった所のようである。


また、これは中核となるkmeansクラスタリングSQLの中に埋め込んでいるので、抽出条件を変えるだけで簡単に母集団を切替える事ができる。
以下の例は、個々のデータのタイムスタンプを用いて、ウィークディ(月曜~金曜)とウィークエンド(土曜、日曜)でそれぞれクラスタリング、図を作ってみたケース。

■平日

■週末

微妙な差ではあるが、週末のケースでは、ダウンタウンに向かう道路がハイウェイと同じ色に分類されている。
平日に比べると交通量が減る分、自動車がスイスイ進むという事だろうか。

閑話休題。今週のPGconf.SVでは、この辺のネタも交えてPL/CUDAメインで話をしてくる事になります。