How to assess the convergence of MCMC algorithms

1. Tổng quan

Markov Chain Monte Carlo (MCMC) là một lớp các thuật toán được dùng khá rộng rãi trong ML và thống kê. Trong ngữ cảnh ML thì:

  • đôi khi ta cần ước lượng một phân bố xác suất nào đó, nhưng lại không thể tính chính xác phân bố đó do công thức quá phức tạp (thông thường vì trong công thức có phép tích phân), và
  • Ta lại biết rằng một chuỗi Markov, khi chạy đủ lâu, thì sẽ hội tụ về một phân bố xác suất cố định, độc lập với trạng thái ban đầu của nó (đây là tính chất đẹp và rất dễ chứng minh, ta sẽ trở lại khi có dịp). Phân bố xác suất này gọi là stationary distribution.

Như vậy tư tưởng của MCMC là cho trước một phân bố xác suất cần ước lượng, ta sẽ xây dựng một chuỗi Markov sao cho stationary distribution của nó chính là phân bố xác suất đang xét. Xuất phát từ trạng thái ngẫu nhiên nào đó, chạy chuỗi Markov một thời gian đủ lâu để nó hội tụ về stationary distribution, thì khi đó mỗi trạng thái của chuỗi Markov có thể xem là được lấy mẫu từ phân bố xác suất đang cần ước lượng.

Ví dụ 1: Một ví dụ kinh điển là các thuật toán suy diễn của undirected graphical model (còn gọi là Markov network). Thông thường joint distribution của các biến trong Markov network được định nghĩa là

\displaystyle P\left(X_1,X_2,\ldots ,X_n\right)=\frac{1}{Z}\phi_1\left(X_{\phi_1}\right)\ldots\phi_k\left(X_{\phi_k}\right)

Trong đó Z thường được gọi là partition function, các \phi_k là potential thứ k thể hiện mối liên hệ giữa các biến nằm trong tập X_{\phi_{k}}\subset\left\{ X_{1},\ldots,X_{n}\right\}X_{\phi_k} thường được gọi là scope của \phi_k. Đặc biệt partition function Z định nghĩa là:

\displaystyle Z =\sum_{X_1,X_2,\ldots ,X_n}\phi_1\left(X_{\phi_1}\right)\ldots\phi_k\left(X_{\phi_k}\right)

nghĩa là ta phải tính tổng trên toàn bộ tất cả các giá trị có thể nhận của n biến trong Markov network. Giả sử các biến là nhị phân (chỉ nhận giá trị 0 hay 1), thì n biến trong mạng có thể tạo thành không gian 2^n phần tử (exponential theo số lượng biến). Rõ ràng việc duyệt qua tất cả các phần tử trong không gian này là bất khả thi khi n tương đối lớn. Như vậy ngay cả việc tính xác suất cho một cấu hình bất kì của n biến cũng có vẻ như là bất khả thi? Thực ra không phải như vậy. Bằng thuật toán đơn giản kiểu như Variable Elimination, ta đã có thể rút độ phức tạp của việc tính toán Z xuống thành đa thức, thay vì exponential, và từ đó hoàn toàn có thể tính xác suất cho một cấu hình bất kì.

Tuy nhiên vấn đề phức tạp là đôi khi ta muốn lấy mẫu (sampling) từ phân phối của Markov network. Rõ ràng là bất khả thi nếu tính xác suất của mọi cấu hình có thể, rồi lấy mẫu dựa trên phân bố xác suất này (vì phân bố xác suất này sẽ có 2^n giá trị trong trường hợp biến nhị phân). Đây là một trường hợp tiêu biểu của việc ứng dụng phương pháp MCMC.

Ví dụ 2: Một ví dụ nữa cũng khá nổi tiếng là trong tracking. Trong hệ thống (Bayesian) tracking, thông thường việc phán đoán trạng thái của đối tượng cần theo vết được thực hiện dựa vào phân bố xác suất hậu nghiệm (posterior distribution) được tính dựa trên mô hình chuyển động và mô hình likelihood. Trong nhiều trường hợp (với nhiều ràng buộc/hạn chế), phân bố xác suất hậu nghiệm là đơn giản và ta có thể sử dụng Kalman filter để thực hiện theo vết. Trong nhiều trường hợp khác, phân bố xác suất hậu nghiệm phức tạp hơn nhiều, do đó cách duy nhất khả dĩ để thực hiện theo vết là sử dụng phương pháp MCMC hoặc một vài biến thể của nó. (Tracking cũng là một bài toán lớn trong Computer Vision mà ta sẽ trở lại khi có dịp).

2. Các thành phần của một thuật toán MCMC

Như vậy ta đã thấy MCMC cung cấp một phương pháp tổng quát để lấy mẫu từ phân bố xác suất bất kì mà thông thường rất khó ước lượng trực tiếp. Nói chung để sử dụng phương pháp MCMC thì các thành phần sau đây cần được xác định:

  • Xây dựng chuỗi Markov sao cho nó có stationary distribution là phân bố xác suất đang xét. Thông thường các trạng thái của chuỗi Markov chính là các trạng thái khả dĩ trong phân bố xác suất đang xét. Do đó việc còn lại cần làm chỉ là tính toán xác suất dịch chuyển (transition matrix) từ trạng thái này sang trạng thái khác trong chuỗi Markov. Hai phương pháp thường dùng nhất là:
    • Gibbs sampling: đơn giản và dễ dùng, nhất là khi phân bố xác suất có điều kiện P\left(X_i \vert \boldsymbol{X}_{-i}\right) có thể được tính dễ dàng.
    • Metropolis-Hastings: phương pháp tương đối tổng quát để thiết kế transition model cho chuỗi Markov, dựa trên một proposal distribution và acceptance ratio để điều khiển dynamic của chuỗi Markov.
  • Chọn lựa chiến lược lấy mẫu từ chuỗi Markov. Trong thực tế các chiến lược lấy mẫu sau hay dùng:
    • Chạy 1 chuỗi Markov từ vị trí ban đầu ngẫu nhiên nào đó, sau khoảng T bước với T đủ lớn (gọi là burn-in phase), ta lấy 1 mẫu. Sau đó chạy lại chuỗi Markov từ vị trí ngẫu nhiên khác, sau T bước lại lấy thêm 1 mẫu. Cứ như vậy M lần ta sẽ được M mẫu độc lập từ phân bố xác suất đang xét (ta gọi là i.i.d samples).
      Tuy nhiên cách làm này khá tốn kém vì thông thường T rất lớn, trong khi ta lại chỉ lấy 1 mẫu sau mỗi bước như vậy. Nếu quả thật sau burn-in phase, chuỗi Markov đã hội tụ về stationary distribution, thì việc chỉ lấy 1 mẫu là khá lãng phí, vì về nguyên tắc ta có thể lấy thêm vài mẫu nữa mà vẫn đảm bảo là lấy từ phân bố xác suất đang xét.
    • Chỉ chạy 1 chuỗi Markov, sau burn-in phase, lấy các mẫu T+1, T+2,...T+M. Tuy nhiên làm theo cách này, rõ ràng M mẫu lấy được là không độc lập.
    • Chỉ chạy 1 chuỗi Markov, sau burn-in phase, lấy mẫu T+d, T+2d, T+3d..., T+Md, với d\in \mathbb{N} nào đó. Cách này vẫn không cho ta M mẫu độc lập, và về lí thuyết thì cũng không hơn gì cách lấy M mẫu liên tiếp ở trên.
    • Chạy một số nhỏ K chuỗi Markov song song, sau đó lấy mẫu đồng thời  từ K chuỗi. Bằng cách này, ta vẫn không đảm bảo là được M mẫu độc lập, nhưng có vẻ như có cơ sở để tin rằng M mẫu lấy được sẽ ít correlated hơn các cách ở trên. Hơn nữa bằng cách này, ta có thể thực hiện một số thống kê để xét tính hội tụ của các chuỗi Markov. Đây là cách thường được dùng trong thực tế.
  • Cuối cùng là cần chọn một chiến lược để đánh giá tính hội tụ của phương pháp MCMC (thực chất là xét xem các chuỗi Markov có hội tụ hay chưa).

Mục tiêu của bài này là nói về việc đánh giá tính hội tụ của phương pháp MCMC. Các chủ đề khác như Gibbs sampling hay Metropolis-Hastings sẽ trở lại khi có dịp.

3. Đánh giá sự hội tụ của thuật toán MCMC

Có nhiều cách để đánh giá tính hội tụ của phương pháp MCMC. Tuy nhiên bài này, dựa theo [1], trình bày một cách khá đơn giản nhưng hiệu quả cho việc này. Ý tưởng chính của các phương pháp này là tính variance của tập mẫu lấy được, và xét tính hội tụ dựa trên variance của tập mẫu.

Một cách hình thức, giả sử ta chạy K chuỗi Markov song song, bắt đầu từ K trạng thái ban đầu tương đối khác nhau. Sau T bước burn-in phase (chạy các chuỗi Markov T bước và không lấy mẫu nào), gọi X_k^m là một mẫu lấy từ chuỗi k tại bước thứ T+m, giả sử rằng ta chỉ lấy M mẫu từ T+1 đến T+M của mỗi chuỗi. Khi đó ta có thể tính variance giữa các chuỗi (B) và trong mỗi chuỗi (W) như sau:

\begin{aligned}\bar{f}_{k} & =\frac{1}{M}\sum_{m=1}^{M}f\left(\mathbf{X}_{k}^{m}\right)\\  \bar{f} & =\frac{1}{K}\sum_{k=1}^{K}\bar{f}_{k}\\  B & =\frac{M}{K-1}\sum_{k=1}^{K}\left(\bar{f}_{k}-\bar{f}\right)^{2}\\  W & =\frac{1}{K}\frac{1}{M-1}\sum_{k=1}^{K}\sum_{m=1}^{M}\left(f\left(\mathbf{X}_{k}^{m}\right)-\bar{f}_{k}\right)^{2}\\  V & =\frac{M-1}{M}W+\frac{1}{M}B  \end{aligned}

Trong đó f\left(\cdot\right) là hàm bất kì (gọi là statistic) có thể tính trên mỗi mẫu X_k^m.

Đến lúc này ta có thể dùng ước lượng \hat{R} = \sqrt{\frac{V}{W}} để đo mức độ tương đồng giữa các chuỗi Markov. Cụ thể \hat{R} có giá trị lớn khi các chuỗi chưa hội tụ về cùng một phân phối xác suất. Nếu \hat{R}=1, điều này có nghĩa là hoặc  các chuỗi đã hội tụ về phân phối đang xét, hoặc là tập K điểm khởi đầu  không đủ đa dạng nên cả K chuỗi hội tụ về cùng một vùng local nào đó trong không gian đang xét. Điều đó có nghĩa là \hat{R} > 1 thì chắc chắn các chuỗi chưa hội tụ, nhưng \hat{R}=1 cũng không có nghĩa là các chuỗi đã hội tụ về phân phối xác suất đang xét.

Ta có thể xem xét kĩ hơn công thức của \hat{R}:

\displaystyle \hat{R}=\sqrt{\frac{V}{W}}=\sqrt{1-\frac{1}{M}\left(1-\frac{B}{W}\right)}

như vậy \hat{R}\approx 1 khi 1) M đủ lớn, hoặc 2) B\approx W. Rõ ràng điều này có ý nghĩa. Nếu M đủ lớn thì tập mẫu của ta đủ lớn để ước lượng phân phối xác suất đang xét, và khi B\approx W thì variance giữa các chuỗi (B) xấp xỉ bằng variance trong mỗi chuối (W), điều này có nghĩa là các chuỗi đã tiến gần về với nhau. Nếu cả 2 điều này cùng đúng thì \hat{R} càng gần về 1.

Trong thực  tế, người ta có thể dùng nhiều hàm f\left(\cdot\right) khác nhau để tăng mức độ confident về tính hội tụ của các chuỗi.

4. Ứng dụng

Trong một đồ án gần đây, người viết bài này đã cài đặt phương pháp này để lấy mẫu từ một Markov network có n=12 biến. Sử dụng Gibbs sampling, chạy đồng thời K=4 chuỗi Markov và lấy M=200 mẫu. Cách làm là cứ sau mỗi M=200 mẫu, ta tính \hat{R} và vẽ lại đồ thị, như trong hình dưới đây.

Giá trị của \hat{R} qua quá trình chạy Gibbs sampling với hơn 2000×200 = 400.000 lần lấy mẫu. Hình này ứng với hàm f\left(\cdot\right) là khoảng cách Euclidean từ X_k^m đến mean của tập M mẫu tương ứng

Giá trị của \hat{R} qua quá trình chạy Gibbs sampling với hơn 2000×200 = 400.000 lần lấy mẫu. Hình này ứng với hàm f\left(\cdot\right) là xác suất của X_k^m

Ta có thể thấy rằng đến hơn 400 nghìn lần lấy mẫu nhưng giá trị của \hat{R} vẫn dao động khá lớn quanh 1. Nhìn kĩ một chút thì \hat{R} dao động trong khoảng 0.998 đến 1.01, nhưng cũng rất khó để nói là các chuỗi đã hội tụ.

Việc này có thể giải thích là bản chất của thuật toán Gibbs sampling khiến cho ta di chuyển trong không gian từng bước rất nhỏ, do đó có thể các chuỗi Markov đều kẹt trong vùng cục bộ nào đó và không thể cùng hội tụ.

5. Kết luận

Bài này không có ý định trình bày chi tiết về MCMC và các thuật toán liên quan, mà chỉ thảo luận một vấn đề nhỏ (nhưng vô cùng quan trọng) trong việc áp dụng phương pháp MCMC, đó là về tính hội tụ. Phần này chủ yếu lấy cảm hứng từ một vài footnote trong [1], nhưng thực sự vẫn còn là vấn đề khó và hình như chưa có câu trả lời nào trọn vẹn.

Một ghi chú ngoài lề, trong trường hợp ta có thể tính được phân bố xác suất tại một điểm bất kì (nhưng không thể tính cho toàn miền giá trị vì chi phí quá cao), thì độ đo relative entropy cũng là một cách để đo mức độ tương đồng (gần gũi) giữa phân phối xác suất thực và phân phối xác suất của tập mẫu (gọi là empirical distribution).

Tham khảo

[1] D. Koller and N. Friedman, Probabilistic Graphical Models: Principles and Techniques. The MIT Press, 2009.

Advertisements

2 comments

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s