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

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β>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(0,). 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

ΛGamma(α,β)Z|Λ=λExp(λ)X|Z=zN(z,σ2)

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à

pX,Z,Λ(x,z,λ)=pX|z(x|z)pZ|λ(z)pΛ(λ)=N(x;z,σ2)Exp(z;λ)Gamma(λ;α,β)

trong đó

N(x;z,σ2)=12πσ2exp[(xz)22σ2]Exp(z;λ)=λexp(λz)1[0,)(z)Gamma(λ;α,β)=βαΓ(α)λα1exp(βλ)1(0,)(λ)

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:

pZ,Λ|x(z,λ)=pX,Z,Λ(x,z,λ)pX(x)

Đặt vấn đề

Việc phải ước lượng cả hai tham số Λ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 Θ=(Θ1,,Θm) có dạng

qΘ(θ)=i=1mqΘi(θi)

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

ELBO(qΘ,pX,Θ)=EqΘ[logpX,Θ(x,Θ)logqΘ(Θ)]

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

pX,Θ(x,Θ)=pX(x)pΘ|x(Θ)

và do đó

EqΘ[logpX,Θ(x,Θ)]=logpX(x)+EqΘ[logpΘ|x(Θ)]

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

EqΘ[logqΘ(Θ)]=qΘ(θ)logqΘ(θ)dθ=(j=1mqΘj(θj))(j=1mlogqΘj(θj))dθ=j=1m(k=1mqΘk(θk))logqΘj(θj)dθ=j=1m(kjqΘk(θk)dθj)(qΘj(θj)logqΘj(θj)dθj)=j=1mEqΘj[logqΘj(Θj)]

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

ELBO(qΘ,pX,Θ)=logpX(x)+EqΘ[logpΘ|x(Θ)]j=1mEqΘj[logqΘj(Θj)]

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[m] bất kì,

EqΘ[logpΘ|x(Θ)]=EqΘkEqΘk[logpΘk|x,Θj(Θk)+logpΘk|x(Θk)]=EqΘkEqΘk[logpΘk|x,Θk(Θk)]+EqΘk[logpΘk|x(Θk)]

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

ELBO(qΘ,pX,Θ)=logpX(x)+EqΘkEqΘk[logpΘk|x,Θk(Θk)]+EqΘk[logpΘk|x(Θk)]EqΘk[logqΘk(Θk)]jkEqΘj[logqΘj(Θj)]

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Θ 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[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

ELBO(qΘk;qΘ,pX,Θ)=cEqΘk{EqΘk[logpΘk|x,Θk(Θk)]logqΘk(Θk)}

Đạo hàm theo qΘk(θk)

qΘk(θk)ELBO(qΘk;qΘ,pX,Θ)=EqΘk[logpΘk|x,Θk(θk)]logqΘk(θk)1

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

logqΘk(θk)=cEqΘk[logpΘk|x,Θk(θk)]

hay

qΘk(θk)exp{EqΘk[logpΘk|x,Θk(θk)]}

pΘk|x,Θ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,Θk(θk)=pX,Θ(x,θk,Θk)pX,Θk(x,Θk)

với pX,Θk(x,Θk) không phụ thuộc vào θk, ta có

qΘk(θk)exp{EqΘk[logpX,Θ(x,θk,Θ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:

  1. Khởi tạo qΘ ngẫu nhiên
  2. Lặp lại cho đến khi hội tụ:

    1. Với mỗi k=1,,m, cập nhật

      qΘk(θk)=exp{EqΘk[logpX,Θ(x,θk,Θ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 ΛZ:

qZ,Λ(z,λ)=qZ(z)qΛ(λ)

Một điều đặc biệt là ta không hề giả định gì về từng thành phần qZqΛ cả, khác với phần trước khi ta giả định qZ 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à

qZ(z)exp{EqΛ[logpX,Z,Λ(x,z,Λ)]}qΛ(λ)exp{EqZ[logpX,Z,Λ(x,Z,λ)]}

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 pX,Z,Λ.

pX,Z,Λ(x,z,λ)=pX|Z,Λ(x|z,λ)pZ|Λ(z|λ)pΛ(λ)=N(x;z,σ2)Exp(z;λ)Gamma(λ;α,β)exp[12σ2(xz)2]λexp(λz)λα1exp(βλ)

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

logpX,Z,Λ(x,z,λ)=c12σ2(xz)2λz+αlogλβλ

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

EqΛ[logpX,Z,Λ(x,z,Λ)]=c12σ2(xz)2EqΛ[Λ]z=c12σ2z2+(xσ2EqΛ[Λ])z=12σ2[z22(xσ2EqΛ[Λ])z]=c12σ2[z(xσ2EqΛ[Λ])]2

Như vậy, cập nhật cho qZ sẽ là

qZ(z)exp{12σ2[z(xσ2EqΛ[Λ])]2}

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

qZ(z)=N(z;xσ2EqΛ[Λ],σ2)

Ở đây, EqΛ[Λ] đ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ó

EqZ[logpX,Z,Λ(x,Z,λ)]=cαlogλ(β+EqZ[Z])λ

Như vậy, cập nhật cho qΛ sẽ là

qΛ(λ)exp{αlogλ(β+EqZ[Z])λ}=λαexp((β+EqZ[Z])λ)

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

qΛ(λ)=Gamma(λ;α+1,β+EqZ[Z])

Lại một lần nữa, EqZ[Z] đang được bỏ ngỏ. Tuy nhiên, dựa vào kết quả ở trên, ta biết qZ là một phân phối chuẩn, và do đó EqZ[Z] có thể được tính toán một cách dễ dàng. Song song đó, ta cũng có thể tính EqΛ[Λ] 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:

  1. Khởi tạo αΛ(0)=α,βΛ(0)=β
  2. Với mỗi t=0,1,:

    1. Cập nhật

      μZ(t+1)=xσ2αΛ(t)βΛ(t)

      và theo đó

      qZ(t+1)(z)=N(z;μZ(t+1),σ2)
    2. Cập nhật

      αΛ(t+1)=αΛ(t)+1βΛ(t+1)=βΛ(t)+μZ(t+1)

      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(qZ,ΛpZ,Λ|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!

Footnotes