Suy diễn Biến phân (Phần 2)

Authors

Affiliations

The-Anh Vu-Le

Published

March 13, 2024

Ở Phần 1, chúng ta đã xem xét một ví dụ đơn giản của suy diễn biến phân. Trong bài viết ngày hôm nay, chúng ta sẽ xem xét một ví dụ phức tạp hơn, đòi hỏi các giả định mạnh hơn để giải quyết.

Bối cảnh

Một bóng đèn có tuổi thọ thực tế $Z$ là một biến ngẫu nhiên tuân theo phân phối lũy thừa với hằng số $Λ$ . Khác với phần trước, ở đây $Λ$ là một biến ngẫu nhiên tuân theo phân phối Gamma với hằng số $α > 0$ và $β > 0$ . Khi sử dụng, tuổi thọ thực nghiệm $X$ của bóng đèn có tồn tại một nhiễu ngẩu nhiên với phân phối chuẩn (Gaussian) có kỳ vọng $0$ và phương sai $σ^{2} \in (0, \infty)$ . Giả sử ta có một bóng đèn với tuổi thọ thực nghiệm là $x > 0$ , ước lượng tuổi thọ thực tế của nó.

Mô hình

Mô hình sinh của tuổi thọ thực nghiệm là một mô hình phân tầng (hierarchical model) được mô tả như sau

\begin{aligned} Λ & \sim Gamma (α, β) \\ Z | Λ = λ & \sim Exp (λ) \\ X | Z = z & \sim N (z, σ^{2}) \end{aligned}

Dựa trên mô hình này, phân phối kết hợp (joint distribution) của tuổi thọ thực nghiệm và tuổi thọ thực tế là

p_{X, Z, Λ} (x, z, λ) = p_{X | z} (x | z) p_{Z | λ} (z) p_{Λ} (λ) = N (x; z, σ^{2}) Exp (z; λ) Gamma (λ; α, β)

trong đó

\begin{aligned} N (x; z, σ^{2}) & = \frac{1}{\sqrt{2 π σ^{2}}} \exp [- \frac{(x - z)^{2}}{2 σ^{2}}] \\ Exp (z; λ) & = λ \exp (- λ z) 1_{[0, \infty)} (z) \\ Gamma (λ; α, β) & = \frac{β^{α}}{Γ (α)} λ^{α - 1} \exp (- β λ) 1_{(0, \infty)} (λ) \end{aligned}

Tương tự như ví dụ trước, ta mong muốn tính được phân phối hậu nghiệm kết hợp (joint posterior distribution) của cả tham số tỉ lệ $Λ$ lẫn tuổi thọ thực tế $Z$ dựa trên dữ liệu quan sát được:

p_{Z, Λ | x} (z, λ) = \frac{p_{X, Z, Λ} (x, z, λ)}{p_{X} (x)}

Đặt vấn đề

Việc phải ước lượng cả hai tham số $Λ$ và $Z$ cùng lúc khiến cho bài toán trở nên khó khăn hơn. Không những thế, hai biến ngẫu nhiên này còn không độc lập! Điều này đồng nghĩa với việc ta không thể phân rã phân phối hậu nghiệm kết hợp thành các thành phần không phụ thuộc vào nhau, và việc tối ưu cũng không còn dễ dàng nữa.

Suy diễn biến phân với trường trung bình

Suy diễn biến phân với trường trung bình (mean-field variational inference) là một trong những cách tiếp cận phổ biến trong các trường hợp như vậy.

Họ phân phối trường trung bình

Ý tưởng chính của phương pháp này là giả định họ phân phối biến phân (variational distribution family) là họ trường trung bình (mean-field family). Đây là họ các phân bố phân rã được (factorizable). Nói cách khác, các biến ngẫu nhiên được giả định là độc lập với nhau, mỗi biến tuân theo một phân phối phần tử (factor) riêng. Cụ thể, các phân phối trong thuộc họ trường trung bình trên không gian vector ngẫu nhiên $m$ -chiều $\vec{Θ} = (Θ_{1}, \dots, Θ_{m})$ có dạng

q_{\vec{Θ}} (\vec{θ}) = \prod_{i = 1}^{m} q_{Θ_{i}} (θ_{i})

Hãy nhớ lại về định nghĩa ELBO:

ELBO (q_{\vec{Θ}}, p_{X, \vec{Θ}}) = E_{q_{\vec{Θ}}} [\log p_{X, \vec{Θ}} (x, \vec{Θ}) - \log q_{\vec{Θ}} (\vec{Θ})]

Trong đó, phân phối kết hợp có thể được phân tách như sau

p_{X, \vec{Θ}} (x, \vec{Θ}) = p_{X} (x) p_{\vec{Θ} | x} (\vec{Θ})

và do đó

\begin{aligned} E_{q_{\vec{Θ}}} [\log p_{X, \vec{Θ}} (x, \vec{Θ})] & = \log p_{X} (x) + E_{q_{\vec{Θ}}} [\log p_{\vec{Θ} | x} (\vec{Θ})] \end{aligned}

Với việc sử dụng họ phân phối trường trung bình,

\begin{aligned} E_{q_{\vec{Θ}}} [\log q_{\vec{Θ}} (\vec{Θ})] & = \int q_{\vec{Θ}} (\vec{θ}) \log q_{\vec{Θ}} (\vec{θ}) d \vec{θ} \\ = \int (\prod_{j = 1}^{m} q_{Θ_{j}} (θ_{j})) (\sum_{j = 1}^{m} \log q_{Θ_{j}} (θ_{j})) d \vec{θ} \\ = \sum_{j = 1}^{m} \int (\prod_{k = 1}^{m} q_{Θ_{k}} (θ_{k})) \log q_{Θ_{j}} (θ_{j}) d \vec{θ} \\ = \sum_{j = 1}^{m} (\int \prod_{k \neq j} q_{Θ_{k}} (θ_{k}) d {\vec{θ}}_{- j}) (\int q_{Θ_{j}} (θ_{j}) \log q_{Θ_{j}} (θ_{j}) d θ_{j}) \\ = \sum_{j = 1}^{m} E_{q_{Θ_{j}}} [\log q_{Θ_{j}} (Θ_{j})] \end{aligned}

Kết hợp 2 điều trên, ELBO trở thành

\begin{array}{r} ELBO (q_{\vec{Θ}}, p_{X, \vec{Θ}}) = \log p_{X} (x) + E_{q_{\vec{Θ}}} [\log p_{\vec{Θ} | x} (\vec{Θ})] - \sum_{j = 1}^{m} E_{q_{Θ_{j}}} [\log q_{Θ_{j}} (Θ_{j})] \end{array}

Một điểm đáng lưu ý ở đây là với việc sử dụng họ phân phối trường trung bình, tại một tọa độ $k \in [m]$ bất kì,

\begin{aligned} E_{q_{\vec{Θ}}} [\log p_{\vec{Θ} | x} (\vec{Θ})] & = E_{q_{Θ_{k}}} E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- j}} (Θ_{k}) + \log p_{{\vec{Θ}}_{- k} | x} ({\vec{Θ}}_{- k})] \\ = E_{q_{Θ_{k}}} E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (Θ_{k})] + E_{q_{{\vec{Θ}}_{- k}}} [\log p_{{\vec{Θ}}_{- k} | x} ({\vec{Θ}}_{- k})] \end{aligned}

Như vậy, tại một tọa độ bất kì, ELBO có thể được tách được thành

\begin{aligned} ELBO (q_{Θ}, p_{X, \vec{Θ}}) & = \log p_{X} (x) \\ + E_{q_{Θ_{k}}} E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (Θ_{k})] + E_{q_{{\vec{Θ}}_{- k}}} [\log p_{{\vec{Θ}}_{- k} | x} ({\vec{Θ}}_{- k})] \\ - E_{q_{Θ_{k}}} [\log q_{Θ_{k}} (Θ_{k})] - \sum_{j \neq k} E_{q_{Θ_{j}}} [\log q_{Θ_{j}} (Θ_{j})] \end{aligned}

Gia tăng từng tọa độ

Việc có thể viết được ELBO dưới dạng trên cho thấy rằng việc tìm $q_{\vec{Θ}}$ tối ưu ELBO có thể được thực hiện thông qua việc tối ưu tuần tự từng thành phần khi giữ cố định các thành phần còn lại. Phương pháp này được gọi là phương pháp gia tăng từng tọa độ (coordinate ascent).

Xét tọa độ $k \in [m]$ . Khi giữ cố định $q_{Θ_{- k}}$ , ta có thể tìm $q_{Θ_{k}}$ để tối ưu ELBO bằng cách xét đạo hàm bằng không.

ELBO khi chỉ xét các thành phần liên quan đến $q_{Θ_{k}}$ là

ELBO (q_{Θ_{k}}; q_{\vec{Θ}}, p_{X, \vec{Θ}}) \overset{c}{=} E_{q_{Θ_{k}}} {E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (Θ_{k})] - \log q_{Θ_{k}} (Θ_{k})}

Đạo hàm theo $q_{Θ_{k}} (θ_{k})$ là

\frac{\partial}{\partial q_{Θ_{k}} (θ_{k})} ELBO (q_{Θ_{k}}; q_{\vec{Θ}}, p_{X, \vec{Θ}}) = E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (θ_{k})] - \log q_{Θ_{k}} (θ_{k}) - 1

Đặt đạo hàm bằng không, ta có

\log q_{Θ_{k}}^{⋆} (θ_{k}) \overset{c}{=} E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (θ_{k})]

hay

q_{Θ_{k}}^{⋆} (θ_{k}) \propto \exp {E_{q_{{\vec{Θ}}_{- k}}} [\log p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (θ_{k})]}

$p_{Θ_{k} | x, {\vec{Θ}}_{- k}}$ được gọi là phân phối điều kiện hoàn toàn (complete conditional) của biến $Θ_{k}$ . Trong một số trường hợp, ta có thể giả định thêm về nó để thuận tiện cho việc tính toán. Tuy nhiên, việc này nằm quá giới hạn của bài viết này và sẽ được thảo luận ở phần sau (nếu có).

Ở đây, ta tiếp cận nó theo một hướng khác, sử dụng việc

p_{Θ_{k} | x, {\vec{Θ}}_{- k}} (θ_{k}) = \frac{p_{X, \vec{Θ}} (x, θ_{k}, {\vec{Θ}}_{- k})}{p_{X, {\vec{Θ}}_{- k}} (x, {\vec{Θ}}_{- k})}

với $p_{X, {\vec{Θ}}_{- k}} (x, {\vec{Θ}}_{- k})$ không phụ thuộc vào $θ_{k}$ , ta có

q_{Θ_{k}}^{⋆} (θ_{k}) \propto \exp {E_{q_{{\vec{Θ}}_{- k}}} [\log p_{X, \vec{Θ}} (x, θ_{k}, {\vec{Θ}}_{- k})]}

Trong một số trường hợp, ta có thể tính được kỳ vọng này một cách dễ dàng.

Thuật toán CAVI

Tóm lại, thuật toán gia tăng từng tọa độ để tối ưu ELBO khi sử dụng họ phân phối trường trung bình, gọi tên là CAVI (Coordinate Ascent Variational Inference), có thể được mô tả như sau:

Khởi tạo $q_{Θ}$ ngẫu nhiên
Lặp lại cho đến khi hội tụ:
1. Với mỗi $k = 1, \dots, m$ , cập nhật
  $q_{Θ_{k}} (θ_{k}) = \exp {E_{q_{{\vec{Θ}}_{- k}}} [\log p_{X, \vec{Θ}} (x, θ_{k}, {\vec{Θ}}_{- k})]}$
2. Chuẩn hóa $q_{Θ}$

Giải quyết bài toán

Họ phân phối trường trung bình

Trong trường hợp này, ta giả định họ phân phối trường trung bình cho cả hai biến ngẫu nhiên $Λ$ và $Z$ :

q_{Z, Λ} (z, λ) = q_{Z} (z) q_{Λ} (λ)

Một điều đặc biệt là ta không hề giả định gì về từng thành phần $q_{Z}$ và $q_{Λ}$ cả, khác với phần trước khi ta giả định $q_{Z}$ là một phân phối chuẩn.

Tính toán các cập nhật

Theo như thuật toán trên, các cập nhật sẽ là

\begin{aligned} q_{Z} (z) & \propto \exp {E_{q_{Λ}} [\log p_{X, Z, Λ} (x, z, Λ)]} \\ q_{Λ} (λ) & \propto \exp {E_{q_{Z}} [\log p_{X, Z, Λ} (x, Z, λ)]} \end{aligned}

Phân phối kết hợp Để tính được bước cập nhật cho từng tọa độ, ta cần biến đổi phân phối kết hợp $p_{X, Z, Λ}$ .

\begin{aligned} p_{X, Z, Λ} (x, z, λ) & = p_{X | Z, Λ} (x | z, λ) p_{Z | Λ} (z | λ) p_{Λ} (λ) \\ = N (x; z, σ^{2}) Exp (z; λ) Gamma (λ; α, β) \\ \propto \exp [- \frac{1}{2 σ^{2}} (x - z)^{2}] λ \exp (- λ z) λ^{α - 1} \exp (- β λ) \end{aligned}

Từ đó, ta có thể tính được

\begin{aligned} \log p_{X, Z, Λ} (x, z, λ) & \overset{c}{=} - \frac{1}{2 σ^{2}} (x - z)^{2} - λ z + α \log λ - β λ \end{aligned}

Cập nhật của $q_{Z}$ Lược bỏ đi các thành phần không liên quan đến $z$ , ta có

\begin{aligned} E_{q_{Λ}} [\log p_{X, Z, Λ} (x, z, Λ)] & \overset{c}{=} - \frac{1}{2 σ^{2}} (x - z)^{2} - E_{q_{Λ}} [Λ] z \\ \overset{c}{=} - \frac{1}{2 σ^{2}} z^{2} + (\frac{x}{σ^{2}} - E_{q_{Λ}} [Λ]) z \\ = - \frac{1}{2 σ^{2}} [z^{2} - 2 (x - σ^{2} E_{q_{Λ}} [Λ]) z] \\ \overset{c}{=} - \frac{1}{2 σ^{2}} {[z - (x - σ^{2} E_{q_{Λ}} [Λ])]}^{2} \end{aligned}

Như vậy, cập nhật cho $q_{Z}$ sẽ là

\begin{aligned} q_{Z} (z) & \propto \exp {- \frac{1}{2 σ^{2}} {[z - (x - σ^{2} E_{q_{Λ}} [Λ])]}^{2}} \end{aligned}

Nếu để ý, ta sẽ nhận ra rằng

q_{Z} (z) = N (z; x - σ^{2} E_{q_{Λ}} [Λ], σ^{2})

Ở đây, $E_{q_{Λ}} [Λ]$ đang được bỏ ngõ, do ta chưa có đủ thông tin về $q_{Λ}$ . Ta sẽ tính nó ngay sau đây.

Cập nhật của $q_{Λ}$ Lược bỏ đi các thành phần không liên quan đến $λ$ , ta có

\begin{aligned} E_{q_{Z}} [\log p_{X, Z, Λ} (x, Z, λ)] & \overset{c}{=} α \log λ - (β + E_{q_{Z}} [Z]) λ \end{aligned}

Như vậy, cập nhật cho $q_{Λ}$ sẽ là

\begin{aligned} q_{Λ} (λ) & \propto \exp {α \log λ - (β + E_{q_{Z}} [Z]) λ} = λ^{α} \exp (- (β + E_{q_{Z}} [Z]) λ) \end{aligned}

Nếu để ý, ta sẽ lại nhận ra rằng

q_{Λ} (λ) = Gamma (λ; α + 1, β + E_{q_{Z}} [Z])

Lại một lần nữa, $E_{q_{Z}} [Z]$ đang được bỏ ngỏ. Tuy nhiên, dựa vào kết quả ở trên, ta biết $q_{Z}$ là một phân phối chuẩn, và do đó $E_{q_{Z}} [Z]$ có thể được tính toán một cách dễ dàng. Song song đó, ta cũng có thể tính $E_{q_{Λ}} [Λ]$ một cách dễ dàng.

Thuật toán

Tóm lại, thuật toán để giải quyết bài toán đầu đề có thể được mô tả như sau:

Khởi tạo $α_{Λ}^{(0)} = α, β_{Λ}^{(0)} = β$
Với mỗi $t = 0, 1, \dots$ :
1. Cập nhật
  $μ_{Z}^{(t + 1)} = x - σ^{2} \frac{α_{Λ}^{(t)}}{β_{Λ}^{(t)}}$
  và theo đó
  $q_{Z}^{(t + 1)} (z) = N (z; μ_{Z}^{(t + 1)}, σ^{2})$
2. Cập nhật
  $\begin{aligned} α_{Λ}^{(t + 1)} & = α_{Λ}^{(t)} + 1 \\ β_{Λ}^{(t + 1)} & = β_{Λ}^{(t)} + μ_{Z}^{(t + 1)} \end{aligned}$
  và theo đó
  $q_{Λ}^{(t + 1)} (λ) = Gamma (λ; α_{Λ}^{(t + 1)}, β_{Λ}^{(t + 1)})$

Minh họa

Nếu quan tâm, mọi người có thể thử tại đây để so sánh giữa phân phối hậu nghiệm chính xác và phân phối biến phân tại các bước của thuật toán.

Hình trên vẽ mặt cắt $h = 0.005$ của các phân phối. Trong đó, theo thứ tự màu đen, đỏ, xanh lam, xanh lục, cam, tím lần lượt là phân phối hậu nghiệm chính xác và các phân phối biến phân theo CAVI tại các bước từ 1 đến 5.

Có thể thấy giả định trường trung bình đã giới hạn khả năng xấp xỉ của phân phối biến phân. Bên cạnh đó, do ta đang tối ưu phân kỳ KL $KL (q_{Z, Λ} ‖ p_{Z, Λ | x})$ , phân phối biến phân không cố gắng che phủ (như trong trường hợp tối ưu phân kỳ KL hướng ngược lại) mà thay vào đó muốn nằm gọn trong phân phối hậu nghiệm chính xác.

Tổng kết và hứa hẹn

Trong bài viết này, chúng ta đã thảo luận về một ví dụ phức tạp hơn của suy diễn biến phân, và cách tiếp cận nó thông qua phương pháp suy diễn biến phân với trường trung bình. Bài viết này có thể sẽ được cập nhật để thêm phần cài đặt và mô phỏng kết quả.

Trong các phần tiếp theo (nếu có), chúng ta sẽ xem xét một số ví dụ khác, đôi chút phức tạp hơn, nhưng vẫn có thể giải quyết được bằng phương pháp này. Hy vọng được gặp lại!

Suy diễn Biến phân (Phần 2)

Authors

Affiliations

Published

Bối cảnh

Mô hình

Đặt vấn đề

Suy diễn biến phân với trường trung bình

Họ phân phối trường trung bình

Gia tăng từng tọa độ

Thuật toán CAVI

Giải quyết bài toán

Họ phân phối trường trung bình

Tính toán các cập nhật

Thuật toán

Minh họa

Tổng kết và hứa hẹn

Footnotes