数値積分のための点過程とサブサンプリング

早川知志（University of Oxford）
2021/8/19
＠実社会に見られる複雑なネットワークと無限粒子系の交差点 II

1 / 28概要♣ 空間 X 上の Borel 確率測度 µ に対して、数値積分n∑∫wi f(xi ) ≈i=1

Xf(x) dµ(x)（wi ∈ R：重み、xi ∈ X ：サンプル点）を考える
この f がある RKHS（再生核 Hilbert 空間）に入っているとき• 誤差を小さくするための「良い点配置」は？
• そのときの収束レートは？という話題について、サンプリングの観点から議論する

2 / 28正定値カーネルと RKHS

♣ 集合 X 上の正定値カーネル k : X × X → R(⇔ 任意の (xi )ni=1 ⊂ X について (k(xi , xj ))ni,j=1 が半正定値対称)ex) Euclid 内積 x⊤ y、Gaussian exp(−|x − y|2 /σ 2 ) など
♣ k を再生核にもつRKHS H ⊂ RX ：• span{k(·, x)}x∈X ⊂ H は H 内稠密• ∀x ∈ X , f ∈ H について f(x) = hk(·, x)i（再生性）
をみたす実 Hilbert 空間

3 / 28

カーネル求積♣ カーネル求積（kernel quadrature/cubature）
X を位相空間、µ を Borel 確率測度として µ の求積公式
Qn (f) =n∑( ∫)f(x) dµ(x) =: µ(f)wi f(xi ) ≈X

i=1の最悪誤差（worst-case error）を次で定める：
wce(Qn ) := sup |Qn (f) − µ(f)|∥f∥≤1（ただし k·k は RKHS のノルム）▷ この最悪誤差を小さくする Qn が欲しい4 / 28

Why kernel quadrature?

♣ カーネル求積を考えるメリットは？
• 古典的な例（Sobolev 空間など）を含んでいる• wce が計算できる（理論保証！）⟨⟩2wce(Qn )2 = sup w⊤ k(·, x) − µ(k(·, ·)), f∥f∥≤1= w⊤ k(x, x)w − 2w⊤ µ(k(·, x)) + µ(µ(k(·, ·)))
• カーネルさえ導入してしまえば、空間を一般化できる• cubature on Wiener space [Lyons and Victoir, 2004]• 時系列データのカーネル [Kiraly and Oberhauser, 2019]5 / 28

Mercer 展開♣ 積分作用素∫K : L2 (µ) → L2 (µ);

f 7→Xk(·, x)f(x) dµ(x)のスペクトル分解を考えると、適当な条件下でk(x, y) =∞∑σm em (x)em (y)m=1と表示できる、ただし• (σm , em ) は K の eigenpair で σ1 ≥ σ2 ≥ · · · ≥ 0√
2∞• (em )∞
m=1 ⊂ L (µ), ( σm em )m=1 ⊂ H はそれぞれ ONB▷ 経験則：最適な Qn は大体 wce(Qn )2 = O(σn )6 / 28

カーネル求積：手法の分類♣ (xi )ni=1 の配置が肝要：大きく分けて二流派がある• 最適化• kernel herding [Chen et al., 2010]：wce への貪欲法• SBQ [Huszár and Duvenaud, 2012]：↑ の発展形

✓安定性・✓一般化・4理論保証・4固有値減衰• サンプリング• optimized sampling [Bach, 2017]：重みつき分布から iid• DPP [Belhadji et al., 2019]：projection DPP• CVS [Belhadji et al., 2020]：点数で条件付けた DPP

4安定性・4一般化・✓理論保証・✓固有値減衰• +α• Mercer + LP [Hayakawa et al., 2021]：後に解説（✓4✓✓）▷ 本講演はサンプリング +α について
7 / 28

Optimized sampling♣ Optimized sampling [Bach, 2017]
初めて kernel quadrature の精度を固有値減衰から示した
Theorem∑σm2λ > 0 に対し、分布 dqλ ∝ ∞m=1 λ+σm em dµ を考える∑∞
σmdλ := m=1 λ+σm , n ≥ 5dλ log(16dλ /δ) の時、(xi )ni=1 ∼iid qλ は()
P wce(Qxn ,opt ) ≤ 4λ ≥ 1 − δをみたす、ただし Qxn ,opt は x = (xi )ni=1 に対して最適な公式()α −α )例えば σm = O m−2α のとき、wce(Qoptn ) = O((log m) m4Mercer 展開が必要・4サンプリングが難しい？
8 / 28

Kernel quadrature with DPPs

♣ Kernel quadrature with DPPs [Belhadji et al., 2019]∑
射影カーネル κn (x, y) := nm=1 em (x)em (y) を用いる
DPP(κn )：p(x) ∝ det(κn (x, x)) で得られる行列式点過程
Theorem∑nrn := ∞m=n+1 σm = o(1/n) のとき (xi )i=1 ∼ DPP(κn ) は[]
E wce(Qxn ,opt )2 ≤ 4nrn + o(nrn )をみたす
4Mercer 展開が必要・4理論保証が弱い（経験的には優秀）

9 / 28

Kernel quadrature with DPPs

♣ 数値実験 [Belhadji et al., 2019, Figure 1](a) Sobolev space, d = 1, s = 1(b) Sobolev space, d = 1, s = 3

10 / 28

Continuous volume sampling♣ Continuous volume sampling [Belhadji et al., 2020]
CVSn (k)：p(x) ∝ det(k(x, x))（ただし点数 n）で得られる点過程
Theorem∑∞
Bn := minm∈{1,...,n} ((n − m + 1)σn+1 )−1 ℓ≥mσℓ についてn(xi )i=1 ∼ CVSn (k) は][E wce(Qxn ,opt )2 ≤ (1 + Bn )σn+1をみたすσm ∼ m−2α (α > 1/2) や σm ∼ β m (β ∈ (0, 1)) なら Bn = O(1)✓Mercer 展開が不要・4サンプリングが非常に難しい
11 / 28

サンプリング：まとめ♣ まとめ手法
OptimizedDPPCVS

wce(Qxn ,opt )2
σ log n⌉∑⌈n/n ∞
m=n+1 σmσn+1Mercer
要要不要サンプリング？
3O(n ) + 棄却 1
MCMC 近似？2♣ 問題意識• サンプリングの実用性• 重み (wi )ni=1 が非負とは限らない
12

Hough et al. [2006,Algorithm
18] & Bardenet and Hardy [2020, Section 2.4]()
Mixing time O n5 log n [Rezaei and Gharan, 2019]12 / 28

+α：正重みカーネル求積♣ 正重みカーネル求積∑n重み (wi )ni=1 を wi ≥ 0,i=1 wi = 1 をみたすようにすると…
• 数値的に安定、k にオーバーフィットしない• 正関数の積分が正になり、定数関数は exact• Qn が確率測度になる• 解釈しやすくなる• カーネル求積の iteration が爆発しない（cf. cubature on Wiener space）

新しい結果の紹介 [Hayakawa et al., 2021]
Hayakawa, S., Oberhauser, H., & Lyons, T. (2021).
Positively Weighted Kernel Quadrature via Subsampling.
(arXiv:2107.09597).
13 / 28

Cubature♣ （古典的な）cubature∑Ω ⊂ Rk と Ω 上の確率測度 ν について、Qn = ni=1 wi δxi が次数t の cubature 公式∫⇔ Qn (f) =f(x) dν(x), ∀f：t 次以下の多項式Ω♣ 一般化 cubature∑
空間 X とその上の確率測度 µ について、Qn = ni=1 wi δxi がテスト関数 φ = (φ1 , . . . , φd )⊤ に対する cubature 公式∫⇔ Qn (φ) =φ(x) dµ(x)X▷ Tchakaloff の定理で正重み cubature の存在が保証
14 / 28

Cubature の構成♣ どうやって cubature を構成するか• 古典的な場合：各論的な代数的構成• 離散測度の場合：recombination [Litterer and Lyons, 2012]• 一般の場合：ランダムサンプリング +LP [Hayakawa, 2021]
Theorem (Hayakawa 2021)
φ = (φ1 , . . . , φd )⊤ ∈ L1 (µ), X1 , X2 , . . . ∼ µ が独立なとき(∫)
Pφ(x) dµ(x) ∈ conv{φ(X1 ), . . . , φ(XN )} → 1 (N → ∞)X15 / 28

Meta-Algorithm♣ 提案する正重みカーネル求積の構成手順(a) n − 1 個のテスト関数 φ = (φ1 , . . . , φn−1 ) ∈ L1 (µ) を RKHSをよく近似するようにとる∫(b) 積分値 X φ(x) dµ(x) を計算する(c) X1 , . . . , XN ∼ µ を独立にサンプルする∑(d) LP を解き、Qn = ni=1 wi δxi であって∫Qn (φ) =φ(x) dµ(x)と wi ≥ 0,∑nXi=1 wi = 1 をみたすものをとる▷ 各ステップ近似で代用可能
16 / 28

Cubature for eigenfunctions♣ 再掲：Mercer 展開k(x, y) =∞∑σm em (x)em (y)m=12（σ1 ≥ σ2 ≥ · · · ≥ 0, (em )∞m=1 は L (µ) の ONB）

Theorem
e1 , . . . , en−1 を exact に積分する n 点正重みカーネル求積公式
Qn について、C := supm≥1 kem k∞ < ∞ のとき
wce(Qn ) ≤ C
22∞∑

σmm=n17 / 28

Cubature for eigenfunctions（改）♣ 実は sup ノルムによる制約は外せる
Theorem
e1 , . . . , en−1 を exact に積分する n 点正重みカーネル求積公式
Qn であって∞∑wce(Qn )2 ≤ 4
σmm=nをみたすものが存在する∑
Meta-Algorithm において (d) の LP で ni=1 wi k(xi , xi ) を最小化するようにすれば、正の確率（P → 1 as N → ∞）で構成可能

18 / 28

Nyström 近似♣ サンプリングに基づく手法では Mercer 展開を使う▷ Nyström 近似で回避可能：z = (z1 , . . . , zℓ )⊤ に対し k(z, z) のs(≤ ℓ) ランクの最良近似を Ws (z) とし、その一般逆を W+s (z)kzs (x, y) ≈ k(x, z)W+s (z)k(z, y)
例えば次が示せる：
Theoremkz のなす積分作用素（in L2 (µ)）を Kz として、z ∼iid µ のとき
ss√()
2 supx k(x, x)1√
1 + 2 log
s+1 +δℓkKz − Kk≤ σsが確率 1 − δ 以上で成り立つ（k·k：作用素ノルム）▷ Nyström 近似で得られる近似固有関数をテスト関数に使用
19 / 28

数値実験：周期的 Sobolev 空間♣ [0, 1] 上の周期関数にノルム(∫ 1)2∫ 122rf(x) dx + (2π)f(r) (x)2 dx
kfk =00を入れる（r は正整数）と、ベルヌーイ多項式を用いてkr (x, y) = 1 +(−1)r−1 (2π)2r
B2r (|x − y|)(2r)!と再生核が表示できる [Wahba, 1990]µ：一様分布に対してσm ∼ m−r ,2−2rwce(Qoptn ) ∼nなどが知られている [Novak, 1988]▷ カーネル求積手法の簡単な評価に使用
20 / 28

数値実験：周期的 Sobolev 空間♣ 実験結果• それぞれの n で 50 回の平均をプロット• ‘Mercer’ では N = 10n, ‘Nyström’ では ℓ = N = 10n• ‘Uniform Grid’ が最適公式であることが既知
00.5

24

1.5

log10(wce)2log10(wce)2

1.0

2.0
2.5
3.08

MercerNyström
Monte Carloiid Bayes
Uniform Grid
0.8

1.06101.2

log10n1.4

1.6

Figure 1: r = 1

1.8

MercerNyström
Monte Carloiid Bayes
Uniform Grid
0.8

1.0

1.2

log10n

1.4

1.6

1.8

Figure 2: r = 321 / 28

数値実験：周期的 Sobolev 空間♣ 実験結果• それぞれの n で 50 回の平均をプロット• ‘Mercer’ では N = 10n, ‘Nyström’ では ℓ = N = 10n• ‘Uniform Grid’ が最適公式であることが既知• r = 3 でBelhadji et al. [2019] と比較
02

log10(wce)2

4
6810

MercerNyström
Monte Carloiid Bayes
Uniform Grid0.8(b) Sobolev space, d = 1, s = 3

1.0

1.2

log10n

1.4

1.6

1.8

22 / 28

数値実験：実データ♣ UCI Machine Learning Repository から Combined Cycle PowerPlant データセット [Kaya et al., 2012; Tüfekci, 2014] を使用• 9568 点からなる 5 次元データ → 等重み離散測度とみなす• パラメータを median heuristics で定めた Gauss カーネル
• 左から wce、X(5) の積分誤差、X(5) 1X(1) ,X(2) ≥0 の積分誤差

11.5

22

2.0

33

2.5

4
56N-reweight
Nyström
Monte Carloiid Bayes
herding
0.8

1.0

1.2

1.4 1.6
log10n

1.8

2.0

2.2

log10(error)2

1log10(error)2log10(wce)2• ℓ = N = 20n

4

3.0

5

3.5

6

4.0

70.8

1.0

1.2

1.4 1.6
log10n

1.8

2.0

2.2

4.5

0.8

1.0

1.2

1.4 1.6
log10n1.8

2.02.223 / 28

まとめ：スパース性から♣ RKHS に入っている関数を精度よく数値積分したい
理想的には、大きなサンプル x = (x1 , . . . , xN )⊤ からスパース最適化問題∫minimize w⊤ k(x, x)w − 2w⊤ X k(x, y) dµ(y)subject to w ≥ 0, 1⊤ w = 1, |w|0 ≤ n,を解いて「良い点配置」を得たい…これは難しいので• 貪欲法による最適化（herding, SBQ）
• そもそものサンプリングを工夫する（optimized, DPP, CVS）
• 良いテスト関数の cubature でサブサンプリング（+α）
24 / 28

まとめ：サンプリング +α

♣ まとめ
手法
OptimizedDPPCVS
提案手法
Nyström

wce(Qxn ,opt )2σ log n⌉∑⌈n/n ∞
m=n+1 σmσn+1∑∞
m=n σm？

Mercer
要要不要
要不要

サンプリング？
3O(n ) + 棄却 3
MCMC 近似？4iid + LP〃• 下二つは正重み• 経験的には DPP・提案手法・Nyström 共に O(σn+1 )34Hough et al. [2006,Algorithm
18] & Bardenet and Hardy [2020, Section 2.4]()
Mixing time O n5 log n [Rezaei and Gharan, 2019]25 / 28

References
Bach, F. (2017). On the equivalence between kernel quadrature rules andrandom feature expansions. The Journal of Machine Learning Research,
18(1):714–751.
Bardenet, R. and Hardy, A. (2020). Monte carlo with determinantal point
processes. The Annals of Applied Probability, 30(1):368–417.
Belhadji, A., Bardenet, R., and Chainais, P. (2019). Kernel quadrature withDPPs. In Advances in Neural Information Processing Systems,
volume 32, pages 12907–12917.
Belhadji, A., Bardenet, R., and Chainais, P. (2020). Kernel interpolation withcontinuous volume sampling. In International Conference on Machine
Learning, pages 725–735. PMLR.
Chen, Y., Welling, M., and Smola, A. (2010). Super-samples from kernel
herding. In Conference on Uncertainty in Artificial Intelligence, pages
109–116.
26 / 28

Hayakawa, S. (2021). Monte Carlo cubature construction. Japan Journal ofIndustrial and Applied Mathematics, 38:561–577.
Hayakawa, S., Oberhauser, H., and Lyons, T. (2021). Positively weighted kernel
quadrature via subsampling. arXiv preprint arXiv:2107.09597.
Hough, J. B., Krishnapur, M., Peres, Y., Virág, B., et al. (2006). Determinantalprocesses and independence. Probability surveys, 3:206–229.
Huszár, F. and Duvenaud, D. (2012). Optimally-weighted herding is Bayesian
quadrature. In Conference on Uncertainty in Artificial Intelligence,
pages 377–386.
Kaya, H., Tüfekci, P., and Gürgen, F. S. (2012). Local and global learning
methods for predicting power of a combined gas & steam turbine. InProceedings of the International Conference on Emerging Trends in
Computer and Electronics Engineering, pages 13–18.
Kiraly, F. J. and Oberhauser, H. (2019). Kernels for sequentially ordered data.
The Journal of Machine Learning Research, 20(31):1–45.
Litterer, C. and Lyons, T. (2012). High order recombination and an application
to cubature on Wiener space. The Annals of Applied Probability,22(4):1301–1327.
27 / 28

Lyons, T. and Victoir, N. (2004). Cubature on Wiener space. Proceedings of theRoyal Society of London Series A, 460:169–198.
Novak, E. (1988). Deterministic and stochastic error bounds in numericalanalysis. Springer.
Rezaei, A. and Gharan, S. O. (2019). A polynomial time MCMC method forsampling from continuous determinantal point processes. In
Chaudhuri, K. and Salakhutdinov, R., editors, Proceedings of the 36th
International Conference on Machine Learning, volume 97 ofProceedings of Machine Learning Research, pages 5438–5447. PMLR.
Tüfekci, P. (2014). Prediction of full load electrical power output of a baseload operated combined cycle power plant using machine learning
methods. International Journal of Electrical Power & Energy Systems,
60:126–140.
Wahba, G. (1990). Spline Models for Observational Data. Society for Industrial
and Applied Mathematics.

28 / 28