# 準 Monte Carlo 法：食い違い量・Koksma–Hlawka の不等式・格子則

> 独立乱数を決定的な点列に置き換える準 Monte Carlo 法を証明つきで扱う。Weyl の判定法、Koksma–Hlawka の不等式、van der Corput 列と Halton 列の食い違い量 O((log N)^d/N)、格子則の誤差と良い生成ベクトルの存在、無作為な平行移動、ANOVA 分解と有効次元。
> https://rikai.mugen-giken.com/computer-science/stochastic-computation/quasi-monte-carlo

## 0. この記事の要点

コンピュータで積分の近似値を求めるとき、乱数を使うと「たまたま偏った点が出る」ことで誤差が残ります。それなら、はじめから偏りが出ないように点を並べておけばよいのではないか、というのがこの記事の出発点です。点がどれだけ均等に並んでいるかを数で測り、その数が小さい並べ方を使うと、乱数より速く誤差が小さくなることを、証明つきで示します。読み終えると、どんな点の並べ方を選べばよいか、なぜそれで誤差が抑えられるか、実際の計算で誤差をどう見積もるかが分かります。

- 点集合が一様分布からどれだけずれているかを**星型食い違い量** $D_N^*$ で測ります。1 次元では $D_N^* \ge 1/(2N)$ で、閉じた式で計算できます（<Ref to="prop-discrepancy-1d" />）。
- 列が一様分布することは、すべての非零の整数ベクトル $h$ について指数和の平均 $N^{-1}\sum_n e^{2\pi i h\cdot x_n}$ が $0$ に収束することと同値です（<Ref to="thm-weyl-criterion" />）。
- 混合偏導関数が連続な被積分関数について、等重みの求積の誤差は「関数の変動 $\times$ 点集合の食い違い量」で抑えられます（<Ref to="thm-koksma-hlawka" />）。関数の性質と点の性質が分かれるのが要点です。
- Halton 列は $N D_N^* \le \prod_j((b_j - 1)\lceil\log_{b_j}N\rceil + 1) = O((\log N)^d)$ を満たします（<Ref to="thm-halton" />）。誤差は $O(N^{-1}(\log N)^d)$ で、Monte Carlo 法の $N^{-1/2}$ より速く減ります。
- 周期関数には**格子則**が向いています。誤差は双対格子（$h\cdot g$ が $N$ の倍数になる整数ベクトル $h$ の全体）の上の [Fourier 係数](/mathematics/fourier-analysis/fourier-series#def-fourier-coefficient)の和に等しく、素数 $N$ について最悪誤差が $O(N^{-\alpha}(\log N)^{\alpha d})$ になる生成ベクトル（格子の点 $\{ng/N\}$ を決める整数ベクトル $g$）があります（<Ref to="thm-lattice-rule-error" />）。
- 決定的な点集合に一様な無作為の平行移動を施すと、推定量が不偏になり、標準誤差を推定できます（<Ref to="prop-random-shift" />）。高次元で準 Monte Carlo 法が効くのは、被積分関数の分散が少数の変数の成分に集中している（有効次元が小さい。分散の大部分が何変数までの成分に収まるかを表す数）ときです（<Ref to="prop-effective-dimension" />）。

<Aside type="tip" title="この記事で使う事前知識">

- **Monte Carlo 積分の標準誤差**: 独立な乱数の平均で積分を近似したとき、誤差の大きさを表す量（$\sigma/\sqrt N$）。→ [標準誤差](/computer-science/stochastic-computation/monte-carlo-integration#def-mc-estimator)
- **Fourier 係数**: 周期関数を $e^{2\pi i h\cdot x}$ の重ね合わせで表したときの、各成分の係数。→ [Fourier 係数](/mathematics/fourier-analysis/fourier-series#def-fourier-coefficient)
- **Fubini の定理**: 多変数の積分を、1 変数ずつの積分の繰り返しに直して順序を入れ替えてよい、という定理。→ [直積測度と Fubini–Tonelli の定理](/mathematics/real-analysis/product-measures-fubini)
- **上限**: 集合の上界のうち最小のもの。食い違い量の定義で $\sup$ として使います。→ [上限](/mathematics/real-analysis/completeness-and-cauchy#def-supremum)
- **分散**: 確率変数が平均のまわりにどれだけ散らばるかを表す量。→ [分散](/mathematics/probability/random-variables#def-variance-covariance)

</Aside>

## 1. 動機

[Monte Carlo 積分](/computer-science/stochastic-computation/monte-carlo-integration)の誤差は $\sigma/\sqrt N$ です。この $N^{-1/2}$ は、独立な乱数が「偶然の偏り」をもつことから来ます。$N$ 個の独立な一様乱数は、ある区間に期待より $\sqrt N$ 個程度多く（または少なく）入ります。一方、積型の求積則は点を規則的に並べるので偏りはありませんが、$d$ 次元では 1 方向に $N^{1/d}$ 点しか置けず、誤差が $N^{-r/d}$ になります（<Ref to="computer-science/stochastic-computation/monte-carlo-integration#thm-mc-vs-quadrature" />）。

準 Monte Carlo 法（quasi-Monte Carlo）は、この中間を狙います。乱数の代わりに、どの箱にも体積にほぼ比例した個数が入るように設計した決定的な点列を使い、Monte Carlo 法と同じ等重みの平均 $N^{-1}\sum_n f(x_n)$ で積分を近似します。点の「むらのなさ」を測る量が食い違い量で、良い点列では $O(N^{-1}(\log N)^d)$ です。$d$ を固定して $N$ を増やせば $N^{-1/2}$ より速く、次元への依存は積型則の $N^{-r/d}$ よりずっと穏やかです。

代償は 2 つあります。第 1 に、誤差評価が被積分関数の滑らかさ（混合偏導関数）を要求すること。第 2 に、決定的な方法なので、そのままでは[標準誤差](/computer-science/stochastic-computation/monte-carlo-integration#def-mc-estimator)が付かないこと。後者は無作為化（§7）で解決します。このテーマの物差しである「推定量の偏りと分散を定理として評価し、報告する数値に誤差を付ける」ことは、準 Monte Carlo 法でも保てます。

## 2. 食い違い量

$d \ge 1$ とし、$[0,1)^d$ の $N$ 個の点 $x_1, \ldots, x_N$（重複があってもよい）の組を**点集合** $P$ と呼びます。$[0,1]^d$ の部分集合 $J$ について、$J$ に入る点の個数を $A(J; P) = \#\{n : x_n \in J\}$ と書きます。$t \in [0,1]^d$ について $[0, t) = \prod_{j=1}^d[0, t_j)$ とします。

<Definition id="def-discrepancy" title="局所食い違い関数・食い違い量">
点集合 $P$ の**局所食い違い関数**を
$$
\delta_P(t) = \frac{A([0,t); P)}{N} - \prod_{j=1}^dt_j \qquad (t \in [0,1]^d)
$$
と定め、
$$
D_N^*(P) = \sup_{t\in[0,1]^d}|\delta_P(t)|, \qquad D_N(P) = \sup_{0 \le a_j < b_j \le 1}\Bigl|\frac{A(\prod_j[a_j, b_j); P)}{N} - \prod_{j=1}^d(b_j - a_j)\Bigr|
$$
を**星型食い違い量**（star discrepancy）、**食い違い量**（extreme discrepancy）という。無限列 $(x_n)$ については、最初の $N$ 項のなす点集合の値を $D_N^*$、$D_N$ と書く。
</Definition>

$D_N^*$ は、原点を頂点とする箱について「入っている点の割合」と「体積」の差の最大値です。$N$ 個の独立な一様乱数では $D_N^*$ は $N^{-1/2}$ の大きさです。

<Lemma id="lem-extreme-vs-star" title="2 つの食い違い量の関係">
$D_N^*(P) \le D_N(P) \le 2^dD_N^*(P)$。
</Lemma>

<Proof of="lem-extreme-vs-star">
原点を頂点とする箱は $a = 0$ の箱なので $D_N^* \le D_N$ です（$t_j = 0$ の成分があるときは箱が空で、差は $0$ です）。逆向きを示します。箱 $J = \prod_j[a_j, b_j)$ の定義関数は
$$
\mathbf{1}_J(x) = \prod_{j=1}^d\bigl(\mathbf{1}_{[0,b_j)}(x_j) - \mathbf{1}_{[0,a_j)}(x_j)\bigr) = \sum_{v\subseteq\{1,\ldots,d\}}(-1)^{|v|}\mathbf{1}_{[0, c^v)}(x)
$$
と展開できます。ここで $c^v$ は、$j \in v$ で $c^v_j = a_j$、$j \notin v$ で $c^v_j = b_j$ とした点です。点の個数 $A(\cdot; P)$ も体積も、定義関数について線形なので、
$$
\frac{A(J; P)}{N} - \prod_j(b_j - a_j) = \sum_v(-1)^{|v|}\delta_P(c^v)
$$
です。右辺は $2^d$ 項で、各項の絶対値は $D_N^*$ 以下です。
</Proof>

1 次元では $D_N^*$ が閉じた形で求まります。

<Proposition id="prop-discrepancy-1d" title="1 次元の星型食い違い量">
$[0,1)$ の $N$ 個の点を小さい順に並べて $x_{(1)} \le \cdots \le x_{(N)}$ とすると
$$
D_N^* = \frac{1}{2N} + \max_{1\le i\le N}\Bigl|x_{(i)} - \frac{2i-1}{2N}\Bigr|.
$$
特に $D_N^* \ge 1/(2N)$ で、等号は $x_{(i)} = (2i-1)/(2N)$（区間を $N$ 等分したときの中点）のときに限る。
</Proposition>

<Proof of="prop-discrepancy-1d">
$x_{(0)} = 0$、$x_{(N+1)} = 1$ とおきます。$[0, t)$ に入る点は $x_n < t$ を満たすものなので、$x_{(i)} < t \le x_{(i+1)}$ のとき $A([0,t)) = i$ です（$i = 0, \ldots, N$。$t = 0$ では $A = 0$）。この範囲で $\delta_P(t) = i/N - t$ は $t$ の 1 次関数なので、$|\delta_P|$ の上限は区間の両端での値 $|i/N - x_{(i)}|$、$|i/N - x_{(i+1)}|$ の大きいほうです（左端は区間に含まれませんが、上限は変わりません）。$x_{(i)} = x_{(i+1)}$ で区間が空のときは寄与がありません。

まず点がすべて異なる場合を考えます。$i = 0$ の区間からは $|0 - x_{(1)}|$、$i = N$ の区間からは $|1 - x_{(N)}|$ が出るので、まとめると
$$
D_N^* = \max_{1\le i\le N}\max\Bigl(\Bigl|x_{(i)} - \frac iN\Bigr|,\ \Bigl|x_{(i)} - \frac{i-1}{N}\Bigr|\Bigr)
$$
です（$x_{(i)}$ は、$i$ 番目の区間の左端と $i-1$ 番目の区間の右端として 2 回現れます）。実数 $x$ と $\alpha < \beta$ について $\max(|x - \alpha|, |x - \beta|) = |x - (\alpha+\beta)/2| + (\beta - \alpha)/2$ なので、$\alpha = (i-1)/N$、$\beta = i/N$ として主張の式になります。

点に重複がある場合。$x_{(p)} < x_{(p+1)} = \cdots = x_{(q)} < x_{(q+1)}$ とし、共通の値を $x$ とします。上の議論で実際に現れるのは $|x - p/N|$（$p$ 番目の区間の右端）と $|x - q/N|$（$q$ 番目の区間の左端）です。$p < i \le q$ について $(i-1)/N$ と $i/N$ は $p/N$ と $q/N$ の間にあるので、$|x - (i-1)/N|$ と $|x - i/N|$ はこの 2 つの大きいほう以下です。したがって、空の区間に対応する項を最大に加えても値は変わらず、同じ式が成り立ちます。

最後の主張は、第 2 項が非負で、$0$ になるのがすべての $i$ で $x_{(i)} = (2i-1)/(2N)$ のときだけであることから従います。
</Proof>

$N$ 点の集合としては $1/(2N)$ が限界です。ただし中点の集合は $N$ を決めないと作れません。$N$ を増やしながら点を足していける無限列で、どの $N$ でも $D_N^*$ が小さいものを作るのが §5 の主題です。

## 3. 一様分布列と Weyl の判定法

<Definition id="def-uniformly-distributed" title="一様分布列">
$[0,1)^d$ の点列 $(x_n)_{n\ge1}$ が**一様分布する**（uniformly distributed modulo 1）とは、すべての箱 $J = \prod_j[a_j, b_j)$（$0 \le a_j < b_j \le 1$）について
$$
\lim_{N\to\infty}\frac{A(J; x_1, \ldots, x_N)}{N} = \prod_{j=1}^d(b_j - a_j)
$$
が成り立つことをいう。
</Definition>

箱ごとの条件を、関数についての条件に言い換えます。$h \in \mathbb{Z}^d$ と $x \in \mathbb{R}^d$ について $h\cdot x = \sum_jh_jx_j$、$e(s) = e^{2\pi is}$ と書きます。

<Theorem id="thm-weyl-criterion" title="Weyl の判定法">
$[0,1)^d$ の点列 $(x_n)_{n\ge1}$ について、次は同値である。

1. $(x_n)$ は一様分布する。
2. すべての連続関数 $f : [0,1]^d \to \mathbb{C}$ について $\displaystyle\lim_{N\to\infty}\frac1N\sum_{n=1}^Nf(x_n) = \int_{[0,1]^d}f(x)\,dx$。
3. すべての $h \in \mathbb{Z}^d \setminus \{0\}$ について $\displaystyle\lim_{N\to\infty}\frac1N\sum_{n=1}^Ne(h\cdot x_n) = 0$。
</Theorem>

証明は 1 ⇒ 2 ⇒ 3 ⇒ 1 の順に進みます。1 ⇒ 2 は、連続関数を階段関数ではさんで箱の条件に帰着させます。2 ⇒ 3 は指数関数が連続であることを使うだけです。3 ⇒ 1 では、指数関数の有限和（三角多項式）で周期的な連続関数を近似し、さらに箱を周期的な連続関数ではさみます。一番長いのは 3 ⇒ 1 です。

<Proof of="thm-weyl-criterion">
$M_N(f) = N^{-1}\sum_{n\le N}f(x_n)$、$I(f) = \int_{[0,1]^d}f$ と書きます。

**1 ⇒ 2.** 実部と虚部に分ければよいので、$f$ は実数値とします。箱の定義関数の有限な一次結合（階段関数）$s$ については、1 と線形性から $M_N(s) \to I(s)$ です。$\varepsilon > 0$ を取ります。$f$ はコンパクト集合 $[0,1]^d$ 上で一様連続なので（<Ref to="mathematics/real-analysis/arzela-ascoli-stone-weierstrass#lem-uniform-continuity-compact" />）、自然数 $m$ を大きく取れば、1 辺 $1/m$ の各立方体 $J_c = \prod_j[c_j/m, (c_j+1)/m)$（$c \in \{0, \ldots, m-1\}^d$）の閉包の上で $f$ の上限と下限の差が $\varepsilon$ 以下になります。$J_c$ の上で、その閉包での $f$ の上限（下限）を値とする階段関数を $s_+$（$s_-$）とすると、$[0,1)^d$ 上で $s_- \le f \le s_+$、$s_+ - s_- \le \varepsilon$ です。点はすべて $[0,1)^d$ にあるので
$$
\limsup_N M_N(f) \le \lim_N M_N(s_+) = I(s_+) \le I(f) + \varepsilon, \qquad \liminf_N M_N(f) \ge I(s_-) \ge I(f) - \varepsilon
$$
です。$\varepsilon$ は任意なので $M_N(f) \to I(f)$ です。

**2 ⇒ 3.** $f(x) = e(h\cdot x)$ は連続で、$h \ne 0$ なら、$h_j \ne 0$ である座標 $x_j$ について先に積分すると $\int_0^1e(h_jx_j)\,dx_j = 0$ なので $I(f) = 0$ です。

**3 ⇒ 1.** まず、$\mathbb{R}^d$ 上の連続関数 $g$ で各変数について周期 $1$ をもつものについて $M_N(g) \to I(g)$ を示します。$K = \{z \in \mathbb{C}^d : |z_j| = 1\}$ はコンパクト Hausdorff 空間で、$g(x) = G(e(x_1), \ldots, e(x_d))$ となる $K$ 上の関数 $G$ が定まります。$x \mapsto (e(x_j))_j$ はコンパクト集合 $[0,1]^d$ から $K$ への連続な全射なので閉写像で、$G$ は連続です。$z^h = \prod_jz_j^{h_j}$（$h \in \mathbb{Z}^d$）の有限な一次結合の全体は、積で閉じ、定数 $1 = z^0$ を含み、$\overline{z^h} = z^{-h}$ から複素共役で閉じ、座標関数 $z_j$ を含むので点を分離します。Stone–Weierstrass の定理（<Ref to="mathematics/real-analysis/arzela-ascoli-stone-weierstrass#thm-stone-weierstrass-complex" />）から、$\varepsilon > 0$ に対して三角多項式 $p(x) = \sum_{h\in H}c_he(h\cdot x)$（$H$ は有限集合）で $\sup_x|g(x) - p(x)| \le \varepsilon$ となるものがあります。3 から、$h \ne 0$ の各項で $M_N(e(h\cdot{})) \to 0 = I(e(h\cdot{}))$、$h = 0$ の項では $M_N = I = c_0$ なので $M_N(p) \to I(p)$ です。$|M_N(g) - M_N(p)| \le \varepsilon$、$|I(g) - I(p)| \le \varepsilon$ から $\limsup_N|M_N(g) - I(g)| \le 2\varepsilon$ で、$\varepsilon$ は任意なので $M_N(g) \to I(g)$ です。

次に箱 $J = \prod_j[a_j, b_j)$ を周期的な連続関数ではさみます。円周 $\mathbb{R}/\mathbb{Z}$ 上の距離を $\mathrm{dist}$ と書き、$0 < \eta < 1/2$ について
$$
\psi_j^+(x) = \max\Bigl(0,\ 1 - \frac{\mathrm{dist}(x, [a_j, b_j])}{\eta}\Bigr), \qquad \psi_j^-(x) = \min\Bigl(1,\ \frac{\mathrm{dist}(x, (\mathbb{R}/\mathbb{Z})\setminus(a_j, b_j))}{\eta}\Bigr)
$$
とおきます。どちらも周期 $1$ の連続関数で、$0 \le \psi_j^- \le \mathbf{1}_{[a_j, b_j)} \le \psi_j^+ \le 1$ です。$\psi_j^+$ は $[a_j, b_j]$ から距離 $\eta$ 以上で $0$、$\psi_j^-$ は $(a_j, b_j)$ の端から距離 $\eta$ 以上内側で $1$ なので、$\int_0^1(\psi_j^+ - \psi_j^-) \le 4\eta$ です（両端の幅 $\eta$ ずつの範囲で差が $1$ 以下）。$g^\pm(x) = \prod_j\psi_j^\pm(x_j)$ とおくと $g^- \le \mathbf{1}_J \le g^+$ で、$[0,1]$ の数 $u_j, v_j$ について $|\prod_ju_j - \prod_jv_j| \le \sum_j|u_j - v_j|$ が成り立つこと（1 因子ずつ取り替える）から
$$
I(g^+) - I(g^-) = \prod_j\int_0^1\psi_j^+ - \prod_j\int_0^1\psi_j^- \le 4d\eta
$$
です。前半から $M_N(g^\pm) \to I(g^\pm)$ なので
$$
\limsup_N\frac{A(J)}{N} \le I(g^+) \le I(g^-) + 4d\eta \le \prod_j(b_j - a_j) + 4d\eta, \qquad \liminf_N\frac{A(J)}{N} \ge I(g^-) \ge \prod_j(b_j - a_j) - 4d\eta
$$
です（$I(g^-) \le \mathrm{vol}(J) \le I(g^+)$ を使いました）。$\eta \to 0$ として 1 を得ます。
</Proof>

この判定法の利点は、箱という扱いにくい集合の代わりに、指数関数という計算しやすい関数を調べればよいことです。最初の応用が次の系です。

<Corollary id="cor-kronecker" title="Kronecker 列の一様分布">
実数 $\alpha_1, \ldots, \alpha_d$ について、$1, \alpha_1, \ldots, \alpha_d$ が有理数体上一次独立ならば、$x_n = (\{n\alpha_1\}, \ldots, \{n\alpha_d\})$（$\{\cdot\}$ は小数部分）は一様分布する。特に $\alpha$ が無理数なら $(\{n\alpha\})$ は $[0,1)$ で一様分布する。
</Corollary>

<Proof of="cor-kronecker">
$h \in \mathbb{Z}^d \setminus \{0\}$ とし、$\theta = h\cdot\alpha = \sum_jh_j\alpha_j$ とおきます。$\theta$ が整数 $k$ なら $k\cdot1 - \sum_jh_j\alpha_j = 0$ が自明でない有理係数の一次関係になり、仮定に反します。よって $\theta \notin \mathbb{Z}$ で $e(\theta) \ne 1$ です。$e(h\cdot x_n) = e(n\theta)$（小数部分を取っても整数の差しかないため）なので、等比数列の和から
$$
\Bigl|\frac1N\sum_{n=1}^Ne(n\theta)\Bigr| = \frac1N\Bigl|\frac{e(\theta)(e(N\theta) - 1)}{e(\theta) - 1}\Bigr| \le \frac{2}{N|e(\theta) - 1|} = \frac{1}{N|\sin\pi\theta|} \to 0
$$
です。<Ref to="thm-weyl-criterion" /> の 3 が成り立つので、一様分布します。
</Proof>

<Example id="ex-golden-kronecker" title="黄金比の Kronecker 列">
$\alpha = (\sqrt5 - 1)/2 = 0.618\ldots$ は無理数なので、$x_n = \{n\alpha\}$ は一様分布します。最初の数項は $0.618, 0.236, 0.854, 0.472, 0.090, 0.708, 0.326, 0.944$ で、新しい点が既にある点の間の広い隙間に入っていきます。無理数回転がエルゴード的であること（<Ref to="mathematics/probability/ergodic-theorems#thm-rotation-ergodic" />）は「ほとんどすべての出発点で時間平均が空間平均に一致する」ことを言いますが、<Ref to="cor-kronecker" /> は出発点を $0$ に固定した列についての主張で、例外集合がありません。黄金比は連分数展開の部分商がすべて $1$ で、有理数で最も近似しにくい無理数です。この列は $D_N^* = O(N^{-1}\log N)$ を満たすことが、連分数を使って証明できます（Kuipers–Niederreiter の第 2 章）。
</Example>

一様分布は $N \to \infty$ の定性的な性質です。数値計算に必要なのは、有限の $N$ での誤差の大きさです。次節で、誤差を $D_N^*$ で抑えます。

## 4. Koksma–Hlawka の不等式

### 4.1. 1 次元

<Theorem id="thm-koksma" title="Koksma の不等式（C¹ 級の場合）">
$f : [0,1] \to \mathbb{R}$ を $C^1$ 級、$P = \{x_1, \ldots, x_N\} \subset [0,1)$ とする。このとき
$$
\frac1N\sum_{n=1}^Nf(x_n) - \int_0^1f(x)\,dx = -\int_0^1\delta_P(t)f'(t)\,dt
$$
（Zaremba の恒等式）が成り立ち、したがって
$$
\Bigl|\frac1N\sum_{n=1}^Nf(x_n) - \int_0^1f(x)\,dx\Bigr| \le D_N^*(P)\int_0^1|f'(t)|\,dt.
$$
</Theorem>

<Proof of="thm-koksma">
微分積分学の基本定理（<Ref to="mathematics/calculus/integration-and-ftc#thm-ftc2" />）から、各点 $x_n$ について
$$
f(x_n) = f(1) - \int_{x_n}^1f'(t)\,dt = f(1) - \int_0^1\mathbf{1}_{[0,t)}(x_n)f'(t)\,dt
$$
です（$\mathbf{1}_{[0,t)}(x_n) = 1$ は $x_n < t$ と同値）。$n$ について平均すると、$N^{-1}\sum_n\mathbf{1}_{[0,t)}(x_n) = A([0,t); P)/N$ なので
$$
\frac1N\sum_nf(x_n) = f(1) - \int_0^1\frac{A([0,t); P)}{N}f'(t)\,dt
$$
です。一方、部分積分（<Ref to="mathematics/calculus/integration-and-ftc#thm-by-parts" />）を $u(t) = f(t)$、$v(t) = t$ で使うと $\int_0^1f(t)\,dt = f(1) - \int_0^1tf'(t)\,dt$ です。差を取ると $-\int_0^1(A([0,t);P)/N - t)f'(t)\,dt = -\int_0^1\delta_P(t)f'(t)\,dt$ で、恒等式を得ます。不等式は $|\delta_P(t)| \le D_N^*(P)$ から従います。
</Proof>

誤差の上界は、関数だけで決まる量 $\int|f'|$（$f$ の全変動）と、点だけで決まる量 $D_N^*$ の積です。<Ref to="prop-discrepancy-1d" /> と合わせると、中点則の誤差が $\int|f'|/(2N)$ 以下であることになります。

### 4.2. 多次元

多次元では、導関数の役を混合偏導関数が果たします。$u \subseteq \{1, \ldots, d\}$ について、$u$ に属する各変数で 1 回ずつ偏微分したものを $\partial_uf = \bigl(\prod_{j\in u}\partial/\partial x_j\bigr)f$ と書きます（$\partial_\emptyset f = f$）。$f$ が $[0,1]^d$ を含む開集合の上で $C^d$ 級なら、$\partial_uf$ は微分の順序によらず定まり連続です（<Ref to="mathematics/calculus/multivariable-differentiation#thm-schwarz-n" /> を繰り返し使います）。$t_u = (t_j)_{j\in u} \in [0,1]^u$ について、$u$ の座標を $t_j$、それ以外の座標を $1$ とした $[0,1]^d$ の点を $(t_u, 1)$ と書きます。

次の定義の直観を $d = 2$ で述べます。このとき $u$ は $\{1\}$、$\{2\}$、$\{1,2\}$ の 3 通りで、変動は
$\int_0^1|\partial_1f(t_1,1)|\,dt_1 + \int_0^1|\partial_2f(1,t_2)|\,dt_2 + \iint|\partial_1\partial_2f(t_1,t_2)|\,dt_1dt_2$
です。はじめの 2 項は正方形の辺（$x_2 = 1$ の辺と $x_1 = 1$ の辺）に沿った関数の変化の大きさ、最後の項は正方形の内部での「ねじれ」の大きさです。1 次元の全変動の多次元版だと考えてください。

<Definition id="def-hardy-krause" title="Hardy–Krause の意味の変動（滑らかな関数の場合）">
$[0,1]^d$ を含む開集合の上で $C^d$ 級の関数 $f$ について
$$
V_{\mathrm{HK}}(f) = \sum_{\emptyset\ne u\subseteq\{1,\ldots,d\}}\int_{[0,1]^u}\bigl|\partial_uf(t_u, 1)\bigr|\,dt_u
$$
を $f$ の **Hardy–Krause の意味の変動**という。
</Definition>

和は $2^d - 1$ 項です。$u$ の項は、$f$ を面 $\{x_j = 1\ (j \notin u)\}$ に制限した $|u|$ 変数の関数の、最高階の混合偏導関数の $L^1$ ノルムです。$d = 1$ では $\int_0^1|f'|$ です。

<Lemma id="lem-anchored-expansion" title="頂点 (1, …, 1) を基点とする展開">
$f$ を $[0,1]^d$ を含む開集合の上で $C^d$ 級の関数とする。すべての $x \in [0,1]^d$ について
$$
f(x) = \sum_{u\subseteq\{1,\ldots,d\}}(-1)^{|u|}\int_{[0,1]^u}\prod_{j\in u}\mathbf{1}_{[0,t_j)}(x_j)\,\partial_uf(t_u, 1)\,dt_u
$$
が成り立つ（$u = \emptyset$ の項は $f(1, \ldots, 1)$）。
</Lemma>

<Proof of="lem-anchored-expansion">
$d$ についての帰納法で示します。$d = 1$ は <Ref to="thm-koksma" /> の証明のはじめの式 $f(x) = f(1) - \int_0^1\mathbf{1}_{[0,t)}(x)f'(t)\,dt$ です。

$d - 1$ 変数で成り立つとし、$x = (x', x_d)$ と書きます。$x'$ を固定して $x_d$ について $d = 1$ の場合を使うと
$$
f(x', x_d) = f(x', 1) - \int_0^1\mathbf{1}_{[0,t_d)}(x_d)\,\partial_df(x', t_d)\,dt_d
$$
です。$x' \mapsto f(x', 1)$ と、各 $t_d$ についての $x' \mapsto \partial_df(x', t_d)$ は、$[0,1]^{d-1}$ を含む開集合の上で $C^{d-1}$ 級です。帰納法の仮定をこの 2 つに使います。$u' \subseteq \{1, \ldots, d-1\}$ について、前者の $\partial_{u'}$ を $(t_{u'}, 1)$ で取った値は $\partial_{u'}f(t_{u'}, 1)$（$d$ 番目の座標も $1$）、後者は $\partial_{u'\cup\{d\}}f$ を「$u'$ の座標が $t_j$、$d$ 番目が $t_d$、残りが $1$」の点で取った値です。したがって
$$
f(x) = \sum_{u'}(-1)^{|u'|}\int\prod_{j\in u'}\mathbf{1}_{[0,t_j)}(x_j)\,\partial_{u'}f(t_{u'}, 1)\,dt_{u'} + \sum_{u'}(-1)^{|u'|+1}\int_0^1\!\!\int\prod_{j\in u'\cup\{d\}}\mathbf{1}_{[0,t_j)}(x_j)\,\partial_{u'\cup\{d\}}f(t_{u'\cup\{d\}}, 1)\,dt_{u'}\,dt_d
$$
です。第 2 の和の被積分関数は立方体の上の有界な可測関数なので、累次積分は $[0,1]^{u'\cup\{d\}}$ 上の重積分に等しくなります（<Ref to="mathematics/real-analysis/product-measures-fubini#thm-fubini" />）。第 1 の和は $d \notin u$ である $u = u'$ の項、第 2 の和は $d \in u$ である $u = u'\cup\{d\}$ の項で、合わせて $\{1, \ldots, d\}$ のすべての部分集合にわたる和になります。
</Proof>

<Theorem id="thm-koksma-hlawka" title="Koksma–Hlawka の不等式">
$f$ を $[0,1]^d$ を含む開集合の上で $C^d$ 級の実数値関数、$P = \{x_1, \ldots, x_N\} \subset [0,1)^d$ とする。このとき
$$
\frac1N\sum_{n=1}^Nf(x_n) - \int_{[0,1]^d}f(x)\,dx = \sum_{\emptyset\ne u\subseteq\{1,\ldots,d\}}(-1)^{|u|}\int_{[0,1]^u}\delta_P(t_u, 1)\,\partial_uf(t_u, 1)\,dt_u
$$
（Hlawka–Zaremba の恒等式）が成り立ち、したがって
$$
\Bigl|\frac1N\sum_{n=1}^Nf(x_n) - \int_{[0,1]^d}f(x)\,dx\Bigr| \le V_{\mathrm{HK}}(f)\,D_N^*(P).
$$
</Theorem>

<Proof of="thm-koksma-hlawka">
<Ref to="lem-anchored-expansion" /> の式で $x = x_n$ とし、$n$ について平均します。$x_n \in [0,1)^d$ なので、$j \notin u$ の座標は自動的に $x_{n,j} < 1$ を満たし、
$$
\frac1N\sum_{n=1}^N\prod_{j\in u}\mathbf{1}_{[0,t_j)}(x_{n,j}) = \frac{A([0, (t_u, 1)); P)}{N}
$$
です。次に同じ式を $x \in [0,1]^d$ について積分します。右辺の各項は有界な可測関数の積分なので、Fubini の定理（<Ref to="mathematics/real-analysis/product-measures-fubini#thm-fubini" />）で $x$ についての積分を先に行うと、$\int_{[0,1]^d}\prod_{j\in u}\mathbf{1}_{[0,t_j)}(x_j)\,dx = \prod_{j\in u}t_j$ で、これは箱 $[0, (t_u, 1))$ の体積です。2 つの式の差を取ると、$u = \emptyset$ の項 $f(1, \ldots, 1)$ は打ち消し合い、$u \ne \emptyset$ の項には
$$
\frac{A([0,(t_u,1));P)}{N} - \prod_{j\in u}t_j = \delta_P(t_u, 1)
$$
が現れて、恒等式を得ます。不等式は $|\delta_P(t_u, 1)| \le D_N^*(P)$ と <Ref to="def-hardy-krause" /> から従います。
</Proof>

<Remark id="rem-bounded-variation" title="有界変動関数への拡張と、不等式の読み方">
不等式は、Hardy–Krause の意味で有界変動な関数（滑らかとは限らない）について同じ形で成り立ちます。そのときの $V_{\mathrm{HK}}$ は、箱への分割についての交代和の上限（Vitali の変動）を各面について足したものです。この一般の場合の証明は Riemann–Stieltjes 積分についての部分積分によるもので、Kuipers–Niederreiter の第 2 章にあります。この記事では $C^d$ 級の場合だけを定理にしました。

実用上の注意が 2 つあります。第 1 に、座標軸に平行でない境界をもつ集合の定義関数（$d \ge 2$）や、$\max(x_1 + x_2 + x_3 - 1, 0)$ のような座標軸に平行でない折れ目のある関数（$d \ge 3$）は $V_{\mathrm{HK}} = \infty$ になり、この不等式は何も言いません（準 Monte Carlo 法が使えないという意味ではなく、保証がないという意味です）。折れ目のある関数では、1 辺 $h$ の小さな立方体のうち折れ目と交わる約 $h^{-(d-1)}$ 個のそれぞれで交代和が $h$ の大きさになるので、変動は $h^{-(d-2)}$ の大きさで、$d = 2$ では有限にとどまり、$d \ge 3$ で発散します。第 2 に、上界は最悪の場合の評価で、個々の $f$ では大きく過大です。$V_{\mathrm{HK}}$ も $D_N^*$ も計算が難しいので、この不等式は誤差の見積もりには使わず、収束の次数の根拠として使います。実際の誤差の見積もりは §7 の無作為化で行います。
</Remark>

## 5. 食い違い量の小さい列

### 5.1. van der Corput 列

<Definition id="def-radical-inverse" title="逆進関数・van der Corput 列・Halton 列・Hammersley 点集合">
$b \ge 2$ を整数とする。非負整数 $n$ の $b$ 進展開を $n = \sum_{k\ge0}a_k(n)b^k$（$a_k(n) \in \{0, \ldots, b-1\}$、有限個を除いて $0$）とするとき
$$
\varphi_b(n) = \sum_{k\ge0}a_k(n)\,b^{-k-1} \in [0, 1)
$$
を基数 $b$ の**逆進関数**（radical inverse）という。$b$ 進展開の桁を小数点について折り返したものである。

- 列 $(\varphi_b(n))_{n\ge0}$ を基数 $b$ の **van der Corput 列**という。
- 整数 $b_1, \ldots, b_d \ge 2$ について、$[0,1)^d$ の列 $x_n = (\varphi_{b_1}(n), \ldots, \varphi_{b_d}(n))$（$n \ge 0$）を **Halton 列**という。
- $N \ge 1$ について、$[0,1)^d$ の $N$ 点 $(n/N, \varphi_{b_1}(n), \ldots, \varphi_{b_{d-1}}(n))$（$n = 0, \ldots, N-1$）を **Hammersley 点集合**という。

以下、これらの列の「最初の $N$ 項」は $n = 0, \ldots, N-1$ の項を指す。
</Definition>

基数 $2$ では $0, \tfrac12, \tfrac14, \tfrac34, \tfrac18, \tfrac58, \tfrac38, \tfrac78, \ldots$ です。区間を半分、そのまた半分と、粗いほうから順に埋めていきます。この「どの段階で止めても均等」という性質を、次の補題が表します。$c\,b^{-m}$ を左端とする長さ $b^{-m}$ の区間 $[c\,b^{-m}, (c+1)b^{-m})$（$0 \le c < b^m$）を、基数 $b$ の**階数 $m$ の基本区間**と呼びます。

<Lemma id="lem-elementary-interval" title="基本区間に入る条件は合同式">
$b \ge 2$、$m \ge 0$ とする。$\{0, \ldots, b^m - 1\}$ からそれ自身への全単射 $c \mapsto r(c)$ があって、すべての $n \ge 0$ と $0 \le c < b^m$ について
$$
\varphi_b(n) \in [c\,b^{-m}, (c+1)b^{-m}) \iff n \equiv r(c) \pmod{b^m}
$$
が成り立つ。
</Lemma>

<Proof of="lem-elementary-interval">
$n = r + b^mq$（$0 \le r < b^m$、$q \ge 0$）と書きます。$n$ の $b$ 進展開の下位 $m$ 桁は $r$ の桁 $a_0, \ldots, a_{m-1}$、それより上は $q$ の桁なので
$$
\varphi_b(n) = \sum_{k=0}^{m-1}a_kb^{-k-1} + b^{-m}\varphi_b(q) = \frac{\rho(r)}{b^m} + \frac{\varphi_b(q)}{b^m}, \qquad \rho(r) = \sum_{k=0}^{m-1}a_kb^{m-1-k}
$$
です。$\rho(r)$ は $r$ の $m$ 桁を逆順に読んだ整数で、$\rho$ は $\{0, \ldots, b^m-1\}$ の全単射です（2 回行うと元に戻ります）。$\varphi_b(q)$ は有限個の桁の和なので $0 \le \varphi_b(q) < \sum_{k\ge0}(b-1)b^{-k-1} = 1$ です。したがって $\varphi_b(n) \in [\rho(r)b^{-m}, (\rho(r)+1)b^{-m})$ で、$\varphi_b(n)$ が $c$ 番目の基本区間に入ることは $\rho(r) = c$、すなわち $n \equiv \rho^{-1}(c) \pmod{b^m}$ と同値です。$r(c) = \rho^{-1}(c)$ とすればよい。
</Proof>

<Theorem id="thm-van-der-corput" title="van der Corput 列の食い違い量">
基数 $b \ge 2$ の van der Corput 列の最初の $N$ 項（$N \ge 1$）について
$$
N\,D_N^* \le (b-1)\bigl(\lfloor\log_bN\rfloor + 1\bigr) \le (b-1)\frac{\log N}{\log b} + b - 1.
$$
</Theorem>

<Proof of="thm-van-der-corput">
**ブロックの評価。** $m \ge 0$、$k \ge 0$ とし、連続する $b^m$ 個の添字 $n = kb^m, \ldots, (k+1)b^m - 1$ の項からなる点集合を $Q$ とします。この添字は $b^m$ を法とする各剰余類をちょうど 1 回ずつ通るので、<Ref to="lem-elementary-interval" /> から、階数 $m$ の各基本区間に $Q$ の点がちょうど 1 個ずつ入ります。$t \in [0,1]$ とし $c = \lfloor tb^m\rfloor$ とおくと、$[0, t)$ は $c$ 個の基本区間 $[0, cb^{-m})$ を含み、$c+1$ 個の基本区間 $[0, (c+1)b^{-m})$（$c = b^m$ のときは $[0,1)$）に含まれます。よって $c \le A([0,t); Q) \le c+1$ で、$c \le tb^m < c + 1$ と合わせて
$$
\bigl|A([0,t); Q) - b^mt\bigr| \le 1
$$
です。

**分割。** $L = \lfloor\log_bN\rfloor$ とし、$N = \sum_{m=0}^La_mb^m$ を $b$ 進展開とします。添字の集合 $\{0, \ldots, N-1\}$ を、先頭から順に、長さ $b^L$ のブロック $a_L$ 個、長さ $b^{L-1}$ のブロック $a_{L-1}$ 個、…、長さ $1$ のブロック $a_0$ 個に分けます。長さ $b^m$ のブロックが始まる位置は $\sum_{i>m}a_ib^i + jb^m$（$0 \le j < a_m$）で、$b^m$ の倍数です。したがって各ブロックは上の $Q$ の形で、ブロックごとに $|A([0,t); Q) - b^mt| \le 1$ です。全部で $\sum_ma_m \le (b-1)(L+1)$ 個のブロックについて足すと
$$
\bigl|A([0,t); P_N) - Nt\bigr| \le (b-1)(L+1)
$$
がすべての $t$ で成り立ちます。$N$ で割って上限を取れば主張の式です。第 2 の不等式は $\lfloor\log_bN\rfloor \le \log N/\log b$ です。
</Proof>

独立な乱数の $D_N^* \sim N^{-1/2}$ に対し、$O(N^{-1}\log N)$ です。<Ref to="thm-koksma" /> と合わせると、$C^1$ 級の関数の積分誤差が $O(N^{-1}\log N)$ になります。

### 5.2. Halton 列

多次元では、座標ごとに異なる基数を使います。基数が互いに素であることが、座標どうしの「独立性」の代わりをします。

<Theorem id="thm-halton" title="Halton 列の食い違い量">
$b_1, \ldots, b_d \ge 2$ をどの 2 つも互いに素な整数とする。基数 $b_1, \ldots, b_d$ の Halton 列の最初の $N$ 項（$N \ge 2$）について
$$
N\,D_N^* \le \prod_{j=1}^d\Bigl((b_j - 1)\lceil\log_{b_j}N\rceil + 1\Bigr).
$$
特に $D_N^* = O(N^{-1}(\log N)^d)$ である。
</Theorem>

証明の方針は van der Corput 列と同じ数え方です。$[0, t)$ を、基本区間の積である箱と、はみ出す小さな箱に分け、どの箱でも点の個数が体積の $N$ 倍から高々 $1$ しかずれないことを示します。箱の数が上界の右辺の積になります。基数が互いに素であることは、各座標の条件をまとめて 1 つの合同式に直す中国剰余定理を使うために要ります。

<Proof of="thm-halton">
最初の $N$ 項を $P_N$ と書きます。

**基本的な箱。** 各 $j$ について基数 $b_j$ の階数 $m_j$ の基本区間 $E_j$ を取り、$E = \prod_jE_j$、$Q = \prod_jb_j^{m_j}$ とします（$\mathrm{vol}(E) = 1/Q$）。<Ref to="lem-elementary-interval" /> から、$x_n \in E$ は連立合同式 $n \equiv r_j \pmod{b_j^{m_j}}$（$j = 1, \ldots, d$）と同値です。$b_j$ がどの 2 つも互いに素なので、法 $b_j^{m_j}$ もそうです。中国剰余定理（<Ref to="mathematics/number-theory/elementary-number-theory#cor-crt-multi" />）から、この連立合同式は $Q$ を法としてただ 1 つの解をもちます。したがって、連続する $Q$ 個の整数の中に $x_n \in E$ となる $n$ がちょうど 1 個あり、$\{0, \ldots, N-1\}$ の中には $\lfloor N/Q\rfloor$ 個または $\lceil N/Q\rceil$ 個あります。よって
$$
\bigl|A(E; P_N) - N\,\mathrm{vol}(E)\bigr| \le 1, \qquad Q \ge N \text{ ならば } A(E; P_N) \le 1.
$$

**区間の分解。** $t \in [0,1]^d$ とし、$L_j = \lceil\log_{b_j}N\rceil$（$b_j^{L_j} \ge N$）とおきます。各 $j$ で、$c_j = \lfloor t_jb_j^{L_j}\rfloor$、$s_j = c_jb_j^{-L_j} \le t_j$ とします。$c_j$ を $c_j = \sum_{k=0}^{L_j}\tau_kb_j^{L_j-k}$（$\tau_0 \in \{0, 1\}$、$k \ge 1$ で $\tau_k \in \{0, \ldots, b_j-1\}$。$\tau_0 = 1$ は $c_j = b_j^{L_j}$ のときだけで、そのとき他の桁は $0$）と展開すると、$[0, s_j)$ は左から順に、階数 $k$ の基本区間 $\tau_k$ 個（$k = 0, 1, \ldots, L_j$）の交わりのない和に分かれます。階数 $k$ の区間の左端は $\sum_{i<k}\tau_ib_j^{-i} + ab_j^{-k}$（$0 \le a < \tau_k$）で $b_j^{-k}$ の倍数なので、確かに基本区間です。個数は $(b_j - 1)L_j$ 以下です（$\tau_0 = 1$ のときは $1$ 個）。残りの $[s_j, t_j)$ は長さが $b_j^{-L_j}$ 未満で、階数 $L_j$ の 1 つの基本区間 $R_j$ に含まれます（$s_j = t_j$ なら空）。こうして $[0, t_j)$ は、高々 $(b_j-1)L_j + 1$ 個の交わりのない区間（基本区間と、高々 1 個の「端切れ」$[s_j, t_j)$）に分かれます。

**箱の分解。** 各座標の区間の積を取ると、$[0, t)$ は高々 $\prod_j((b_j-1)L_j + 1)$ 個の交わりのない箱 $B$ に分かれます。各 $B$ で $|A(B; P_N) - N\,\mathrm{vol}(B)| \le 1$ を示します。

- $B$ の因子がすべて基本区間のとき。$B$ は基本的な箱なので、はじめの評価そのものです。
- $B$ の因子に端切れがあるとき。端切れを、それを含む階数 $L_j$ の基本区間 $R_j$ で置き換えた箱を $E'$ とすると $B \subseteq E'$ で、$E'$ は基本的な箱です。端切れのある座標 $j_0$ について $Q' \ge b_{j_0}^{L_{j_0}} \ge N$ なので $A(B; P_N) \le A(E'; P_N) \le 1$、また $N\,\mathrm{vol}(B) \le N\,\mathrm{vol}(E') = N/Q' \le 1$ です。$A(B; P_N)$ と $N\,\mathrm{vol}(B)$ はどちらも $[0,1]$ にあるので、差は $1$ 以下です。

点の個数も体積も交わりのない和について加法的なので、
$$
\bigl|A([0,t); P_N) - N\textstyle\prod_jt_j\bigr| \le \sum_B\bigl|A(B; P_N) - N\,\mathrm{vol}(B)\bigr| \le \prod_{j=1}^d\bigl((b_j-1)L_j + 1\bigr)
$$
です。$t$ について上限を取れば主張を得ます。
</Proof>

<Ref to="thm-koksma-hlawka" /> と合わせると、$C^d$ 級の関数の積分誤差は $O(N^{-1}(\log N)^d)$ です。基数は、小さいほうから $d$ 個の素数 $2, 3, 5, \ldots$ を使うのが普通です。Hammersley 点集合は、第 1 座標を $n/N$ にすることで対数の指数が 1 つ下がり、$D_N^* = O(N^{-1}(\log N)^{d-1})$ になります（証明は上と同じ数え方です。$N$ を先に決める必要があります）。

<Remark id="rem-discrepancy-lower-bounds" title="下界と、ほかの構成">
食い違い量はどこまで小さくできるか。$d \ge 2$ の任意の $N$ 点集合について $D_N^* \ge c_d\,N^{-1}(\log N)^{(d-1)/2}$ が成り立ちます（K. F. Roth, *Mathematika* 1 (1954)）。$d = 2$ では $D_N^* \ge c\,N^{-1}\log N$、したがって 1 次元の無限列では無限個の $N$ で $D_N^* \ge c\,N^{-1}\log N$ です（W. M. Schmidt, *Acta Arith.* 21 (1972)）。van der Corput 列は 1 次元の列として次数が最良です。$d \ge 3$ の点集合で $N^{-1}(\log N)^{d-1}$ が最良かどうかは未解決です。これらの下界の証明は直交関数系を使う本格的な議論で、この記事では定理にしません。

Halton 列の定数 $\prod_j(b_j - 1)/\log b_j$ は次元とともに急激に大きくなります。定数が次元について良い構成に、基数 $2$ の Sobol' 列と、その一般化である $(t, m, s)$ ネット・$(t, s)$ 列（Niederreiter の本の第 4 章）があります。これらに桁ごとの無作為な置換を施す Owen のスクランブルは、滑らかな被積分関数に対して分散を $O(N^{-3}(\log N)^{d-1})$ まで下げます（A. B. Owen, *Ann. Statist.* 25 (1997)）。
</Remark>

<Example id="ex-halton-correlation" title="大きな基数の Halton 列の相関">
基数 $17$ と $19$ の Halton 列の 2 次元の射影を考えます。$n < 17$ では $\varphi_{17}(n) = n/17$、$\varphi_{19}(n) = n/19$ なので、最初の 17 点はすべて直線 $x_2 = (17/19)x_1$ の上に並びます。箱 $J = [0, 1/2) \times [1/2, 1)$ に入るには $n/19 \ge 1/2$、すなわち $n \ge 10$ が要りますが、そのとき $n/17 \ge 10/17 > 1/2$ なので、$J$ には 1 点も入りません。体積は $1/4$ なので $D_{17} \ge 1/4$ です。<Ref to="thm-halton" /> の上界は $N = 17$ で $(16 + 1)(18 + 1)/17 = 19$ となり、何も言っていません。上界が意味をもつのは $N$ が $17 \times 19 = 323$ よりずっと大きくなってからです。

高次元の Halton 列では、後ろのほうの座標（大きな素数を基数とする座標）の組にこの種の相関が残ります。対策は、桁に置換を施すこと（スクランブルした Halton 列）や、Sobol' 列のように全座標で基数 $2$ を使う構成に替えることです。「準 Monte Carlo 法は高次元で効かない」と言われるときの原因の 1 つがこれで、もう 1 つの面は §7.2 で扱います。
</Example>

<Figure caption="[0,1)² の 89 点。左から、独立な一様乱数（種 20261008）、基数 2・3 の Halton 列の最初の 89 項、N = 89・生成ベクトル (1, 55) の Fibonacci 格子。乱数には塊と隙間があり、Halton 列と格子では点が互いを避けて並ぶ。">
<svg viewBox="0 0 510 190" width="100%" role="img" aria-label="単位正方形の 89 点を 3 通りに置いた図。乱数は塊と隙間があり、Halton 列と Fibonacci 格子は均等に並ぶ">
  <rect x="10" y="10" width="150" height="150" fill="none" stroke="currentColor" stroke-width="1" />
  <g fill="currentColor"><circle cx="28.8" cy="14.4" r="1.8"/><circle cx="158.7" cy="100.5" r="1.8"/><circle cx="124.3" cy="143.6" r="1.8"/><circle cx="60.5" cy="133.0" r="1.8"/><circle cx="150.6" cy="41.7" r="1.8"/><circle cx="152.5" cy="81.6" r="1.8"/><circle cx="67.1" cy="146.0" r="1.8"/><circle cx="56.5" cy="123.2" r="1.8"/><circle cx="84.6" cy="150.8" r="1.8"/><circle cx="109.9" cy="70.6" r="1.8"/><circle cx="78.3" cy="32.3" r="1.8"/><circle cx="25.5" cy="131.5" r="1.8"/><circle cx="92.0" cy="151.4" r="1.8"/><circle cx="42.8" cy="43.8" r="1.8"/><circle cx="158.9" cy="54.3" r="1.8"/><circle cx="108.6" cy="133.7" r="1.8"/><circle cx="100.5" cy="137.0" r="1.8"/><circle cx="126.0" cy="19.0" r="1.8"/><circle cx="127.5" cy="63.0" r="1.8"/><circle cx="121.1" cy="81.9" r="1.8"/><circle cx="112.2" cy="77.1" r="1.8"/><circle cx="39.5" cy="74.1" r="1.8"/><circle cx="97.9" cy="47.8" r="1.8"/><circle cx="43.9" cy="66.1" r="1.8"/><circle cx="159.4" cy="115.8" r="1.8"/><circle cx="108.2" cy="145.2" r="1.8"/><circle cx="13.2" cy="95.5" r="1.8"/><circle cx="46.7" cy="39.2" r="1.8"/><circle cx="72.9" cy="148.2" r="1.8"/><circle cx="61.7" cy="154.1" r="1.8"/><circle cx="33.4" cy="81.4" r="1.8"/><circle cx="158.9" cy="114.3" r="1.8"/><circle cx="26.4" cy="67.4" r="1.8"/><circle cx="66.3" cy="76.7" r="1.8"/><circle cx="158.1" cy="140.3" r="1.8"/><circle cx="96.2" cy="150.5" r="1.8"/><circle cx="50.0" cy="41.4" r="1.8"/><circle cx="32.4" cy="114.0" r="1.8"/><circle cx="83.4" cy="154.3" r="1.8"/><circle cx="157.2" cy="35.0" r="1.8"/><circle cx="23.5" cy="36.4" r="1.8"/><circle cx="34.7" cy="74.8" r="1.8"/><circle cx="108.1" cy="150.3" r="1.8"/><circle cx="59.5" cy="132.6" r="1.8"/><circle cx="12.8" cy="47.9" r="1.8"/><circle cx="19.3" cy="137.6" r="1.8"/><circle cx="20.2" cy="37.5" r="1.8"/><circle cx="124.2" cy="98.4" r="1.8"/><circle cx="31.0" cy="46.8" r="1.8"/><circle cx="122.6" cy="69.5" r="1.8"/><circle cx="82.6" cy="19.9" r="1.8"/><circle cx="67.3" cy="42.7" r="1.8"/><circle cx="84.5" cy="153.6" r="1.8"/><circle cx="52.5" cy="28.5" r="1.8"/><circle cx="55.8" cy="120.9" r="1.8"/><circle cx="147.0" cy="110.7" r="1.8"/><circle cx="76.0" cy="95.3" r="1.8"/><circle cx="120.5" cy="52.5" r="1.8"/><circle cx="35.3" cy="146.3" r="1.8"/><circle cx="65.0" cy="110.9" r="1.8"/><circle cx="91.4" cy="139.9" r="1.8"/><circle cx="138.9" cy="77.9" r="1.8"/><circle cx="116.1" cy="29.8" r="1.8"/><circle cx="59.6" cy="149.7" r="1.8"/><circle cx="27.7" cy="150.9" r="1.8"/><circle cx="38.6" cy="14.1" r="1.8"/><circle cx="118.9" cy="57.7" r="1.8"/><circle cx="139.1" cy="62.8" r="1.8"/><circle cx="96.4" cy="97.4" r="1.8"/><circle cx="36.0" cy="88.4" r="1.8"/><circle cx="49.5" cy="96.0" r="1.8"/><circle cx="89.5" cy="151.6" r="1.8"/><circle cx="74.2" cy="140.3" r="1.8"/><circle cx="55.9" cy="100.5" r="1.8"/><circle cx="49.3" cy="38.6" r="1.8"/><circle cx="76.0" cy="151.3" r="1.8"/><circle cx="69.6" cy="59.1" r="1.8"/><circle cx="64.5" cy="11.9" r="1.8"/><circle cx="116.3" cy="31.8" r="1.8"/><circle cx="10.5" cy="55.5" r="1.8"/><circle cx="138.4" cy="142.8" r="1.8"/><circle cx="126.8" cy="138.3" r="1.8"/><circle cx="95.8" cy="88.4" r="1.8"/><circle cx="17.7" cy="54.3" r="1.8"/><circle cx="19.9" cy="121.0" r="1.8"/><circle cx="115.1" cy="99.5" r="1.8"/><circle cx="19.2" cy="44.0" r="1.8"/><circle cx="62.8" cy="55.9" r="1.8"/><circle cx="152.3" cy="59.2" r="1.8"/></g>
  <text x="85" y="180" font-size="12" text-anchor="middle" fill="currentColor">独立な一様乱数</text>
  <rect x="180" y="10" width="150" height="150" fill="none" stroke="currentColor" stroke-width="1" />
  <g fill="currentColor"><circle cx="180.0" cy="160.0" r="1.8"/><circle cx="255.0" cy="110.0" r="1.8"/><circle cx="217.5" cy="60.0" r="1.8"/><circle cx="292.5" cy="143.3" r="1.8"/><circle cx="198.8" cy="93.3" r="1.8"/><circle cx="273.8" cy="43.3" r="1.8"/><circle cx="236.2" cy="126.7" r="1.8"/><circle cx="311.2" cy="76.7" r="1.8"/><circle cx="189.4" cy="26.7" r="1.8"/><circle cx="264.4" cy="154.4" r="1.8"/><circle cx="226.9" cy="104.4" r="1.8"/><circle cx="301.9" cy="54.4" r="1.8"/><circle cx="208.1" cy="137.8" r="1.8"/><circle cx="283.1" cy="87.8" r="1.8"/><circle cx="245.6" cy="37.8" r="1.8"/><circle cx="320.6" cy="121.1" r="1.8"/><circle cx="184.7" cy="71.1" r="1.8"/><circle cx="259.7" cy="21.1" r="1.8"/><circle cx="222.2" cy="148.9" r="1.8"/><circle cx="297.2" cy="98.9" r="1.8"/><circle cx="203.4" cy="48.9" r="1.8"/><circle cx="278.4" cy="132.2" r="1.8"/><circle cx="240.9" cy="82.2" r="1.8"/><circle cx="315.9" cy="32.2" r="1.8"/><circle cx="194.1" cy="115.6" r="1.8"/><circle cx="269.1" cy="65.6" r="1.8"/><circle cx="231.6" cy="15.6" r="1.8"/><circle cx="306.6" cy="158.1" r="1.8"/><circle cx="212.8" cy="108.1" r="1.8"/><circle cx="287.8" cy="58.1" r="1.8"/><circle cx="250.3" cy="141.5" r="1.8"/><circle cx="325.3" cy="91.5" r="1.8"/><circle cx="182.3" cy="41.5" r="1.8"/><circle cx="257.3" cy="124.8" r="1.8"/><circle cx="219.8" cy="74.8" r="1.8"/><circle cx="294.8" cy="24.8" r="1.8"/><circle cx="201.1" cy="152.6" r="1.8"/><circle cx="276.1" cy="102.6" r="1.8"/><circle cx="238.6" cy="52.6" r="1.8"/><circle cx="313.6" cy="135.9" r="1.8"/><circle cx="191.7" cy="85.9" r="1.8"/><circle cx="266.7" cy="35.9" r="1.8"/><circle cx="229.2" cy="119.3" r="1.8"/><circle cx="304.2" cy="69.3" r="1.8"/><circle cx="210.5" cy="19.3" r="1.8"/><circle cx="285.5" cy="147.0" r="1.8"/><circle cx="248.0" cy="97.0" r="1.8"/><circle cx="323.0" cy="47.0" r="1.8"/><circle cx="187.0" cy="130.4" r="1.8"/><circle cx="262.0" cy="80.4" r="1.8"/><circle cx="224.5" cy="30.4" r="1.8"/><circle cx="299.5" cy="113.7" r="1.8"/><circle cx="205.8" cy="63.7" r="1.8"/><circle cx="280.8" cy="13.7" r="1.8"/><circle cx="243.3" cy="156.3" r="1.8"/><circle cx="318.3" cy="106.3" r="1.8"/><circle cx="196.4" cy="56.3" r="1.8"/><circle cx="271.4" cy="139.6" r="1.8"/><circle cx="233.9" cy="89.6" r="1.8"/><circle cx="308.9" cy="39.6" r="1.8"/><circle cx="215.2" cy="123.0" r="1.8"/><circle cx="290.2" cy="73.0" r="1.8"/><circle cx="252.7" cy="23.0" r="1.8"/><circle cx="327.7" cy="150.7" r="1.8"/><circle cx="181.2" cy="100.7" r="1.8"/><circle cx="256.2" cy="50.7" r="1.8"/><circle cx="218.7" cy="134.1" r="1.8"/><circle cx="293.7" cy="84.1" r="1.8"/><circle cx="199.9" cy="34.1" r="1.8"/><circle cx="274.9" cy="117.4" r="1.8"/><circle cx="237.4" cy="67.4" r="1.8"/><circle cx="312.4" cy="17.4" r="1.8"/><circle cx="190.5" cy="145.2" r="1.8"/><circle cx="265.5" cy="95.2" r="1.8"/><circle cx="228.0" cy="45.2" r="1.8"/><circle cx="303.0" cy="128.5" r="1.8"/><circle cx="209.3" cy="78.5" r="1.8"/><circle cx="284.3" cy="28.5" r="1.8"/><circle cx="246.8" cy="111.9" r="1.8"/><circle cx="321.8" cy="61.9" r="1.8"/><circle cx="185.9" cy="11.9" r="1.8"/><circle cx="260.9" cy="159.4" r="1.8"/><circle cx="223.4" cy="109.4" r="1.8"/><circle cx="298.4" cy="59.4" r="1.8"/><circle cx="204.6" cy="142.7" r="1.8"/><circle cx="279.6" cy="92.7" r="1.8"/><circle cx="242.1" cy="42.7" r="1.8"/><circle cx="317.1" cy="126.0" r="1.8"/><circle cx="195.2" cy="76.0" r="1.8"/></g>
  <text x="255" y="180" font-size="12" text-anchor="middle" fill="currentColor">Halton 列（基数 2, 3）</text>
  <rect x="350" y="10" width="150" height="150" fill="none" stroke="currentColor" stroke-width="1" />
  <g fill="currentColor"><circle cx="350.0" cy="160.0" r="1.8"/><circle cx="351.7" cy="67.3" r="1.8"/><circle cx="353.4" cy="124.6" r="1.8"/><circle cx="355.1" cy="31.9" r="1.8"/><circle cx="356.7" cy="89.2" r="1.8"/><circle cx="358.4" cy="146.5" r="1.8"/><circle cx="360.1" cy="53.8" r="1.8"/><circle cx="361.8" cy="111.1" r="1.8"/><circle cx="363.5" cy="18.4" r="1.8"/><circle cx="365.2" cy="75.7" r="1.8"/><circle cx="366.9" cy="133.0" r="1.8"/><circle cx="368.5" cy="40.3" r="1.8"/><circle cx="370.2" cy="97.6" r="1.8"/><circle cx="371.9" cy="154.9" r="1.8"/><circle cx="373.6" cy="62.2" r="1.8"/><circle cx="375.3" cy="119.6" r="1.8"/><circle cx="377.0" cy="26.9" r="1.8"/><circle cx="378.7" cy="84.2" r="1.8"/><circle cx="380.3" cy="141.5" r="1.8"/><circle cx="382.0" cy="48.8" r="1.8"/><circle cx="383.7" cy="106.1" r="1.8"/><circle cx="385.4" cy="13.4" r="1.8"/><circle cx="387.1" cy="70.7" r="1.8"/><circle cx="388.8" cy="128.0" r="1.8"/><circle cx="390.4" cy="35.3" r="1.8"/><circle cx="392.1" cy="92.6" r="1.8"/><circle cx="393.8" cy="149.9" r="1.8"/><circle cx="395.5" cy="57.2" r="1.8"/><circle cx="397.2" cy="114.5" r="1.8"/><circle cx="398.9" cy="21.8" r="1.8"/><circle cx="400.6" cy="79.1" r="1.8"/><circle cx="402.2" cy="136.4" r="1.8"/><circle cx="403.9" cy="43.7" r="1.8"/><circle cx="405.6" cy="101.0" r="1.8"/><circle cx="407.3" cy="158.3" r="1.8"/><circle cx="409.0" cy="65.6" r="1.8"/><circle cx="410.7" cy="122.9" r="1.8"/><circle cx="412.4" cy="30.2" r="1.8"/><circle cx="414.0" cy="87.5" r="1.8"/><circle cx="415.7" cy="144.8" r="1.8"/><circle cx="417.4" cy="52.1" r="1.8"/><circle cx="419.1" cy="109.4" r="1.8"/><circle cx="420.8" cy="16.7" r="1.8"/><circle cx="422.5" cy="74.0" r="1.8"/><circle cx="424.2" cy="131.3" r="1.8"/><circle cx="425.8" cy="38.7" r="1.8"/><circle cx="427.5" cy="96.0" r="1.8"/><circle cx="429.2" cy="153.3" r="1.8"/><circle cx="430.9" cy="60.6" r="1.8"/><circle cx="432.6" cy="117.9" r="1.8"/><circle cx="434.3" cy="25.2" r="1.8"/><circle cx="436.0" cy="82.5" r="1.8"/><circle cx="437.6" cy="139.8" r="1.8"/><circle cx="439.3" cy="47.1" r="1.8"/><circle cx="441.0" cy="104.4" r="1.8"/><circle cx="442.7" cy="11.7" r="1.8"/><circle cx="444.4" cy="69.0" r="1.8"/><circle cx="446.1" cy="126.3" r="1.8"/><circle cx="447.8" cy="33.6" r="1.8"/><circle cx="449.4" cy="90.9" r="1.8"/><circle cx="451.1" cy="148.2" r="1.8"/><circle cx="452.8" cy="55.5" r="1.8"/><circle cx="454.5" cy="112.8" r="1.8"/><circle cx="456.2" cy="20.1" r="1.8"/><circle cx="457.9" cy="77.4" r="1.8"/><circle cx="459.6" cy="134.7" r="1.8"/><circle cx="461.2" cy="42.0" r="1.8"/><circle cx="462.9" cy="99.3" r="1.8"/><circle cx="464.6" cy="156.6" r="1.8"/><circle cx="466.3" cy="63.9" r="1.8"/><circle cx="468.0" cy="121.2" r="1.8"/><circle cx="469.7" cy="28.5" r="1.8"/><circle cx="471.3" cy="85.8" r="1.8"/><circle cx="473.0" cy="143.1" r="1.8"/><circle cx="474.7" cy="50.4" r="1.8"/><circle cx="476.4" cy="107.8" r="1.8"/><circle cx="478.1" cy="15.1" r="1.8"/><circle cx="479.8" cy="72.4" r="1.8"/><circle cx="481.5" cy="129.7" r="1.8"/><circle cx="483.1" cy="37.0" r="1.8"/><circle cx="484.8" cy="94.3" r="1.8"/><circle cx="486.5" cy="151.6" r="1.8"/><circle cx="488.2" cy="58.9" r="1.8"/><circle cx="489.9" cy="116.2" r="1.8"/><circle cx="491.6" cy="23.5" r="1.8"/><circle cx="493.3" cy="80.8" r="1.8"/><circle cx="494.9" cy="138.1" r="1.8"/><circle cx="496.6" cy="45.4" r="1.8"/><circle cx="498.3" cy="102.7" r="1.8"/></g>
  <text x="425" y="180" font-size="12" text-anchor="middle" fill="currentColor">Fibonacci 格子 g = (1, 55)</text>
</svg>
</Figure>

## 6. 格子則

Halton 列の誤差 $O(N^{-1}(\log N)^d)$ は、関数が $C^d$ 級より滑らかでも改善しません。周期的で滑らかな関数に対しては、滑らかさに応じて速く収束する点集合があります。1 次元の周期関数に対する台形則が代数的な次数を超えて速く収束すること（<Ref to="mathematics/numerical-analysis/extrapolation-and-romberg#thm-periodic-trapezoid" />）の多次元版です。

<Definition id="def-lattice-rule" title="階数 1 の格子則・双対格子">
$N \ge 1$、$g = (g_1, \ldots, g_d) \in \mathbb{Z}^d$ とする。点集合 $x_n = \{ng/N\} = (\{ng_1/N\}, \ldots, \{ng_d/N\})$（$n = 0, \ldots, N-1$）による求積
$$
Q_N(f) = \frac1N\sum_{n=0}^{N-1}f\Bigl(\Bigl\{\frac{ng}{N}\Bigr\}\Bigr)
$$
を、**生成ベクトル** $g$ の**階数 1 の格子則**という。
$$
L^\perp = L^\perp(g, N) = \{h \in \mathbb{Z}^d : h\cdot g \equiv 0 \pmod N\}
$$
をその**双対格子**という。
</Definition>

$d = 1$、$g = 1$ は等間隔の点での平均（周期関数に対する台形則）です。$d \ge 2$ では、$N$ 点が 1 本の「巻きついた直線」の上に等間隔に並びます。双対格子の直観は次の定理の 1 で分かります。格子則の誤差は、双対格子の点 $h$ に対応する Fourier 係数だけの和で決まるので、双対格子に $0$ に近い点（$r(h)$ の小さい点）が無いほど誤差は小さくなります。生成ベクトル $g$ を選ぶことは、この双対格子の形を選ぶことです。関数は $\mathbb{R}^d$ 上で各変数について周期 $1$ とし、Fourier 係数を $\hat f(h) = \int_{[0,1]^d}f(x)e(-h\cdot x)\,dx$（$h \in \mathbb{Z}^d$）と書きます。$r(h) = \prod_{j=1}^d\max(1, |h_j|)$ とおきます。

<Lemma id="lem-character-sum" title="指標の直交性">
整数 $k$ と $N \ge 1$ について、$\dfrac1N\displaystyle\sum_{n=0}^{N-1}e\Bigl(\frac{nk}{N}\Bigr)$ は $N \mid k$ のとき $1$、そうでないとき $0$ である。
</Lemma>

<Proof of="lem-character-sum">
$N \mid k$ なら各項が $1$ です。そうでなければ $w = e(k/N) \ne 1$、$w^N = e(k) = 1$ で、等比数列の和は $(w^N - 1)/(w - 1) = 0$ です。
</Proof>

<Theorem id="thm-lattice-rule-error" title="格子則の誤差と、良い生成ベクトルの存在">
<Ref to="def-lattice-rule" /> の記号で次が成り立つ。

1. $f$ が各変数について周期 $1$ の連続関数で $\sum_{h\in\mathbb{Z}^d}|\hat f(h)| < \infty$ ならば
   $$
   Q_N(f) - \int_{[0,1]^d}f(x)\,dx = \sum_{h\in L^\perp\setminus\{0\}}\hat f(h).
   $$
2. $\alpha > 1$、$C > 0$ とし、周期 $1$ の連続関数で $|\hat f(h)| \le C\,r(h)^{-\alpha}$（すべての $h$）を満たすものの全体を $E_\alpha(C)$ とする。このとき
   $$
   \sup_{f\in E_\alpha(C)}\Bigl|Q_N(f) - \int_{[0,1]^d}f\Bigr| = C\,P_\alpha(g, N), \qquad P_\alpha(g, N) = \sum_{h\in L^\perp\setminus\{0\}}r(h)^{-\alpha}.
   $$
3. $N$ を素数で $N \ge 3$ かつ $1 + 1/\ln N \le \alpha$ を満たすものとする。このとき $g \in \{0, 1, \ldots, N-1\}^d$ で
   $$
   P_\alpha(g, N) \le (2e)^\alpha\,\frac{(3 + 2\ln N)^{\alpha d}}{N^\alpha}
   $$
   を満たすものが存在する。
</Theorem>

1 と 2 は Fourier 級数と指標の直交性（<Ref to="lem-character-sum" />）から直接出ます。3 は「良い $g$ を具体的に作る」のではなく、「すべての $g$ について $P_\beta$ を平均すると十分小さい。最小値は平均以下なので、平均以下の $g$ が少なくとも 1 つある」という論法で存在を示します。まず $\beta$ で平均を評価し、次に不等式で $\beta$ から $\alpha$ に移ります。

<Proof of="thm-lattice-rule-error">
**1.** $S(x) = \sum_h\hat f(h)e(h\cdot x)$ とおきます。$|\hat f(h)e(h\cdot x)| = |\hat f(h)|$ の和が有限なので、この級数は一様に絶対収束し、$S$ は周期 $1$ の連続関数です。一様収束する級数は項別に積分できるので、$S$ の Fourier 係数は $\hat f(h)$ です。$f - S$ は Fourier 係数がすべて $0$ の連続関数なので、$f = S$ です（<Ref to="mathematics/fourier-analysis/fourier-series#prop-torus-uniqueness" />。変数を $\theta = 2\pi x$ と取り替えて使います）。$x_n = \{ng/N\}$ を代入すると $e(h\cdot x_n) = e(n(h\cdot g)/N)$ で、絶対収束から和の順序を入れ替えて
$$
Q_N(f) = \sum_{h\in\mathbb{Z}^d}\hat f(h)\cdot\frac1N\sum_{n=0}^{N-1}e\Bigl(\frac{n(h\cdot g)}{N}\Bigr) = \sum_{h\in L^\perp}\hat f(h)
$$
です（<Ref to="lem-character-sum" /> を $k = h\cdot g$ で使いました）。$h = 0$ の項は $\hat f(0) = \int f$ です。

**2.** $\sum_{h\in\mathbb{Z}^d}r(h)^{-\alpha} = \prod_{j=1}^d\sum_{h_j\in\mathbb{Z}}\max(1, |h_j|)^{-\alpha} = (1 + 2\zeta(\alpha))^d < \infty$ です（$\zeta(\alpha) = \sum_{k\ge1}k^{-\alpha}$）。したがって $f \in E_\alpha(C)$ は 1 の仮定を満たし、$|Q_N(f) - \int f| \le \sum_{h\in L^\perp\setminus\{0\}}|\hat f(h)| \le C\,P_\alpha(g, N)$ です。等号は $f^*(x) = C\sum_hr(h)^{-\alpha}e(h\cdot x)$ で成り立ちます。この級数は一様収束するので $f^*$ は周期 $1$ の連続関数で、Fourier 係数は $C\,r(h)^{-\alpha}$ です。$f^* \in E_\alpha(C)$ で、1 から誤差はちょうど $C\,P_\alpha(g, N)$ です。

**3.** まず $1 < \beta \le \alpha$ である $\beta$ について、$g$ を $\{0, \ldots, N-1\}^d$ の全体で動かした $P_\beta(g, N)$ の平均 $M_\beta$ を評価します。
$$
M_\beta = \frac{1}{N^d}\sum_g\sum_{h\ne0}\mathbf{1}[h\cdot g \equiv 0]\,r(h)^{-\beta} = \sum_{h\ne0}r(h)^{-\beta}\cdot\frac{\#\{g : h\cdot g \equiv 0 \pmod N\}}{N^d}
$$
です。$h$ の成分のうちに $N$ で割り切れないもの $h_{j}$ があるとき、$N$ が素数なので $h_j$ は $N$ を法として逆元をもち、他の成分 $g_i$（$i \ne j$）を任意に決めると $g_j$ が $N$ を法としてただ 1 つ決まります。解は $N^{d-1}$ 個で、割合は $1/N$ です。すべての成分が $N$ で割り切れるときは、すべての $g$ が解で、割合は $1$ です。後者の $h$ は $h = Nk$（$k \in \mathbb{Z}^d\setminus\{0\}$）と書け、$k_j \ne 0$ の成分で $\max(1, N|k_j|) = N|k_j|$、その他の成分で $1$ なので $r(Nk) \ge N\,r(k)$ です。よって
$$
M_\beta \le \frac1N\sum_{h\ne0}r(h)^{-\beta} + \sum_{k\ne0}\bigl(N\,r(k)\bigr)^{-\beta} \le \frac{2}{N}\bigl(1 + 2\zeta(\beta)\bigr)^d
$$
です（$N^{-\beta} \le N^{-1}$）。最小値は平均以下なので、$P_\beta(g, N) \le M_\beta$ を満たす $g$ があります。

次に $\beta$ から $\alpha$ へ移ります。非負の数 $a_i$ と $p \ge 1$ について $\sum_ia_i^p \le (\sum_ia_i)^p$ です（$S = \sum_ia_i > 0$ として $a_i/S \le 1$ から $(a_i/S)^p \le a_i/S$ を足します）。$a_h = r(h)^{-\beta}$、$p = \alpha/\beta$ とすると、同じ $g$ について
$$
P_\alpha(g, N) = \sum_{h\in L^\perp\setminus\{0\}}\bigl(r(h)^{-\beta}\bigr)^{\alpha/\beta} \le P_\beta(g, N)^{\alpha/\beta} \le M_\beta^{\alpha/\beta}
$$
です。ここで $\beta = 1 + 1/\ln N$ と選びます（仮定から $\beta \le \alpha$）。$\zeta(\beta) \le 1 + \int_1^\infty x^{-\beta}dx = 1 + 1/(\beta - 1) = 1 + \ln N$ なので
$$
M_\beta \le y, \qquad y = \frac{2(3 + 2\ln N)^d}{N} \ge \frac1N
$$
です。$y^{\alpha/\beta} \le e^\alpha y^\alpha$ を示します。$y \ge 1$ なら $\alpha/\beta \le \alpha$ から $y^{\alpha/\beta} \le y^\alpha$ です。$y < 1$ なら $y^{\alpha/\beta} = y^\alpha\cdot y^{-\alpha(1 - 1/\beta)}$ で、$1 - 1/\beta = (\beta-1)/\beta \le \beta - 1 = 1/\ln N$ と $y \ge 1/N$ から $y^{-\alpha(1-1/\beta)} \le N^{\alpha/\ln N} = e^\alpha$ です。以上から
$$
P_\alpha(g, N) \le e^\alpha y^\alpha = (2e)^\alpha\frac{(3 + 2\ln N)^{\alpha d}}{N^\alpha}.
$$
</Proof>

3 の主張は「良い $g$ がある」という存在定理で、平均が小さいことから導いています。ほとんどの $g$ が悪くないということでもあります。実際に $g$ を見つけるには、$P_\alpha$ を計算して探します。$\alpha$ が偶数のとき $P_\alpha$ は Bernoulli 多項式で閉じた形に書け（次の例）、1 つの $g$ について $O(Nd)$ の手間で計算できます。成分を 1 つずつ決めていく構成法（component-by-component 構成）は、Dick–Kuo–Sloan の総説にあります。

収束の次数は $\alpha$、すなわち関数の滑らかさで決まります。$f$ の混合偏導関数が各変数について $\alpha$ 階まで連続かつ周期的なら、部分積分から $|\hat f(h)| \le C\,r(h)^{-\alpha}$ です。滑らかな周期関数ほど速く収束し、Halton 列のような $N^{-1}$ の壁がありません。周期的でない関数に格子則をそのまま使うと、境界での跳びが Fourier 係数の減衰を $r(h)^{-1}$ に落とし、利点が消えます。

次の例は、2 次元で具体的に良い $g$ の一つを挙げるものです。$g = (1, F_{m-1})$ の比 $F_{m-1}/F_m$ は黄金比に近く、§3 で見たとおり黄金比は有理数で最も近似しにくい無理数です。そのため双対格子の点で $r(h)$ が小さいものが現れにくくなる、というのがこの選び方の動機で、例ではそれを実際に確かめます。

<Example id="ex-fibonacci-lattice" title="Fibonacci 格子（d = 2）">
Fibonacci 数を $F_1 = F_2 = 1$、$F_{k+1} = F_k + F_{k-1}$ とし、$N = F_m$、$g = (1, F_{m-1})$ とします。双対格子は $\{h : h_1 + F_{m-1}h_2 \equiv 0 \pmod{F_m}\}$ です。恒等式 $F_{m-1}F_k - F_mF_{k-1} = (-1)^{k-1}F_{m-k}$（$1 \le k \le m-1$。$k$ についての帰納法で確かめられます）から、
$$
h^{(k)} = \bigl((-1)^kF_{m-k},\ F_k\bigr) \in L^\perp, \qquad r(h^{(k)}) = F_kF_{m-k}
$$
です。$F_kF_{m-k} + F_{k-1}F_{m-k-1} = F_{m-1}$ から $r(h^{(k)}) \ge F_{m-1}/2$ で、これらの双対格子の点はどれも $r(h)$ が $N$ に比例する大きさです。双対格子の $0$ でない点すべてについて $r(h) \ge cN$（$c > 0$ は $m$ によらない定数）が成り立つことが連分数の理論から分かり、それを使うと $P_\alpha = O(N^{-\alpha}\log N)$ が出ます（Niederreiter の本の第 5 章、Sloan–Joe の本）。2 次元では対数が 1 個で済みます。

$\alpha = 2$ では、$\sum_{h\in\mathbb{Z}}\max(1,|h|)^{-2}e(hx) = 1 + 2\pi^2B_2(\{x\})$（$B_2(x) = x^2 - x + 1/6$ は Bernoulli 多項式）なので、<Ref to="thm-lattice-rule-error" /> の 2 の $f^*$（$C = 1$）は $\prod_j(1 + 2\pi^2B_2(\{x_j\}))$ で、
$$
P_2(g, N) = \frac1N\sum_{n=0}^{N-1}\prod_{j=1}^d\Bigl(1 + 2\pi^2B_2\Bigl(\Bigl\{\frac{ng_j}{N}\Bigr\}\Bigr)\Bigr) - 1
$$
と有限和で計算できます。§8 で Fibonacci 格子についてこの値を計算し、$N^2P_2/\ln N$ がほぼ一定であることを確かめます。
</Example>

## 7. 無作為化と有効次元

### 7.1. 無作為な平行移動

準 Monte Carlo 法の点集合は決定的なので、1 回の計算では誤差の大きさが分かりません。<Ref to="thm-koksma-hlawka" /> の上界は計算できないうえに過大です。点集合全体を無作為にずらすと、統計的な誤差評価が戻ってきます。

<Definition id="def-random-shift" title="無作為な平行移動（Cranley–Patterson）">
点集合 $P = \{x_1, \ldots, x_N\} \subset [0,1)^d$ と、$[0,1)^d$ 上の一様分布に従う確率変数 $U$ について、点集合 $\{x_n + U\}$（成分ごとに小数部分を取る）を $P$ の**無作為な平行移動**といい、
$$
Q_U(f) = \frac1N\sum_{n=1}^Nf(\{x_n + U\})
$$
をそれによる推定量という。
</Definition>

<Proposition id="prop-random-shift" title="無作為な平行移動による推定量は不偏">
$f : [0,1)^d \to \mathbb{R}$ を可積分な関数とし、$I = \int_{[0,1)^d}f$ とする。

1. 各 $x \in [0,1)^d$ について $\{x + U\}$ は $[0,1)^d$ 上の一様分布に従う。したがって $\mathbb{E}[Q_U(f)] = I$。
2. $f$ が 2 乗可積分とする。独立な平行移動 $U_1, \ldots, U_M$（$M \ge 2$）について、$\bar Q = M^{-1}\sum_kQ_{U_k}(f)$ は $\mathbb{E}[\bar Q] = I$、$\mathrm{Var}(\bar Q) = \mathrm{Var}(Q_U(f))/M$ を満たし、標本分散 $s^2 = (M-1)^{-1}\sum_k(Q_{U_k}(f) - \bar Q)^2$ について $\mathbb{E}[s^2/M] = \mathrm{Var}(\bar Q)$ である。
</Proposition>

<Proof of="prop-random-shift">
**1.** $[0,1)^d$ 上の有界な可測関数 $\phi$ を、各変数について周期 $1$ に延長した関数も $\phi$ と書きます。$\phi(\{x + u\}) = \phi(x + u)$ です。1 変数では、周期 $1$ の関数の長さ $1$ の区間での積分は区間の位置によらないので $\int_0^1\phi(x + u)\,du = \int_x^{x+1}\phi(v)\,dv = \int_0^1\phi(v)\,dv$ です。$d$ 変数では、Fubini の定理（<Ref to="mathematics/real-analysis/product-measures-fubini#thm-fubini" />）で座標ごとにこの置き換えを行って
$$
\mathbb{E}[\phi(\{x + U\})] = \int_{[0,1)^d}\phi(x + u)\,du = \int_{[0,1)^d}\phi(v)\,dv
$$
です。$\phi$ を可測集合の定義関数にすれば、$\{x + U\}$ の分布が一様分布であることになります。したがって可積分な $f$ について $\mathbb{E}[f(\{x_n + U\})] = I$ で、$n$ について平均して $\mathbb{E}[Q_U(f)] = I$ です。

**2.** $Y_k = Q_{U_k}(f)$ は独立同分布です。$(N^{-1}\sum_na_n)^2 \le N^{-1}\sum_na_n^2$ と 1 から $\mathbb{E}[Y_k^2] \le N^{-1}\sum_n\mathbb{E}[f(\{x_n + U\})^2] = \int f^2 < \infty$ なので分散は有限です。独立同分布な標本の平均と標本分散についての主張（<Ref to="computer-science/stochastic-computation/monte-carlo-integration#thm-mc-error" /> の 1〜3）から結論を得ます。
</Proof>

平行移動しても、格子則の点集合は格子のままです（双対格子は変わりません）。一般の点集合でも、箱を平行移動すると高々 $2^d$ 個の箱に分かれるので、食い違い量は定数倍しか悪くなりません。つまり、点集合の良さを保ったまま無作為性が入ります。使い方は次のとおりです。全部で $NM$ 回の関数評価を、$M$ 個（10〜50 個）の独立な平行移動に分け、$M$ 個の値の平均を答え、$s/\sqrt M$ を標準誤差とします。$M$ を大きくしすぎると $N$ が小さくなって準 Monte Carlo 法の利点が減るので、$M$ は誤差の推定に必要な最小限にします。

### 7.2. ANOVA 分解と有効次元

<Ref to="thm-halton" /> の上界 $(\log N)^d/N$ が $N^{-1/2}$ より小さくなるのは、$N$ が $d$ について指数的に大きくなってからです。それでも、金融の経路依存型の積分のような数百次元の問題で、準 Monte Carlo 法は実際に Monte Carlo 法より良い結果を出します。その理由は、被積分関数が「見かけの次元」ほど高次元でないことにあります。これを分散の分解で表します。

考え方は次のとおりです。多変数の関数 $f$ を、「$x_1$ だけに依存する部分」「$x_2$ だけに依存する部分」「$x_1, x_2$ の組に依存する部分」…というように、依存する変数の組ごとの部分に分けます。各部分の分散を足すと全体の分散になるようにできます。そのとき、依存する変数の数が少ない部分に分散のほとんどがあるなら、見かけの次元が高くても実質的には低次元の積分とみなせます。次の定義と命題はこれを正確に述べたものです。

$D = \{1, \ldots, d\}$ とし、$f \in L^2([0,1]^d)$、$u \subseteq D$ について、$u$ に属さない変数で積分した関数を
$$
(\Pi_uf)(x_u) = \int_{[0,1]^{D\setminus u}}f(x)\,dx_{D\setminus u}
$$
と書きます（$\Pi_\emptyset f = I = \int f$、$\Pi_Df = f$）。

<Definition id="def-anova" title="ANOVA 分解・有効次元">
$f \in L^2([0,1]^d)$ について、$|u|$ の小さい順に
$$
f_u = \Pi_uf - \sum_{v\subsetneq u}f_v
$$
と定め、$f = \sum_{u\subseteq D}f_u$ を $f$ の **ANOVA 分解**という。$\sigma_u^2 = \int_{[0,1]^d}f_u^2$（$u \ne \emptyset$）、$\sigma^2 = \int(f - I)^2$ とする。$0 < p < 1$ について、$\sum_{0<|u|\le s}\sigma_u^2 \ge p\,\sigma^2$ となる最小の $s$ を、$f$ の（重ね合わせの意味の）**有効次元**という（$p = 0.99$ がよく使われる）。
</Definition>

<Proposition id="prop-effective-dimension" title="ANOVA 分解の直交性と、無作為化した準 Monte Carlo 法の誤差">
$f \in L^2([0,1]^d)$ とし、<Ref to="def-anova" /> の記号を使う。

1. $f_u$ は $x_u$ だけの関数で、$j \in u$ ならば $\int_0^1f_u\,dx_j = 0$。
2. $u \ne v$ ならば $\int_{[0,1]^d}f_uf_v\,dx = 0$。したがって $\sigma^2 = \sum_{\emptyset\ne u\subseteq D}\sigma_u^2$。
3. $1 \le s \le d$ とし、$f_{\le s} = \sum_{|u|\le s}f_u$ とおく。点集合 $P \subset [0,1)^d$ の無作為な平行移動による推定量 $Q_U$ について
   $$
   \sqrt{\mathbb{E}\bigl[(Q_U(f) - I)^2\bigr]} \le \sqrt{\mathbb{E}\bigl[(Q_U(f_{\le s}) - I)^2\bigr]} + \Bigl(\sum_{|u|>s}\sigma_u^2\Bigr)^{1/2}.
   $$
</Proposition>

<Proof of="prop-effective-dimension">
**1.** $\Pi_uf$ は $x_u$ だけの関数で、$|u|$ についての帰納法から $f_u$ もそうです（$v \subsetneq u$ の $f_v$ は $x_v$ だけの関数）。2 乗可積分性も同様です（$\Pi_u$ は平均を取る操作で、Cauchy–Schwarz の不等式から $\int(\Pi_uf)^2 \le \int f^2$）。積分が $0$ になることを $|u|$ についての帰納法で示します。$j \in u$ とします。Fubini の定理（<Ref to="mathematics/real-analysis/product-measures-fubini#thm-fubini" />）から $\int_0^1\Pi_uf\,dx_j = \Pi_{u\setminus\{j\}}f$ です。$v \subsetneq u$ について、$j \notin v$ なら $f_v$ は $x_j$ によらないので $\int_0^1f_v\,dx_j = f_v$、$j \in v$ なら帰納法の仮定（$|v| < |u|$）から $\int_0^1f_v\,dx_j = 0$ です。$j \notin v$ かつ $v \subsetneq u$ となる $v$ は、$u\setminus\{j\}$ の部分集合の全体なので
$$
\int_0^1f_u\,dx_j = \Pi_{u\setminus\{j\}}f - \sum_{v\subseteq u\setminus\{j\}}f_v = \Pi_{u\setminus\{j\}}f - f_{u\setminus\{j\}} - \sum_{v\subsetneq u\setminus\{j\}}f_v = 0
$$
です。最後の等号は $f_{u\setminus\{j\}}$ の定義です。

**2.** $u \ne v$ なら、一方だけに属する添字があります。$j \in u\setminus v$ とします（逆の場合も同じです）。$f_uf_v$ は可積分で、$f_v$ は $x_j$ によらないので、Fubini の定理で $x_j$ について先に積分すると、1 から
$$
\int f_uf_v\,dx = \int\Bigl(\int_0^1f_u\,dx_j\Bigr)f_v\,dx_{D\setminus\{j\}} = 0
$$
です。$f_D$ の定義から $f = \sum_{u\subseteq D}f_u$、$f_\emptyset = I$ なので $f - I = \sum_{u\ne\emptyset}f_u$ で、直交性から $\sigma^2 = \sum_{u\ne\emptyset}\sigma_u^2$ です。

**3.** $g = f - f_{\le s} = \sum_{|u|>s}f_u$ とおくと、2 から $\int g^2 = \sum_{|u|>s}\sigma_u^2$ です。$Q_U$ は線形なので $Q_U(f) - I = (Q_U(f_{\le s}) - I) + Q_U(g)$ で、2 乗平均のノルムについての三角不等式（Minkowski の不等式、<Ref to="mathematics/real-analysis/lp-spaces#thm-minkowski" />）から
$$
\sqrt{\mathbb{E}[(Q_U(f) - I)^2]} \le \sqrt{\mathbb{E}[(Q_U(f_{\le s}) - I)^2]} + \sqrt{\mathbb{E}[Q_U(g)^2]}
$$
です。$(N^{-1}\sum_na_n)^2 \le N^{-1}\sum_na_n^2$ と、$\{x_n + U\}$ が一様分布に従うこと（<Ref to="prop-random-shift" /> の 1）から
$$
\mathbb{E}[Q_U(g)^2] \le \frac1N\sum_{n=1}^N\mathbb{E}\bigl[g(\{x_n + U\})^2\bigr] = \int g^2 = \sum_{|u|>s}\sigma_u^2
$$
です。
</Proof>

3 の右辺の第 1 項は、高々 $s$ 個の変数の関数の和 $f_{\le s}$ に対する誤差です。$f_u$ の積分には点集合の $u$ 座標への射影しか効かず、射影の食い違い量は $|u|$ 次元のものなので、$f_u$ が滑らかなら誤差は <Ref to="thm-koksma-hlawka" /> から $O(N^{-1}(\log N)^{|u|})$ です。第 2 項は高次の成分の大きさで、$N$ によらず、しかも Monte Carlo 法の誤差 $\sigma/\sqrt N$ と違って $N$ で減りませんが、有効次元が小さければはじめから小さい量です。つまり、準 Monte Carlo 法が高次元で効く条件は「$d$ が小さいこと」ではなく、「分散の大部分が少数の変数の成分にあり、点集合の低次元の射影が良いこと」です。変数の順序や取り方を変えて有効次元を下げる工夫（Brown 運動を Brown 橋で構成して、重要な変動をはじめの座標に集めるなど）が効くのもこのためです。

<Example id="ex-product-anova" title="積の形の関数の ANOVA 分解">
$f(x) = \prod_{j=1}^d(1 + a_j(x_j - 1/2))$ とします。$\int_0^1(x_j - 1/2)\,dx_j = 0$ なので、積を展開した各項 $\prod_{j\in u}a_j(x_j - 1/2)$ は <Ref to="prop-effective-dimension" /> の 1 の性質をもち、これが $f_u$ です（分解の一意性は <Ref to="exr-anova-xy" /> と同じ計算で確かめられます）。$I = 1$ で、
$$
\sigma_u^2 = \prod_{j\in u}\frac{a_j^2}{12}, \qquad \sigma^2 = \prod_{j=1}^d\Bigl(1 + \frac{a_j^2}{12}\Bigr) - 1
$$
です。$d = 5$、$a_j = 1/j$ では $\sigma^2 = 0.12570$、1 変数の成分の和は $\sum_ja_j^2/12 = 0.12197$ で全体の $97.0\%$、2 変数の成分まで入れると $99.97\%$ です。見かけは 5 次元でも、有効次元は $2$（$p = 0.99$）です。$a_j = 1$ で $d = 100$ とすると、$\sigma^2 = (13/12)^{100} - 1 \approx 2.99\times10^3$ に対し 1 変数の成分の和は $100/12 = 8.3$ で、分散のほとんどが高次の成分にあります。同じ形の関数でも、係数の減り方で難しさがまったく違います。
</Example>

## 8. 数値実験

コードは Python 3.9.6・NumPy 2.0.2・SciPy 1.13.1 で実行し、出力をそのまま載せています。乱数は `numpy.random.default_rng(20261008)` です。

### 8.1. 食い違い量の厳密な値

1 次元は <Ref to="prop-discrepancy-1d" /> の式で計算します。2 次元では、点の座標の作る格子の各長方形の内部で $A([0,t); P)$ が一定で、$\delta_P(t) = \text{定数} - t_1t_2$ は各変数について単調なので、上限は長方形の頂点（格子点）での極限値、すなわち格子点での値（境界の点を含める数え方と含めない数え方の両方）の最大です。これを全格子点について計算します。

```python
import numpy as np

def radical_inverse(n, b):
    x, f = 0.0, 1.0 / b
    while n > 0:
        n, a = divmod(n, b)
        x += a * f
        f /= b
    return x

def star_disc_1d(x):
    x = np.sort(x); n = len(x)
    return 0.5 / n + np.abs(x - (2 * np.arange(1, n + 1) - 1) / (2 * n)).max()

def star_disc_2d(p):
    # 格子 (X ∪ {1}) × (Y ∪ {1}) の全点で、閉じた箱と開いた箱の点の個数を数える
    n = len(p)
    xs, ys = np.append(np.sort(p[:, 0]), 1.0), np.append(np.sort(p[:, 1]), 1.0)
    ix = np.searchsorted(xs[:-1], p[:, 0]); iy = np.searchsorted(ys[:-1], p[:, 1])
    cnt = np.zeros((n + 1, n + 1)); np.add.at(cnt, (ix, iy), 1.0)
    closed = cnt.cumsum(axis=0).cumsum(axis=1)             # #{x <= xs[i], y <= ys[j]}
    opened = np.zeros_like(closed); opened[1:, 1:] = closed[:-1, :-1]   # #{x < xs[i], y < ys[j]}
    vol = np.outer(xs, ys)
    return max((closed / n - vol).max(), (vol - opened / n).max())

print("     N   N*D(vdC,b=2)  bound   N*D(vdC,b=3)  bound   N*D(Halton 2,3)  bound")
for N in (10, 100, 1000, 4000):
    v2 = np.array([radical_inverse(n, 2) for n in range(N)])
    v3 = np.array([radical_inverse(n, 3) for n in range(N)])
    ilog = lambda N, b: len(np.base_repr(N, b)) - 1          # floor(log_b N)
    clog = lambda N, b: len(np.base_repr(N - 1, b))          # ceil(log_b N)  (N >= 2)
    b2, b3 = (ilog(N, 2) + 1), 2 * (ilog(N, 3) + 1)
    bh = (clog(N, 2) + 1) * (2 * clog(N, 3) + 1)
    print(f"{N:6d}   {N * star_disc_1d(v2):8.3f}    {b2:4d}    {N * star_disc_1d(v3):8.3f}    {b3:4d}    "
          f"{N * star_disc_2d(np.column_stack([v2, v3])):8.3f}      {bh:5d}")
```

```text
     N   N*D(vdC,b=2)  bound   N*D(vdC,b=3)  bound   N*D(Halton 2,3)  bound
    10      1.750       4       1.889       6       2.667         35
   100      2.313       7       2.951      10       5.024         88
  1000      2.453      10       3.476      14       6.848        165
  4000      2.453      12       4.749      16       9.136        221
```

`bound` は <Ref to="thm-van-der-corput" /> の $(b-1)(\lfloor\log_bN\rfloor + 1)$ と <Ref to="thm-halton" /> の $\prod_j((b_j-1)\lceil\log_{b_j}N\rceil + 1)$ です。$N D_N^*$ は上界の中にあり、$N$ とともにゆっくり（対数的に）増えています。独立な乱数なら $ND_N^*$ は $\sqrt N$ の大きさ（$N = 4000$ で $\sqrt N \approx 63$）です。2 次元の上界は実際の値の $13$〜$24$ 倍で、上界の定数が粗いことも分かります。

### 8.2. 5 次元の積分

$d = 5$ で、2 つの被積分関数（どちらも積分は $1$）を使います。

- $f_A(x) = \prod_j(1 + (x_j - 1/2)/j)$。滑らかですが周期的ではありません（<Ref to="ex-product-anova" /> の関数）。
- $f_B(x) = \prod_j(1 + 2\pi^2\gamma_jB_2(x_j))$、$\gamma_j = 1/j^2$。周期的で、Fourier 係数は $\prod_{h_j\ne0}\gamma_j|h_j|^{-2}$ です（<Ref to="ex-fibonacci-lattice" /> の関数に重みを付けたもの）。

Monte Carlo 法、Halton 列（基数 $2, 3, 5, 7, 11$）、Sobol' 列（`scipy.stats.qmc`）、格子則を $N = 2^6, 2^8, \ldots, 2^{16}$ で比べます。格子則の生成ベクトルは Korobov 型 $g = (1, a, a^2, a^3, a^4) \bmod N$ とし、$a$ は無作為に選んだ高々 200 個の奇数の候補のうち、$f_B$ に対する誤差（<Ref to="thm-lattice-rule-error" /> の 2 から、重みつきの $P_2$ に等しい）が最小のものにします。

```python
import numpy as np
from scipy.stats import qmc

d = 5
a = 1.0 / np.arange(1, d + 1)
f_a = lambda x: np.prod(1.0 + a * (x - 0.5), axis=1)                    # 非周期。積分 1
b2 = lambda x: x * x - x + 1.0 / 6.0
gam = 1.0 / np.arange(1, d + 1) ** 2
f_b = lambda x: np.prod(1.0 + gam * 2 * np.pi**2 * b2(x), axis=1)     # 周期的。積分 1

def lattice(N, g):
    return (np.arange(N)[:, None] * g[None, :] % N) / N

def korobov(N, rng, tries=200):
    # g = (1, a, a^2, ...) mod N のうち P_2 が最小のもの（a は奇数の候補から探す）
    best = None
    for a_ in np.unique(rng.integers(1, N // 2, size=tries) * 2 + 1):
        g = np.array([pow(int(a_), j, N) for j in range(d)])
        p2 = f_b(lattice(N, g)).mean() - 1.0
        if best is None or p2 < best[0]:
            best = (p2, g)
    return best[1]

rng = np.random.default_rng(20261008)
halton = qmc.Halton(d=d, scramble=False).random(2**16)
sobol = qmc.Sobol(d=d, scramble=False).random_base2(16)
for f, label in ((f_a, "f_A"), (f_b, "f_B")):
    print(label)
    print("   log2N   MC err   (MC se)     Halton      Sobol     lattice    gen. a")
    errs = []
    for k in range(6, 17, 2):
        N = 2**k
        y = f(rng.random((N, d)))
        g = korobov(N, rng)
        e = [abs(y.mean() - 1), abs(f(halton[:N]).mean() - 1), abs(f(sobol[:N]).mean() - 1), abs(f(lattice(N, g)).mean() - 1)]
        errs.append(e)
        print(f"   {k:3d}   {e[0]:.2e}  ({y.std(ddof=1) / np.sqrt(N):.2e})   {e[1]:.2e}   {e[2]:.2e}   {e[3]:.2e}   {g[1]:6d}")
    ks = np.arange(6, 17, 2)
    sl = [np.polyfit(ks, np.log2(np.array(errs)[:, j]), 1)[0] for j in range(4)]
    print("   slope    %.2f                  %.2f       %.2f       %.2f" % tuple(sl))
```

```text
f_A
   log2N   MC err   (MC se)     Halton      Sobol     lattice    gen. a
     6   4.46e-02  (3.93e-02)   4.12e-02   1.79e-02   4.67e-03       23
     8   6.70e-03  (2.23e-02)   1.24e-02   4.60e-03   3.68e-03       69
    10   7.76e-03  (1.15e-02)   3.22e-03   1.10e-03   1.87e-03      469
    12   5.85e-03  (5.58e-03)   1.01e-03   2.81e-04   7.42e-04     3739
    14   1.83e-03  (2.77e-03)   2.88e-04   6.81e-05   9.76e-05     7185
    16   1.20e-03  (1.39e-03)   7.51e-05   1.74e-05   2.08e-05    10759
   slope    -0.46                  -0.91       -1.00       -0.80
f_B
   log2N   MC err   (MC se)     Halton      Sobol     lattice    gen. a
     6   4.68e-01  (2.42e-01)   8.31e-02   7.22e-02   7.42e-02       23
     8   2.20e-01  (9.33e-02)   2.27e-02   1.30e-02   1.02e-02       69
    10   1.17e-02  (5.12e-02)   2.81e-03   2.45e-03   1.27e-03      555
    12   2.37e-02  (2.57e-02)   7.60e-06   5.34e-04   1.68e-04     3543
    14   1.16e-02  (1.31e-02)   1.55e-05   1.39e-04   1.82e-05     8603
    16   2.29e-03  (6.46e-03)   1.52e-05   3.73e-05   2.29e-06     7531
   slope    -0.72                  -1.46       -1.09       -1.50
```

`slope` は $\log_2(\text{誤差})$ を $\log_2N$ に最小二乗であてはめた傾きです。6 点からの値で、決定的な方法の誤差は $N$ について滑らかに減るとは限らないので、目安として読んでください。

- **Monte Carlo 法**の誤差は標準誤差 `(MC se)` の $0.2$〜$2.4$ 倍の範囲にあり、標準誤差は $N$ が 4 倍になるごとに半分になっています（傾き $-1/2$）。
- **Halton 列と Sobol' 列**は $f_A$ で傾きがほぼ $-1$ で、<Ref to="thm-koksma-hlawka" /> と <Ref to="thm-halton" /> の $O(N^{-1}(\log N)^d)$ と整合します。$N = 2^{16}$ での誤差は Monte Carlo 法の標準誤差の約 $1/20$（Halton 列）と $1/80$（Sobol' 列）です。$f_B$ での Halton 列の誤差は $N = 2^{12}$ で偶然小さくなっていて、傾きの値はこの 1 点に引きずられています。
- **格子則**は、周期的な $f_B$ では $N = 2^{16}$ で誤差 $2.3\times10^{-6}$ と最も良く、傾きは $-1.5$ です。<Ref to="thm-lattice-rule-error" /> の 3 の上界 $N^{-2}(\log N)^{10}$ の、この範囲での見かけの傾きは $-2 + 10/\ln N$（$N = 2^{16}$ で $-1.1$）なので、$-2$ より緩いのは対数の因子の分です。周期的でない $f_A$ では格子則に特別な利点はなく、Sobol' 列と同程度かそれ以下です。

### 8.3. Fibonacci 格子と無作為な平行移動

```python
import numpy as np

b2 = lambda x: x * x - x + 1.0 / 6.0
print("Fibonacci lattice (d = 2): P_2 = (1/N) sum prod (1 + 2 pi^2 B_2) - 1")
print("     N      P_2        N^2 P_2 / ln N")
F = [1, 1]
while F[-1] < 200000:
    F.append(F[-1] + F[-2])
for m in range(10, len(F), 3):
    N, g = F[m], np.array([1, F[m - 1]])
    x = (np.arange(N)[:, None] * g[None, :] % N) / N
    p2 = np.prod(1.0 + 2 * np.pi**2 * b2(x), axis=1).mean() - 1.0
    print(f"{N:7d}   {p2:.3e}   {N * N * p2 / np.log(N):8.2f}")

# 無作為な平行移動: d = 5, N = 4096, 生成ベクトル (1, a, a^2, a^3, a^4) mod N, a = 3543
d, N, a_gen, M = 5, 4096, 3543, 32
g = np.array([pow(a_gen, j, N) for j in range(d)])
pts = (np.arange(N)[:, None] * g[None, :] % N) / N
a = 1.0 / np.arange(1, d + 1)
gam = 1.0 / np.arange(1, d + 1) ** 2
f_a = lambda x: np.prod(1.0 + a * (x - 0.5), axis=1)
f_b = lambda x: np.prod(1.0 + gam * 2 * np.pi**2 * b2(x), axis=1)
rng = np.random.default_rng(20261008)
print("random shifts (M = 32, N = 4096)")
print("        mean - 1     std.err    |mean - 1|/se   rms of single-shift error")
for f, label in ((f_a, "f_A"), (f_b, "f_B")):
    q = np.array([f((pts + rng.random(d)) % 1.0).mean() for _ in range(M)])
    se = q.std(ddof=1) / np.sqrt(M)
    print(f"  {label}  {q.mean() - 1:+.3e}   {se:.3e}   {abs(q.mean() - 1) / se:6.2f}        {np.sqrt(np.mean((q - 1) ** 2)):.3e}")
```

```text
Fibonacci lattice (d = 2): P_2 = (1/N) sum prod (1 + 2 pi^2 B_2) - 1
     N      P_2        N^2 P_2 / ln N
     89   1.603e-02      28.29
    377   1.139e-03      27.28
   1597   7.713e-05      26.67
   6765   5.060e-06      26.25
  28657   3.244e-07      25.96
 121393   2.044e-08      25.73
random shifts (M = 32, N = 4096)
        mean - 1     std.err    |mean - 1|/se   rms of single-shift error
  f_A  +3.797e-06   2.195e-05     0.17        1.223e-04
  f_B  +2.932e-07   4.140e-06     0.07        2.305e-05
```

Fibonacci 格子では $N^2P_2/\ln N$ が $26$〜$28$ でほぼ一定で、$P_2 = O(N^{-2}\log N)$ が見えます（<Ref to="ex-fibonacci-lattice" />）。$N$ を $1364$ 倍にすると誤差は $8\times10^5$ 分の 1 になります。

無作為な平行移動では、32 個の平均の真値からのずれが、32 個のばらつきから求めた標準誤差の $0.17$ 倍と $0.07$ 倍でした（どちらも 1 倍以内です。2 例だけなので、標準誤差が正しい大きさであることの確認としては弱く、矛盾がないことを見たにとどまります）。真値を知らなくても誤差の大きさが分かるというのが、<Ref to="prop-random-shift" /> の内容です。全体の関数評価は $32 \times 4096 \approx 1.3\times10^5$ 回で、同じ回数の Monte Carlo 法の標準誤差は $f_A$ で約 $9.8\times10^{-4}$（§8.2 の $N = 2^{16}$ の値の $1/\sqrt2$）なので、$f_A$ でも約 $45$ 分の 1 です。$f_B$ では、平行移動しない格子則の誤差 $1.68\times10^{-4}$（§8.2 の $N = 2^{12}$）より、平行移動 1 個あたりの誤差の 2 乗平均の平方根 $2.3\times10^{-5}$ のほうが小さくなっています。$f_B$ の Fourier 係数がすべて正なので、平行移動なし（$U = 0$）が双対格子の上の係数の位相がそろう最悪の場合に当たるためです。

## 9. 演習

<Exercise id="exr-equidistant-discrepancy" difficulty="易">
$[0,1)$ の $N$ 点 $\{0, 1/N, \ldots, (N-1)/N\}$ の $D_N^*$ を求めよ。これから、基数 $2$ の van der Corput 列の最初の $N = 2^m$ 項について $ND_N^* = 1$ であることを示せ。
<Solution>
$x_{(i)} = (i-1)/N$ なので $|x_{(i)} - (2i-1)/(2N)| = 1/(2N)$ で、<Ref to="prop-discrepancy-1d" /> から $D_N^* = 1/(2N) + 1/(2N) = 1/N$ です。van der Corput 列の最初の $2^m$ 項は、<Ref to="lem-elementary-interval" /> の証明の式で $q = 0$ とした $\varphi_2(r) = \rho(r)/2^m$（$0 \le r < 2^m$）で、$\rho$ が全単射なので集合として $\{c/2^m : 0 \le c < 2^m\}$ に等しく、$ND_N^* = 1$ です。§8.1 の表で $N$ が $2$ の冪でないときの値が $1$ より大きいのは、最後のブロックが途中で切れているためです。
</Solution>
</Exercise>

<Exercise id="exr-zaremba-check" difficulty="標準">
$f(x) = x^2$ と点集合 $P = \{0, 1/N, \ldots, (N-1)/N\}$ について、<Ref to="thm-koksma" /> の恒等式の両辺を直接計算して一致を確かめ、不等式の上界と比べよ。
<Solution>
左辺は $N^{-1}\sum_{i=0}^{N-1}(i/N)^2 - 1/3 = \dfrac{(N-1)(2N-1)}{6N^2} - \dfrac13 = -\dfrac{1}{2N} + \dfrac{1}{6N^2}$ です。右辺について、$i/N < t \le (i+1)/N$ で $A([0,t)) = i + 1$ なので $\delta_P(t) = (i+1)/N - t$ で、
$$
-\int_0^1\delta_P(t)\,2t\,dt = -\sum_{i=0}^{N-1}\int_{i/N}^{(i+1)/N}2t\Bigl(\frac{i+1}{N} - t\Bigr)dt
$$
です。$t = i/N + s$ とおくと、被積分関数は $2(i/N + s)(1/N - s)$ で、$\int_0^{1/N}(1/N - s)\,ds = 1/(2N^2)$、$\int_0^{1/N}s(1/N - s)\,ds = 1/(6N^3)$ から、各区間の積分は $i/N^3 + 1/(3N^3)$ です。$i$ について足すと $(N-1)/(2N^2) + 1/(3N^2) = 1/(2N) - 1/(6N^2)$ で、符号を付けると左辺に一致します。上界は $D_N^*\int_0^1|f'| = (1/N)\cdot1 = 1/N$ で、実際の誤差の約 2 倍です。$\delta_P \ge 0$ で $f' \ge 0$ なので打ち消しがなく、上界に近い値になっています。
</Solution>
</Exercise>

<Exercise id="exr-small-lattice" difficulty="標準">
$d = 2$、$N = 5$、$g = (1, 2)$ の格子則について、双対格子の $0$ でない点のうち $r(h)$ が最小のものを求めよ。$f_1(x) = \cos(2\pi(x_1 + x_2))$ と $f_2(x) = \cos(2\pi(x_1 + 2x_2))$ のそれぞれについて $Q_5(f)$ を求めよ。
<Solution>
双対格子は $h_1 + 2h_2 \equiv 0 \pmod 5$ です。$r(h) = 1$ の点 $(\pm1, 0)$、$(0, \pm1)$、$(\pm1, \pm1)$ は $h\cdot g = \pm1, \pm2, \pm3, \pm1$ で条件を満たしません。$r(h) = 2$ の点のうち $(1, 2)$ は $1 + 4 = 5$、$(2, -1)$ は $2 - 2 = 0$ で条件を満たします（符号を反転した点も）。よって最小値は $2$ です。

$f_1 = (e(h\cdot x) + e(-h\cdot x))/2$、$h = (1, 1)$ で、$\pm h \notin L^\perp$ なので <Ref to="thm-lattice-rule-error" /> の 1 から $Q_5(f_1) = \int f_1 = 0$（厳密）です。$f_2$ は $h = (1, 2) \in L^\perp$ で、$Q_5(f_2) = \hat f_2(h) + \hat f_2(-h) = 1/2 + 1/2 = 1$ です。実際、点 $x_n = (n/5, \{2n/5\})$ では $x_1 + 2x_2 \equiv 5n/5 = n \equiv 0 \pmod 1$ なので、$f_2$ はすべての点で $1$ です。積分は $0$ なので誤差は $1$ で、格子則は双対格子の上の振動数を定数と区別できません。
</Solution>
</Exercise>

<Exercise id="exr-anova-xy" difficulty="標準">
$f(x_1, x_2) = x_1x_2$ の ANOVA 分解を定義に従って求め、$\sigma^2 = \sigma_{\{1\}}^2 + \sigma_{\{2\}}^2 + \sigma_{\{1,2\}}^2$ を確かめよ。
<Solution>
$f_\emptyset = I = 1/4$。$\Pi_{\{1\}}f = x_1/2$ から $f_{\{1\}} = x_1/2 - 1/4 = (x_1 - 1/2)/2$、同様に $f_{\{2\}} = (x_2 - 1/2)/2$。
$$
f_{\{1,2\}} = x_1x_2 - \frac14 - \frac{x_1 - 1/2}{2} - \frac{x_2 - 1/2}{2} = \Bigl(x_1 - \frac12\Bigr)\Bigl(x_2 - \frac12\Bigr)
$$
です。$\int_0^1(x - 1/2)^2dx = 1/12$ から $\sigma_{\{1\}}^2 = \sigma_{\{2\}}^2 = 1/48$、$\sigma_{\{1,2\}}^2 = 1/144$ です。一方 $\sigma^2 = \int x_1^2x_2^2 - 1/16 = 1/9 - 1/16 = 7/144$ で、$3/144 + 3/144 + 1/144 = 7/144$ と一致します。1 変数の成分が分散の $6/7 \approx 86\%$ を占めます。$x_1x_2 = (1/2 + (x_1 - 1/2))(1/2 + (x_2 - 1/2))$ と書いて展開しても同じ分解が得られます。<Ref to="ex-product-anova" /> はこの形の一般化です。
</Solution>
</Exercise>

## 参考文献

1. L. Kuipers, H. Niederreiter, *Uniform Distribution of Sequences*, Wiley, 1974 — 第 1 章（Weyl の判定法）、第 2 章（食い違い量、Koksma–Hlawka の不等式、van der Corput 列・Halton 列）。
2. H. Niederreiter, *Random Number Generation and Quasi-Monte Carlo Methods*, SIAM, 1992 — 第 2〜5 章（食い違い量、低食い違い量列、ネット、格子則）。
3. I. H. Sloan, S. Joe, *Lattice Methods for Multiple Integration*, Oxford University Press, 1994 — 格子則の誤差、$P_\alpha$、良い格子の存在。
4. J. Dick, F. Y. Kuo, I. H. Sloan, "High-dimensional integration: the quasi-Monte Carlo way", *Acta Numerica* 22 (2013), 133–288 — 無作為化、重みつき空間、成分ごとの構成法、有効次元。
5. K. F. Roth, "On irregularities of distribution", *Mathematika* 1 (1954), 73–79 — 食い違い量の下界。
6. A. B. Owen, "Scrambled net variance for integrals of smooth functions", *The Annals of Statistics* 25 (1997), 1541–1562 — スクランブルしたネットの分散。