# Chuỗi Sobol và Quasi-Monte Carlo


## 1. Giới thiệu

Trong [Mô phỏng Monte Carlo](/posts/math/stats/monte-carlo-simulation/), chúng ta thường sử dụng các số giả ngẫu nhiên (pseudo-random numbers) để tính toán tích phân số học hoặc mô phỏng các hệ thống phức tạp. Tuy nhiên, các bộ sinh số giả ngẫu nhiên truyền thống có một nhược điểm lớn: các điểm sinh ra thường có xu hướng tạo thành các **cụm (clumps)** và để lại những **khoảng trống (gaps)** trong không gian. 

Điều này dẫn đến việc ước lượng không hiệu quả, vì một số khu vực bị lấy mẫu quá nhiều, trong khi một số khu vực khác lại bị bỏ qua hoàn toàn. Để giải quyết vấn đề này, phương pháp **Quasi-Monte Carlo (QMC)** ra đời, thay thế số giả ngẫu nhiên bằng các chuỗi số được thiết kế đặc biệt để phân bố đều đặn hơn.

## 2. Chuỗi độ sai lệch thấp (Low-discrepancy sequence)

Trái tim của phương pháp QMC chính là các **chuỗi độ sai lệch thấp (low-discrepancy sequences)**. "Sai lệch" (discrepancy) ở đây là một độ đo toán học về mức độ không đồng đều của các điểm trong một không gian. 

Một chuỗi có độ sai lệch thấp không thực sự "ngẫu nhiên". Thay vào đó, mục tiêu của nó là lấp đầy không gian một cách **đồng đều nhất có thể**. Nếu bạn tưởng tượng việc rắc hạt tiêu lên một quả trứng ốp la, số giả ngẫu nhiên sẽ tạo ra những chỗ cay xè và những chỗ nhạt nhẽo, trong khi chuỗi độ sai lệch thấp sẽ đảm bảo mỗi miếng cắn đều có vị cay vừa vặn.

### 2.1. Trước hết: bạn cần loại điểm nào?

PCG64DXSM và scrambled Sobol đều có thể trả về các số trong $[0,1)$, nhưng chúng tối ưu cho hai yêu cầu khác nhau:

<link rel="stylesheet" href="/css/rng-diagrams.css">

<figure class="rng-diagram rng-decision" aria-label="Sơ đồ chọn giữa PCG64DXSM và scrambled Sobol">
  <div class="rng-diagram__node rng-decision__question">Bạn muốn gì?</div>
  <div class="rng-diagram__arrow" aria-hidden="true"><span class="rng-diagram__arrow-symbol">↓</span></div>
  <div class="rng-decision__branches">
    <section class="rng-decision__branch rng-decision__branch--iid">
      <div class="rng-diagram__criterion">Các số vẫn ngẫu nhiên độc lập,<br>nhưng về lâu dài tuân theo Uniform(0,1)</div>
      <div class="rng-diagram__arrow" aria-hidden="true"><span class="rng-diagram__arrow-symbol">↓</span></div>
      <div class="rng-diagram__answer">PCG64DXSM</div>
    </section>
    <section class="rng-decision__branch rng-decision__branch--qmc">
      <div class="rng-diagram__criterion">Một tập hữu hạn các điểm phủ không gian<br>càng đều, càng ít khoảng trống càng tốt</div>
      <div class="rng-diagram__arrow" aria-hidden="true"><span class="rng-diagram__arrow-symbol">↓</span></div>
      <div class="rng-diagram__answer">Scrambled Sobol</div>
    </section>
  </div>
</figure>

- Nếu mô hình cần các lần rút gần với giả thiết **độc lập và cùng phân phối** (independent and identically distributed, IID), hãy dùng một PRNG hiện đại như [PCG64DXSM](/posts/math/stats/pcg-to-pcg64dxsm/).
- Nếu mục tiêu là xấp xỉ tích phân bằng một **tập hữu hạn điểm phủ miền thật đều**, scrambled Sobol thường thích hợp hơn. Các điểm của nó có phụ thuộc được thiết kế có chủ ý; chính sự phụ thuộc này giúp giảm hiện tượng vón cụm và khoảng trống.

Xáo trộn không biến Sobol thành một chuỗi IID. Nó ngẫu nhiên hoá mạng điểm để ta có thể lặp thí nghiệm và ước lượng sai số, trong khi vẫn giữ cấu trúc phân tầng của Quasi-Monte Carlo.

## 3. Chuỗi Sobol

Được nhà toán học người Nga Ilya M. Sobol giới thiệu trong một [công trình năm 1967](https://www.mathnet.ru/eng/zvmmf7334), **chuỗi Sobol (Sobol sequence)** là một trong những chuỗi độ sai lệch thấp phổ biến nhất.

Điểm khác biệt cốt lõi không nằm ở việc Sobol dùng phép toán bit còn bộ sinh số giả ngẫu nhiên (PRNG) dùng phép cộng hay phép nhân — nhiều PRNG hiện đại cũng dùng phép XOR. Khác biệt nằm ở **mục tiêu thiết kế**: PRNG cố tái hiện các đặc trưng thống kê của những mẫu ngẫu nhiên độc lập, còn Sobol chủ động sắp điểm theo một cấu trúc nhị phân để những khối đầu tiên gồm $2^m$ điểm phủ không gian cân bằng.

Ta có thể hình dung mỗi điểm mới đang đi vào một vùng còn thưa. Tuy nhiên, thuật toán không thực sự tìm “khoảng trống lớn nhất” sau mỗi bước. Vị trí của điểm đã được quyết định trước bởi ba thành phần: biểu diễn nhị phân của chỉ số, **mã Gray (Gray code)** và các **số định hướng (direction numbers)**.

### 3.1. Bắt đầu từ một chiều

Trên đoạn $[0,1)$, tám điểm đầu tiên của chuỗi Sobol là

$$
0,\quad \frac12,\quad \frac34,\quad \frac14,\quad
\frac38,\quad \frac78,\quad \frac58,\quad \frac18.
$$

Thứ tự này thoạt nhìn hơi lạ, nhưng có một quy luật rõ ràng. Hai điểm đầu rơi vào hai nửa của đoạn; bốn điểm đầu rơi vào bốn khoảng dài $1/4$; tám điểm đầu rơi vào tám khoảng dài $1/8$. Tổng quát hơn, $2^m$ điểm đầu tiên cho đúng một điểm trong mỗi khoảng nhị phân

$$
\left[\frac{k}{2^m},\frac{k+1}{2^m}\right),
\qquad k=0,1,\ldots,2^m-1.
$$

Để sinh điểm thứ $n$, trước hết ta đổi chỉ số $n$ sang mã Gray:

$$
g(n)=n\mathbin{\mathtt{XOR}}\left\lfloor\frac n2\right\rfloor.
$$

Hai mã Gray liên tiếp chỉ khác nhau đúng một bit. Trong một chiều, các số định hướng là

$$
v_1=0.1_2=\frac12,\quad
v_2=0.01_2=\frac14,\quad
v_3=0.001_2=\frac18,\quad\ldots
$$

Nếu $g_k(n)$ là bit thứ $k$ của mã Gray, tính từ phải sang trái, thì

$$
x_n=(g_1v_1)\mathbin{\mathtt{XOR}}(g_2v_2)
\mathbin{\mathtt{XOR}}\cdots.
$$

Phép XOR ở đây tác động trên các chữ số nhị phân và không có nhớ. Chẳng hạn, $n=4=100_2$ cho $g(4)=110_2$. Hai bit thứ hai và thứ ba đang bật, nên

$$
x_4=v_2\mathbin{\mathtt{XOR}}v_3
=0.010_2\mathbin{\mathtt{XOR}}0.001_2
=0.011_2=\frac38.
$$

Vì mã Gray chỉ đổi một bit ở mỗi bước, ta không cần tính lại toàn bộ tổng XOR. Muốn đi từ $x_{n-1}$ tới $x_n$, thuật toán chỉ cần XOR điểm cũ với đúng một số định hướng. Đây là lý do chuỗi Sobol có thể được sinh rất nhanh.

### 3.2. Từ một chiều lên hai chiều

Trong hai chiều, ta vẫn dùng cùng mã Gray của chỉ số $n$, nhưng mỗi trục có một họ số định hướng riêng:

$$
x_{n,j}=\bigoplus_{k\geq 1} g_k(n)v_{j,k},
\qquad j=1,2.
$$

Với hai chiều đầu tiên, vài số định hướng là

| $k$ | $v_{1,k}$ cho trục $x$ | $v_{2,k}$ cho trục $y$ |
|---:|:---:|:---:|
| 1 | $1/2$ | $1/2$ |
| 2 | $1/4$ | $3/4$ |
| 3 | $1/8$ | $5/8$ |

Chẳng hạn, $n=2$ có mã Gray $011_2$, nên

$$
(x_{2,1},x_{2,2})
=(v_{1,1}\mathbin{\mathtt{XOR}}v_{1,2},
  v_{2,1}\mathbin{\mathtt{XOR}}v_{2,2})
=\left(\frac34,\frac14\right).
$$

Nhờ trục $y$ không lặp lại nguyên xi các số định hướng của trục $x$, các điểm không nằm trên một đường chéo. Với cấu hình hai chiều minh hoạ dưới đây, $2^m$ điểm đầu cho đúng một điểm trong mỗi hình chữ nhật nhị phân có diện tích $2^{-m}$. Ở số chiều cao hơn, việc chọn số định hướng trở nên quan trọng; các thư viện thực tế thường dùng những bộ đã được tìm kiếm và kiểm nghiệm, chẳng hạn [bộ số định hướng Joe–Kuo](https://web.maths.unsw.edu.au/~fkuo/sobol/).

<link rel="stylesheet" href="/css/sobol-explorer.css">

<figure class="sobol-explorer" data-sobol-construction>
  <div class="sobol-explorer__header">
    <p class="sobol-explorer__eyebrow">Minh hoạ tương tác</p>
    <h4>Cách chuỗi Sobol thêm từng điểm</h4>
    <p>Đổi số chiều hoặc thêm từng điểm để quan sát cấu trúc nhị phân được hoàn thiện dần.</p>
  </div>
  <div class="sobol-explorer__controls">
    <div class="sobol-segment" role="group" aria-label="Chọn số chiều">
      <button type="button" data-dimension="1" aria-pressed="false">1D</button>
      <button type="button" data-dimension="2" aria-pressed="true">2D</button>
    </div>
    <label class="sobol-range">
      <span>Số điểm: <strong data-count-output>8</strong></span>
      <input type="range" min="1" max="32" value="8" step="1" data-count-input>
    </label>
    <div class="sobol-actions" role="group" aria-label="Điều khiển chuỗi điểm">
      <button type="button" data-action="previous" aria-label="Bớt một điểm">−1</button>
      <button type="button" data-action="next">+1 điểm</button>
      <button type="button" data-action="play">Tự chạy</button>
      <button type="button" data-action="reset">Về đầu</button>
    </div>
  </div>
  <div class="sobol-explorer__stage">
    <canvas data-construction-canvas role="img" aria-label="Minh hoạ các điểm Sobol trong hai chiều"></canvas>
  </div>
  <div class="sobol-readouts">
    <div><span>Chỉ số và mã Gray</span><strong data-gray-output></strong></div>
    <div><span>Bit vừa đổi</span><strong data-direction-output></strong></div>
    <div><span>Điểm vừa thêm</span><strong data-point-output></strong></div>
  </div>
  <figcaption data-construction-caption aria-live="polite"></figcaption>
  <noscript>Cần bật JavaScript để dùng minh hoạ này.</noscript>
</figure>

Minh hoạ dùng các số định hướng chuẩn cho hai chiều đầu. Hãy dừng ở $2$, $4$, $8$, $16$ hoặc $32$ điểm: đó là lúc một tầng phân hoạch nhị phân vừa được lấp đầy. Ở các bước nằm giữa hai lũy thừa của $2$, tầng mới vẫn đang trong quá trình hoàn thiện.

### 3.3. Xáo trộn chuỗi Sobol

Chuỗi Sobol nguyên bản hoàn toàn tiền định: cùng số chiều và cùng bộ số định hướng luôn cho cùng một dãy điểm. Điều này tốt cho khả năng tái lập, nhưng một ước lượng duy nhất không cung cấp các lần lặp ngẫu nhiên độc lập để đánh giá sai số. **Xáo trộn (scrambling)** đưa ngẫu nhiên vào vị trí các điểm mà vẫn giữ cấu trúc phân tầng của mạng Sobol.

Cách đơn giản nhất là **dịch số nhị phân (digital shift)**. Với mỗi chiều $j$, ta sinh một mặt nạ bit ngẫu nhiên $\Delta_j$ và đặt

$$
\widetilde{x}_{n,j}=x_{n,j}\mathbin{\mathtt{XOR}}\Delta_j.
$$

Đây không phải phép cộng một độ dời hình học thông thường. XOR hoán vị các ô nhị phân: mỗi ô cũ đi tới đúng một ô mới, nên số điểm trong các ô vẫn được bảo toàn.

**Owen scrambling** đi xa hơn. Hãy xem các chữ số sau dấu phẩy nhị phân như một cây: bit đầu chọn nửa trái hay nửa phải; bit thứ hai chọn một trong hai phần tư bên trong nửa ấy; rồi đến phần tám, phần mười sáu, v.v. Ở mỗi nút cây, Owen scrambling chọn ngẫu nhiên giữ nguyên hay đổi chỗ hai nhánh. Phép đổi ở một tầng có thể phụ thuộc vào các bit đứng trước, vì vậy hai điểm ở hai nhánh khác nhau không nhất thiết nhận cùng một mặt nạ. Kết quả trông bớt cứng nhắc hơn phép dịch XOR nhưng vẫn bảo toàn sự cân bằng của mạng trên các ô nhị phân.

<figure class="sobol-explorer" data-sobol-scramble>
  <div class="sobol-explorer__header">
    <p class="sobol-explorer__eyebrow">Minh hoạ tương tác</p>
    <h4>Xáo trộn trong 1D và 2D</h4>
    <p>So sánh các điểm gốc với cùng một mạng sau khi xáo trộn. Mỗi seed tạo ra một phép hoán vị khác.</p>
  </div>
  <div class="sobol-explorer__controls">
    <div class="sobol-segment" role="group" aria-label="Chọn số chiều cho minh hoạ xáo trộn">
      <button type="button" data-scramble-dimension="1" aria-pressed="false">1D</button>
      <button type="button" data-scramble-dimension="2" aria-pressed="true">2D</button>
    </div>
    <label class="sobol-select">
      <span>Phép xáo trộn</span>
      <select data-scramble-method>
        <option value="shift">Dịch số XOR</option>
        <option value="owen" selected>Owen lồng</option>
      </select>
    </label>
    <label class="sobol-range sobol-range--short">
      <span>Số điểm: <strong data-scramble-count-output>8</strong></span>
      <input type="range" min="2" max="5" value="3" step="1" data-scramble-count-input>
    </label>
    <button class="sobol-seed-button" type="button" data-new-seed>Đổi seed</button>
  </div>
  <div class="sobol-explorer__stage">
    <canvas data-scramble-canvas role="img" aria-label="So sánh mạng Sobol gốc và mạng sau xáo trộn trong hai chiều"></canvas>
  </div>
  <div class="sobol-readouts">
    <div><span>Seed</span><strong data-seed-output></strong></div>
    <div><span>Phương pháp</span><strong data-method-output></strong></div>
    <div><span>Kiểm tra phân tầng</span><strong data-balance-output></strong></div>
  </div>
  <figcaption data-scramble-caption aria-live="polite"></figcaption>
  <noscript>Cần bật JavaScript để dùng minh hoạ này.</noscript>
</figure>

Trong 1D, các đường nối cho biết từng điểm đã đi từ vị trí gốc tới vị trí mới. Trong 2D, hai hình vuông cho thấy toàn bộ mạng trước và sau khi xáo trộn. Hãy đổi seed nhiều lần: hình dạng thay đổi, nhưng ở $N=2^m$ số điểm trong các ô nhị phân vẫn cân bằng.

Một chi tiết thực hành quan trọng: `scipy.stats.qmc.Sobol(scramble=True)` dùng **left linear matrix scramble kết hợp với digital random shift (LMS+shift)**, không phải Owen scrambling theo nghĩa hẹp. Cả hai đều thuộc họ Randomized Quasi-Monte Carlo (RQMC). Bằng cách chạy nhiều seed độc lập, ta có thể quan sát độ biến thiên giữa các lần chạy và xây dựng ước lượng sai số. Theo [tài liệu SciPy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.Sobol.html), nên lấy $N=2^m$ điểm bằng `random_base2(m)` để giữ các tính chất cân bằng; bỏ điểm đầu, lấy cách quãng hoặc dùng số điểm tùy ý có thể làm hỏng chúng.

<script defer src="/js/sobol-explorer.js"></script>

## 4. So sánh với PRNG truyền thống

Các bộ sinh giả ngẫu nhiên truyền thống như [LCG](/posts/math/stats/linear-congruential-generator/) cố làm mờ cấu trúc để đầu ra trông giống các mẫu ngẫu nhiên độc lập. Sobol đi theo hướng khác: nó chấp nhận một cấu trúc nhìn thấy được để đổi lấy độ phủ không gian cân bằng.

Sự khác biệt lớn nhất nằm ở **tốc độ hội tụ** khi ước lượng tích phân:

- Với Monte Carlo tiêu chuẩn, sai số căn trung bình thường giảm theo bậc $O(N^{-1/2})$. Muốn giảm sai số 10 lần, ta thường phải tăng số mẫu lên khoảng 100 lần.
- Với Quasi-Monte Carlo, cận sai số cho các hàm đủ đều có dạng gần $O(N^{-1})$, kèm theo một thừa số logarit phụ thuộc số chiều. Trong thực tế, lợi thế rõ nhất khi hàm đủ trơn và bài toán có **chiều hiệu dụng** không quá cao; đây không phải bảo đảm rằng mọi tích phân đều hội tụ đúng theo $1/N$.

{{< echarts >}}
{
  "title": { "text": "So sánh tốc độ hội tụ (Log-Log)", "left": "center" },
  "tooltip": { "trigger": "item" },
  "legend": { "data": ["Monte Carlo (O(1/√N))", "Quasi-Monte Carlo (O(1/N))"], "top": "bottom" },
  "xAxis": { "type": "log", "name": "Số điểm mẫu (N)", "nameLocation": "middle", "nameGap": 25 },
  "yAxis": { "type": "log", "name": "Sai số ước lượng", "nameLocation": "middle", "nameGap": 35 },
  "series": [
    {
      "name": "Monte Carlo (O(1/√N))",
      "type": "line",
      "showSymbol": false,
      "data": [[10, 0.3162], [100, 0.1], [1000, 0.0316], [10000, 0.01], [100000, 0.00316], [1000000, 0.001]]
    },
    {
      "name": "Quasi-Monte Carlo (O(1/N))",
      "type": "line",
      "showSymbol": false,
      "data": [[10, 0.1], [100, 0.01], [1000, 0.001], [10000, 0.0001], [100000, 0.00001], [1000000, 0.000001]]
    }
  ]
}
{{< /echarts >}}
*Biểu đồ so sánh mức độ giảm sai số theo số lượng điểm mẫu (thang logarit). Khoảng cách giữa hai đường ngày càng rộng ra khi số lượng mẫu (N) tăng, thể hiện ưu thế vượt trội của Quasi-Monte Carlo.*

Hiệu quả này cực kỳ rõ rệt khi số lượng mẫu lớn, giúp tiết kiệm đáng kể tài nguyên máy tính.

## 5. Ứng dụng

Với ưu điểm về tốc độ hội tụ, Quasi-Monte Carlo và chuỗi Sobol được ứng dụng rộng rãi trong:

- **Tài chính định lượng:** Định giá các công cụ phái sinh phức tạp (pricing derivatives), nơi cần tính toán tích phân ở số chiều rất cao.
- **Đồ họa máy tính (Rendering):** Lấy mẫu ánh sáng (ray tracing) để tạo ra hình ảnh mượt mà, ít nhiễu (noise) hơn so với việc dùng số ngẫu nhiên thuần túy.
- **Vật lý và Kỹ thuật:** Mô phỏng các hệ thống vật lý và tối ưu hóa hệ thống.

## 6. Ví dụ Python

Hãy cùng xem sự khác biệt giữa hàm sinh số ngẫu nhiên mặc định và chuỗi Sobol trong Python. Chúng ta sẽ sử dụng `scipy.stats.qmc.Sobol`.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import qmc

# Số điểm cần sinh
N = 256

# Sinh số ngẫu nhiên (Pseudo-random)
np.random.seed(42)
points_random = np.random.rand(N, 2)

# Sinh chuỗi Sobol đã xáo trộn (Randomized Quasi-Monte Carlo)
sampler = qmc.Sobol(d=2, scramble=True, seed=42)
points_sobol = sampler.random_base2(m=8)  # 2^8 = 256 điểm

# Vẽ đồ thị so sánh
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

ax1.scatter(points_random[:, 0], points_random[:, 1], alpha=0.6, s=15, color='blue')
ax1.set_title("Pseudo-Random (numpy.random)")
ax1.set_xlim(0, 1)
ax1.set_ylim(0, 1)
ax1.set_aspect('equal')

ax2.scatter(points_sobol[:, 0], points_sobol[:, 1], alpha=0.6, s=15, color='red')
ax2.set_title("Sobol Sequence (scipy.stats.qmc)")
ax2.set_xlim(0, 1)
ax2.set_ylim(0, 1)
ax2.set_aspect('equal')

plt.tight_layout()
plt.show()
```

Khi chạy đoạn code trên, bạn sẽ thấy biểu đồ bên trái (ngẫu nhiên) có các điểm tập trung thành từng cụm và để lại nhiều mảng trắng. Biểu đồ bên phải phủ hình vuông cân bằng hơn. Ở đây `random_base2(m=8)` chủ động lấy trọn một khối $2^8$ điểm; `scramble=True` dùng phép xáo trộn LMS+shift đã trình bày ở trên.

## 7. Tổng kết

Chuỗi Sobol và phương pháp Quasi-Monte Carlo là một bước tiến lớn giúp tối ưu hóa quá trình mô phỏng bằng cách thay thế sự "ngẫu nhiên" bằng sự "đồng đều". Khi đối mặt với các bài toán đòi hỏi tính toán tích phân nhiều chiều hay mô phỏng với chi phí tính toán cao, việc chuyển từ số giả ngẫu nhiên sang chuỗi độ sai lệch thấp như Sobol thường mang lại sự cải thiện hiệu suất vượt bậc.