ขั้นตอนวิธีการของหวังและลันเดา

ขั้นตอนวิธีของหวังและลันเดา (อังกฤษ: Wang and Landau algorithm) เป็นขั้นตอนวิธีที่ถูกเพิ่มตามขั้นตอนวิธีการแบบ Metropolis Monte Carlo ซึ่งถูกเสนอขึ้นโดยฟุเกา หวัง (Fugao Wang) และ เดวิด พี.ลันเดา (David P. Landau) ซึ่งขั้นตอนวิธีการนี้ถูกออกแบบเพื่อให้สามารถคำนวณความหนาแน่นของสภาวะของสารต่าง ๆ ที่ถูกจำลองโดยระบบคณิตกร (คอมพิวเตอร์) เช่น แบบจำลองไอซิ่งของแก้วที่ถูกหมุน หรือ แบบจำลองอะตอมในสนามของแรงในโมเลกุลต่าง ๆ แนวคิดหลักของขั้นตอนวิธีการนี้คือการสร้างความหนาแน่นของสภาวะของสารโดยใช้ความจริงที่ว่า ระบบนั้นเมื่อถูกจำลองแล้ว จะสร้างความหนาแน่นของสภาวะของสารขึ้น ขั้นตอนวิธีการนี้เป็นวิธีการทั่วไปที่สามารถทำให้รู้ได้ถึงความหนาแน่นของสภาวะของสาร ซึ่งเป็นที่ยอมรับโดยทั่วไป วิธีการแบบ Metropolis Monte Carlo ถูกตั้งขึ้นตามความคิดของโบลท์แมนที่ว่า เมื่อเพิ่มอุณหภูมิให้กับโมเลกุลที่กระจัดกระจายในระหว่างช่วงพลังงานที่สูง หรือเรียกอีกอย่างว่าในสถานะที่ไม่น่าพึงพอใจ กับ ช่วงพลังงานต่ำ หรือเรียกอีกอย่างว่าในสถานะที่น่าพึงพอใจ ซึ่งมีความน่าจะเป็นที่พลังงานนั้นจะมีผลต่อความหนาแน่นของสภาวะของสาร ดังนั้นถ้าหากว่ามีสารใดสารหนึ่งอยู่ในภาวะที่พลังงานต่ำ แล้วได้รับพลังงานมากเพียงพอที่จะทำให้สถานะของสารเปลี่ยนแปลงไป ซึ่งทำให้สภาวะความหนาแน่นของสารที่มีพลังงานมากถูกทำลายไป ยกตัวอย่างเข่นการที่ขี้ผึ้งละลาย ในสภาวะอุณหภูมิต่ำ เฉพาะสถานะของแข็งเท่านั้นที่สามารถจับต้องได้ เมื่ออุณหภูมิเพิ่มขึ้น สถานนะที่ไม่น่าพึงพอใจจะเพิ่มมากขึ้น จนในที่สุดทำให้ของแข็งนั้นหลอมเหลวกลายเป็นของเหลว ในทฤษฎีของขั้นตอนวิธีการของหวังและลันเดาสามารถนำไปใช้กับระบบใด ๆ ที่ถูกกำหนดสถานะโดยพลังงาน (หรือค่าอื่น ๆ) ในปัจจุบัน ขั้นตอนวิธีการนี้ถูกนำไปใช้กับการแก้ปัญหาปริพันธ์ต่าง ๆ และการหดตัวของโปรตีน

ขั้นตอนวิธี

ขั้นแรก ถ้าเกิดว่าการจัดสรรแบบพลวัตไม่ถูกนำมาใช้ จำนวนน้อยสุดกับมากสุดที่เป็นไปได้ของสภาวะของระบบจะต้องถูกคำนวณก่อน ในขณะเดียวกัน ถ้าหากเราพิจารณาแบบจำลองไอซิ่งในขณะที่อยู่ในสถานะที่น่าพอใจเป็นอย่างยิ่งและอยู่รอบจุดหมุนที่มีทิศทางเดียวกัน และเมื่ออยู่ในสถานะที่น่าพึงพอใจน้อยลงมา จะสมมติการใช้แบบแผนแบบตารางหมากรุก ช่วงของค่าจะถูกนำมาหารเป็นค่าที่แยกจากกันโดยสิ้นเชิงในกราฟแสดงความถี่เชิงสถิติ (Discrete Histogram) เราจะทำการกำหนดค่ามากสุดและน้อยสุดของกระบวนการคำนวณนี้เป็น Emax และ Emin ซึ่งจำเป็นต้องใช้ค่าความต่างระหว่างพลังงานมากสุดและน้อยสุดที่แม่นยำ (Δ) จะได้สมการว่า

$N={\frac {E_{\rm {max}}-E_{\rm {min}}}{\Delta }}$

หลังจากนั้นอาเรย์ซึ่งมีจำนวน N ช่องจะต้องเก็บค่าความหนาแน่นของสภาวะของสาร ในขั้นตอนนี้จะต้องใส่ใจเรื่องความถูกต้องของ Δ เพราะว่ายิ่งค่าNเพิ่มขึ้นเวลาที่ต้องใช้ในการจำลองรูปแบบก็จะเพิ่มอย่างรวดเร็ว ในขั้นต้นค่าความหนาแน่นของสภาวะของสาร g (E) นั้นยังไม่สามารถทราบได้ ดังนั้นค่าทุกช่องในอาเรย์จะต้องถูกตั้งเป็นค่า ๆ หนึ่งซึ่งเป็นค่าปริยาย เนื่องจากค่าความหนาแน่นโดยทั่วไปนั้นจะมีค่าระหว่างพหุคูณของมิติ ดังนั้นเป็นการเก็บข้อมูลควรเก็บเป็นค่าลอการิททึ่มของค่าความหนาแน่นของสภาวะของสาร คือ ln [g (E)] ยิ่งไปกว่านั้น ค่าของกราฟแสดงความถี่เชิงสถิติ H (E) จะต้องถูกทำให้คงที่ไว้ ค่าเริ่มต้นของจุดที่ยังไม่ได้ถูกพบของ ln[g (E)] และ H (E) จะต้องถูกตั้งเป็น0 ในจุดแต่ละจุดถูกแต้มสีเหมือนกับการจำลองแบบMonte-Carlo ซึ่งอาจจะถูกยอมรับโดยพิจารณาจากค่าความหนาแน่นของสภาวะของสารแทนที่จะเป็นการพิจารณาจากเกณฑ์ของMetropolis ซึ่งการเลือกจุดจะถูกยอมรับถ้า

$p<\min \left\{1,{\frac {g(E)}{g(E')}}\right\}$

โดยpเป็นค่า ๆ หนึ่งซึ่งมากกว่าเท่ากับ0แต่น้อยกว่า1 Eคือค่าพลังงานของสถานะในปัจจุบัน และE' เป็นค่าพลังงานที่จะต้องเปลี่ยนแปลงไป หลังจากที่การเคลื่อนย้ายของขุดถูกยอมรับหรือถูกปฏิเสธ ค่าการถูกพบในกราฟแสดงความถี่เชิงสถิติ H (E) จะถูกเพิ่มขึ้น1 และค่าความหนาแน่นของสภาวะของสารในกราฟแสดงความถี่เชิงสถิติ H (E) จะถูกคูณโดยค่าคงที่ค่าหนึ่ง

$g(E)\rightarrow g(E)\times f$

ซึ่งค่าคงที่มักถูกแทนที่ด้วย $f=c\simeq 2.72$ (ตัวเลขของออยเลอร์) เมื่อกราฟแสดงความถี่เชิงสถิติ H (E) ถูกค้นพบจนทั่วแล้ว g (E) จะมีค่าประมาณที่แม่นยำของค่าความหนาแน่นของสภาวะของสาร ซึ่งอยู่ในขอบเขตของการเปลี่ยนแปลงของค่าคงที่ ในขั้นตอนนี้กราฟแสดงความถี่เชิงสถิติ H (E) จะถูกตั้งค่าเป็น0และค่าการเปลี่ยนแปลงของค่าคงที่จะถูกลดดังนี้

$f_{i+1}\rightarrow {\sqrt {f_{i}}}$

การทำขั้นตอนข้างบนซ้ำ ๆ จะทำให้เกิดปัญหาที่ว่า g (E) จะต้องเก็บข้อมูลจำนวนมหาศาลตามที่ได้กล่าวไว้ในเบื้องต้น ถ้าต้องการหลีกเลี่ยงปัญหานี้ ค่าของ

$\ln[g(E)]\rightarrow \ln[g(E)]+\ln[f]$

จะต้องถูกนำมาใช้ ซึ่งขั้นตอนนี้ค่า f = e ที่ได้ตั้งไว้ก่อนหน้านี้จะมีประโยชน์เพราะ ln (e) = 1 และค่าคงที่จะถูกเปลี่ยนไปดังนี้

$\ln(f_{i+1})\rightarrow (1/2)\ln(f_{i})$

ขั้นตอนวิธีการนี้สามารถใช้งานได้เพราะการเปลี่ยนแปลงค่าในขั้นตอนต่าง ๆ นั้นเป็นเหมือนการจำลองค่าความหนาแน่นของสภาวะของสารจริง ๆ ดังนั้นถ้าค่าความหนาแน่นของสภาวะของสารถูกสร้างโดยกราฟแสดงความถี่เชิงสถิติถูกยอมรับได้เมื่อเป็นค่าจำนวนเลขที่กลับกัน ค่าที่ได้จะเป็นค่าประมาณที่มีความแม่นยำ ซึ่งในความเป็นจริงแล้วจริงแล้ว H (E) จะไม่สามารถไปสู่จุดที่ราบอย่างสมบูรณ์แบบได้ ดังนั้นบางเกณฑ์อาจกล่าวว่า ค่าที่บางที่สุดที่เป็นไปได้จะอยู่ในช่วงระหว่าง10%ของจุดสูงสุดและจุดต่ำสุด การเปลี่ยนแปลงของค่าคงที่โดยใช้รากที่สองอาจทำให้เกิดความผิดพลาดเพราะการอิ่มตัวได้ ดังนั้นวิธีที่จะสามารถหลีกเลี่ยงปัญหานี้ได้คือใช้ค่าคงที่เป็น1 / t โดยที่ tคือเวลาในการจำลอง ในแบบจำลองง่าย ๆ เช่นแบบจำลองไอซิ่ง ปัญหานี้อาจไม่ได้มีปัญหามากนัก แต่หากเป็นระบบที่ซับซ้อนเช่นโปรตีน เยื่อหุ้มเซล และการหาปริพันธ์หลายมิติจะเป็นปัญหาใหญ่มาก เพราะต้องใช้เวลาคำนวณทีเยอะมาก การหาค่าความหนาแน่นของสภาวะของสารนั้นไม่ขึ้นต่ออุณหภูมิ แต่ว่าจะสามารถพิจารณาผลลัพธ์เมื่อมีค่าของอุณหภูมิมาเกี่ยวข้องได้

การทดสอบระบบ

เราต้องการหาค่าพลังงานศักย์ของการแกว่งแบบฮาร์โมนิคซึ่งมีสูตรดังนี้ E (x) = x2 จากการพิจารณาค่าความหนาแน่นของสาร ได้ว่า

$g(E)=\int \delta (E(x)-E_{0})dx=\int \delta (x^{2}-E_{0})dx$

เมื่อทำการแก้ปริพันธ์แล้วจะได้

$g(E)\propto E^{-1/2}$

โดยทั่วไปนั้นการหาค่าความหนาแน่นของสภาวะของสารสำหรับการแกว่งแบบฮาร์โมนิคหลายมิติ จะถูกกำหนดโดยพลังงานE ซึ่งเลขชี้กำลังจะถูกกำหนดตามมิติของระบบ

รหัสเทียม

รหัสเทียมนี้เป็นรหัสแสดงขั้นตอนวิธีการของหวังและลันเดาในภาษาการเขียนโปรแกรม D

real[real] log_g; // อาเรย์นี้เก็บค่าล็อกการิททึ่มของค่าความหนาแน่นของสภาวะของสาร
uint[real] H; // อาเรย์เก็บค่ากราฟแสดงความถี่ทางสถิติ

real E_old = the_system.energy (); // เก็บค่าพลังงานเก่า
real F = 1; //เก็บค่าคงที่

while (F > epsilon) {
the_system.evolve (); // สร้างการพัฒนาการของระบบ
real E_new = the_system.energy (); // คำนวณค่าพลังงานของระบบปัจจุบัน

if (log_g[E_new] < log_g[E_old]) // ตรวจสอบว่าข้อเสนอสามารถถูกยอมรับได้หรือไม่
E_old = E_new;
else if (random () < exp (log_g[E_old] - log_g[E_new]) )
E_old = E_new;
else
the_system.reject_and_undo ();

H[E_old]++; // เปลี่ยนแปลงค่ากราฟแสดงความถี่ทางสถิติและค่าความหนาแน่นของสภาวะของสาร
log_g[E_old] += F;

if (is_flat (H) ) { // ถ้ากราฟแสดงความถี่มีลักษณะแบนราบ
H[] = 0; // ให้ตั้งค่ากราฟแสดงความถี่เป็น0และลดค่าคงที่
F *= 0.5;
}
}

return log_g;

รหัสเทียมเป็นตัวอักษร

เริ่มต้นจาก ตั้งค่า f=e

do{ do{

Metropolis แก้ไขค่าซึ่งขึ้นกับความน่าจะเป็น W[i->j] = min[1, r (Ei) / r (Ej)]

ปรับตัว r (E) ในแต่ละขั้น : r (E) <- r (E) x f

เพิ่มค่ากราฟแสดงความถี่เชิงสถิติ H (E) ++ }

จนกว่า H (E) มีลักษณะแบนราบ (ประมาณ ร้อยละ20) }

ลด f } จนกว่า f-1 มีค่าประมาณ 10^ (-8)

ปัญหาที่เกี่ยวข้อง

ความหนาแน่นในสภาวะของสาร Density of states (DOS) ของระบบนั้นสามารถเขียนแทนด้วย จำนวนขั้น ต่อ ช่วงระดับพลังงานแต่ละค่าที่สร้างสามารถมีอยู่ได้ ซึ่งแตกต่างกับ isolated systems (เช่น อะตอม หรือ โมเลกุลแก๊ส) การกระจายตัวของ ความหนาแน่น นั้นมีความต่อเนื่อง DOS ที่สูงนั้นหมายถึง มีหลาย state ที่เป็นไปได้ โดยทั่วไปนั้น DOS จะเป็นค่าเฉลี่ยบนที่ว่างและเวลาของระบบ การคำนวณ DOS

ปกติระบบที่เราสนใจนั้นจะมีความซับซ้อน เช่น สารประกอบ, biomolecules, โพลิเมอร์ ฯลฯ ซึ่งส่วนใหญ่ระบบที่มีความซับซ้อน การคำนวณ DOS เทียบจะเป็นไปได้ยากหรือเป็นไปไม่ได้เลย จึงมีการใช้คอมพิวเตอร์ simulate algorithms เพื่อประเมิน DOS ให้ได้ถูกต้องใกล้เคียง หนึ่งใน algorithms ที่มีการใช้ คือ Wang and Landau algorithm

ใน Wang and Landau algorithm นั้น ข้อมูล DOS ก่อนหน้าจำเป็นต้องรู้ cost function ซึ่งเทียบได้กับ พลังงานของระบบนั้นถูก discretized ซึ่งแต่ละรอบ bin i จะ update DOS ของ state g (i) ใน histogram ด้วย g (i) -> g (i) +f โดย f เรียกว่า modification factor

เมื่อแต่ละ bin ใน histogram ถูก update เป็นจำนวน 10-15 รอบ ค่า f จะถูกลดทอนลงด้วย f (n+1) =1/2f (n) โดยที่ n คือรอบที่ n ของการ update

การ simulate นั้นจะทำไปเรื่อย ๆ จนกว่า modification factor จะน้อยกว่าค่า threshold ที่ได้กำหนดไว้ เช่น fn < 10 − 8.

ข้อดีของ ขั้นตอนวิธีการของหวังและลันเดา ที่เหนือกว่าขั้นตอนวิธีอื่น ๆ (เช่น multicanonical simulations and Parallel tempering เป็นต้น) นั่นคือ การจำลองนั้นให้ผลลัพธ์หลักคือ DOS นอกจากนี้ การจำลองเป็นอิสระต่ออุณหภูมิ การคำนวณ DOS ของระบบโดยใช้ขั้นตอนวิธีนี้มักใช้ในการคำนวณ DOS ของโปรตีน

ตัวอย่างการใช้งาน

การติดตามการเคลื่อนไหวอย่างทันทีทันใด

ขั้นตอนวิธีการของหวังและลันเดาสามารถนำไปใช้กับการติดตามการเคลื่อนไหวอย่างทันทีทันใด

เราสามารถตั้งสมมติฐานขั้นตอนวิธีที่ใช้ติดตามการเคลื่อนไหวในทันทีทันใดโดยใช้พื้นฐานจากขั้นตอนวิธีของหวังและลันเดาเป็นตัวอย่าง ซึ่งมีประสิทธิภาพในการรองรับภาพที่เคลื่อนไหวในทันทีทันใด การเคลื่อนไหวในทันทีทันใดมักจะทำให้การติดตามการเคลื่อนไหวธรรมดาล้มเหลว เพราะว่าความราบรื่นของกิริยาท่าทางมีน้อยเกินไป เราจึงสามารถนำขั้นตอนวิธีการของหวังและลันเดามาใช้เพราะมีความสามารถในการตรวจสอบสถานะทางกายภาพ และนำขั้นตอนวิธีนี้มาใช้ร่วมกับขั้นตอนวิธีแบบ Markov Chain Monte Carlo ซึ่งเป็นพื้นฐานในการติดตามการเคลื่อนไหวปกติ

วิธีการนี้ทำให้ความไม่ราบรื่นของการเคลื่อนไหวลดลง โดยใช้ความคล้ายคลึงกับลักษณะความหนาแน่นของสภาวะของสาร ซึ่งสามารถหาค่าโดยประมาณได้จากขั้นตอนวิธีการของหวังและลันเดา ซึ่งการใช้ขั้นตอนวิธีนี้ได้ผล ทำให้ความแม่นยำในการติดตามการเคลื่อนไหวราบรื่นขึ้น เพราะว่าลักษณะความหนาแน่นของสภาวะของสารสามารถจับการเคลื่อนไหวอย่างทันทีทันใดได้อย่างแน่นหนา

ผลการทดลองพบว่าการใช้ขั้นตอนวิธีนี้ได้ผลเป็นอย่างดี ตัวอย่างที่ใช้ทดสอบพบว่าลักษณะความหนาแน่นของสภาวะของสารสามารถควบคุมการติดตามการเคลื่อนไหวโดยที่ไม่เสียเวลาเลย สรุปได้ว่า การใช้ขั้นตอนวิธีของหวังและลันเดาในการติดตามภาพที่เคลื่อนไหวในทันทีทันใดสามารถติดตามได้อย่างแม่นยำและถูกต้องเป็นอย่างมาก

การใช้งานกับHeisenberg Chain

L=10 ให้ค่า Heisenberg chain L = 250

พลังงานเสรีกับค่าความยุ่งเหยิงสามารถนำมาใช้ได้

การจำลองหนึ่งครั้งต่อทุกอุณหภูมิ

T>0.05จูล

ได้ผลดังนี้

อื่น ๆ

การเพิ่มสมรรถภาพของขั้นตอนวิธีของหวังและลันเดาโดยใช้ขอบเขตที่สามารถเปลี่ยนแปลงได้

ขั้นตอนวิธีของหวังและลันเดาในการจำลองระบบขนาดใหญ่ จำเป็นต้องทำการหารช่วงของพลังงานลงเป็น ขอบเขต หลาย ๆ ขอบเขต แล้วนำมารวมกันเพื่อหาคำตอบ ผลของความหนาแน่นของสถานะของสารที่มีระบบขนาดใหญ่มากเกินไป ทำให้เกิดปัญหาการจำลองพลังงานแลททิซของโพลีเมอร์และแบบจำลอง 5ขั้นตอนของ Potts เพราะการจำลองของทั้งสองอย่างนั้นต้องใช้การจำลองที่มีขนาดใหญ่มาก

วิธีการแก้ปัญหาที่เกิดจากระบบมีขนาดใหญ่มากเกินไปทำได้โดยใช้ ขอบเขตที่สามารถเปลี่ยนแปลงได้แทนที่จะใช้ขอบเขตที่มีขนาดคงที่ ค่าของขอบเขตนั้นขึ้นกับค่าของพลังงานในกลุ่มในกราฟแสดงความถี่เชิงสถิติ (Discrete Histogram) มีลักษณะแบนราบหรือไม่เป็นลำดับขั้นของการจำลอง

การเปลี่ยนแปลงขนาดของขอบเขตในแต่ละครั้งของการเปลี่ยนค่าf เราจะทำการตัดขอบของค่า ที่มีผลต่อการจำลองซึ่งใช้ในการจำลองโครงสร้างของขอบเขตปกติ ขอบเขตแบบเปลี่ยนแปลงได้สามารถเพิ่มช่วงของระบบที่จะนำมาจำลองสภานะสภาวะของสารได้อย่างมีประสิทธิภาพ

ความคิดเห็นของผู้เขียน

ขั้นตอนวิธีนี้เป็นขั้นตอนวิธีที่สามารถใช้ประโยชน์ได้มาก นอกจากที่จะใช้ในเชิงฟิสิกส์แล้ว ยังสามารถใช้ในเชิงพาณิชย์ได้ เช่นการทำซอฟต์แวร์ การติดตามการเคลื่อนที่แบบทันทีทันใด สามารถใช้ในการตรวจสอบการเคลื่อนไหวของนักเทนนิส อาจดัดแปลงไปเป็นการทำงานอื่น ๆ ที่มีประโยชน์ต่อมวลมนุษยชาติก็เป็นได้

อ้างอิง

• Wang, Fugao and Landau, D. P. (Mar 2001). "Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States". Phys. Rev. Lett. (American Physical Society) 86 (10) : 2050–2053. Bibcode 2001PhRvL..86.2050W. doi:10.1103/PhysRevLett.86.2050. PMID 11289852.

• R. E. Belardinelli and S. Manzi and V. D. Pereyra (Dec 2008). "Analysis of the convergence of the 1∕t and Wang-Landau algorithms in the calculation of multidimensional integrals". Phys. Rev. E (American Physical Society) 78: 067701. doi:10.1103/PhysRevE.78.067701.

• P. Ojeda and M. Garcia and A. Londono and N.Y. Chen (Feb 2009). "Monte Carlo Simulations of Proteins in Cages: Influence of Confinement on the Stability of Intermediate States". Biophys. Jour. (Biophysical Society) 96 (3) : 1076–1082. Bibcode 2009BpJ....96.1076O. doi:10.1529/biophysj.107.125369.

• Belardinelli, R. E. and Pereyra, V. D. (2007). "Wang-Landau algorithm: A theoretical analysis of the saturation of the error". Jour. Chem. Phys. 127 (18) : 184105. doi:10.1063/1.2803061.

http://cv.snu.ac.kr/research/~wlmctracker/index.html

http://www.mendeley.com/research/improving-wanglandau-sampling-with-adaptive-windows/

http://www.comp-phys.org/lugano04/Talks/QWL.pdf เก็บถาวร 2009-01-08 ที่ เวย์แบ็กแมชชีน