ก้าวแรก HPC #05: อัลกอริทึมอย่างง่ายและการปรับปรุงประสิทธิภาพ

เกริ่นนำ

บทความนี้เป็นบทความแรกตั้งแต่เริ่มทำเว็ปนี้ที่เขียนเป็นโปรแกรมจริงๆ  โปรแกรมที่เราจะเขียนกันนั้นก็มาจากทฤษฎีบทก่อนหน้าเรื่องการหาค่า Bernoulli Number นั่นเอง ซึ่งเป็นแบบง่ายที่สุด เป้าหมายที่แท้จริงของบทความนี้คงไม่ใช่อยู่ตัวอัลกอริทึมแต่อย่างใด เพราะไม่มีอะไรซับซ้อนมากมายตรงไปตรงมา แต่หากอยู่ที่อยู่มุมมองเรื่องของประสิทธิภาพ การเขียนโปรแกรมให้ทำงานได้ถูกต้องนั้นเป็นแค่เพียงขั้นแรกเท่านั้น ขั้นต่อไปเราต้องปรับปรุงประสิทธิภาพให้ดีด้วย ในบทความนี้จะใช้ Python เป็นตัวทำ Prototype เมื่อเรียบร้อยแล้วจึงไปใช้ C++ ในกระบวนการสุดท้าย ท่านจะได้เห็น workflow การทำงานจริง เชิญสดับ

อัลกอริทึมในการหา Bernoulli Number อย่างง่าย

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

เมื่อเข้าไปใน Wikipedia เรื่องของ Bernoulli Number  ที่ http://en.wikipedia.org/wiki/Bernoulli_number#Recursive_definition จะพบนิยาม Recursive ซึ่งผมขอปรับปรุงให้อ่านง่ายขึ้นดังนี้ครับ

B_{0}=1
B_{m}=1-\sum_{k=0}^{m-1}\binom{m}{k}\frac{B_{k}}{m-k+1}

วิธีนี้เป็นการนิยาม Bernoulli Number โดยวิธีการเรียกตัวเองหรือ recursive นั่นเอง ซึ่งผมได้อธิบายแนวคิดของการเรียกตัวเองไว้แล้วในบทความ Functional Programming #2 : ฟังก์ชันพิสุทธิ์และการเรียกตัวเอง นิยามดังกล่าวสามารถถอดเป็นภาษา Python ได้ดังนี้ครับ

 ผมถอดนิยามคณิตศาสตร์ออกมาเป็นฟังก์ชันแบบตรงไปตรงมา เข้าใจไม่ยาก ผมตั้งใจที่จะซ้อนฟังก์ชัน binomial(n, r) ไว้ภายใต้ฟังก์ชัน B(m) เพราะอยากที่จะประกาศตัว B(m) สู่ภายนอกเพียงตัวเดียว ภาษาอย่าง Python ทำได้ครับ นับได้ว่าเป็นข้อดีมากข้อหนึ่ง ว่ากันไปแล้วมันก็คือฟังก์ชันมีระดับ ที่สามารถซ้อนฟังก์ชันเข้าไปได้นั่นเอง นิยามการเรียกตัวเองที่เห็นนี้ แท้ที่จริงแล้วเป็นการย้อนเกล็ดการหาค่าในบทความ ก้าวแรก HPC #04 : ตำนาน Note G นั่นเอง อาจจะลองไล่ย้อนดูได้ครับถ้าท่านสนใจ แต่นิยามเหล่านั้นไม่ใช่ประเด็นหลักที่ผมจะกล่าวถึงในบทความนี้ครับ ดังนั้นท่านที่ไม่สนใจเรื่องนิยามทางคณิตศาสตร์ก็ข้ามมันไปได้ ถือว่าเรามีอัลกอริทึมแล้วก็แล้วกัน 

การทดสอบอัลกอริทึม

รู้ได้อย่างไรว่าโปรแกรมถูกต้อง ไม่น่าถามใช่ไหมครับก็เอาโปรแกรมไปรันดู เทียบกับคำตอบที่มี ก็จะรู้เองว่าถูกต้องหรือไม่ใช่ไหมครับ แต่ก่อนจะทดสอบก็ต้องเขียนโปรแกรมขึ้นมาก่อน ผมมี source code แปะให้ในบทความนี้อยู่แล้วครับ วิธีพื้นฐานคือ สร้างแฟ้มใหม่แล้วพิมพ์ทีละตัวลอกจากหน้าเว็ปของผม ยุคเร่งรีบอย่างนี้ผมไม่คิดว่ามีใครจะทำนะครับ ส่วนมากใช้วิธีที่สองคือสร้างแฟ้มเปล่าไว้แล้วตัดแปะโปรแกรมจากหน้าเว็ปไปปล่อยเลย แบบนี้ก็ทำได้ครับ แต่วันนี้ผมขอเสนอวิธียอดนิยม สะดวกและไม่เสี่ยงว่าตัดแปะแล้วเพี้ยนทำงานไม่ได้ ภาษาบางภาษาถ้าเยื้องย่อหน้าเพี้ยนไปแม้เพียงเล็กน้อยก็พาลประท้วงไม่ยอมทำงานอย่างภาษา Python เป็นต้น หรือแม้กระทั่งโปรแกรมของผมมีการปรับปรุงในเครื่องคอมพิวเตอร์แล้ว แต่ไม่ได้ปรับปรุงที่เว็บ ท่านก็จะได้โปรแกรมรุ่นเก่าไป  วิธีที่ 3 ที่ผมอยากให้ท่านใช้แก้ปัญหาของสองวิธีแรกได้ครับ นั่นคือการใช้ Git ครับ

Git ห้านาที

ผมขอสองนาทีคุยเรื่องที่มาของ Git หน่อยครับ (ตามธรรมเนียม) คือ Kernel ของ Linux นั้นเขียนขึ้นจากภาษา C นับได้ว่ามีขนาดใหญ่มาก แค่รุ่นแรกก็เกือบสองแสนบรรทัดแล้ว ปัจจุบันก็เข้าใกล้ 20 ล้านบรรทัดเข้าทุกที เทียบแแล้วใหญ่ขึ้นประมาณ 100 เท่า ใหญ่ขนาดนี้เขียนคนเดียวไม่ได้หรอกครับ  กะเอาว่าน่าจะผัดเปลี่ยนหมุนเวียนกันหลายพันคนอยู่ แน่นอนครับ คนมากขนาดนี้จากทั่วโลก ร่วมกันเขียนโปรแกรมเดียวกัน มันเป็นเรื่องหอคอยบาเบลดีๆ นี่เองครับ ดังนั้นจึงต้องมีระบบจัดการ source code ที่ดีที่เราเรียกว่า version control system (VCS) ในอดีตสิบกว่าปีที่แล้วเขาใช้โปรแกรม BitKeeper แต่ตอนหลังเกิดปัญหาเรื่องลิขสิทธิ์กัน ผู้สร้าง Linux คือ Linus Tovalds เลยเขียน VCS ขึ้นมาเองทดแทน ตั้งชื่อว่า Git  ปัจจุบันทีมพัฒนา Git ก็มีหลายคน มีรุ่นใหม่มาเรื่อยๆ Git นั้นมีความพิเศษอย่างหนึ่งครับ คือฐานข้อมูลที่เก็บ source code ที่เขาเรียกว่า repository นั้นมีให้บริการบน cloud ครับ และบางบริษัทถ้าเรา opensouce โปรแกรมของเรา เราก็สามารถใช้บริการได้ฟรีเช่น Github หรือ Bitbucket เป็น  Github ก็เป็นผู้บริการรายใหญ่รายหนึ่งครับ  source code  เว็ปนี้ใช้บริการของ Github ในการเก็บ source code ทั้งหมด ดังนั้นผมแนะนำให้ท่านดึง Source code ของผมจาก Github ครับ เพราะท่านจะได้รับ source code รุ่นล่าสุดเสมอ ท่านสามารถ download source code ของทั้งเว็บนี้ใจาก Github ในขั้นตอนเดียวเลยเพียงท่านใช้ Git แต่ถ้าท่านยังไม่มี Git  สำหรับ windows ให้ไป download ที่ http://git-scm.com/downloads แล้ว install แต่ถ้าเป็น *buntu ให้พิมพ์

เมื่อท่านได้ Git แล้ว ไปที่ Shell ครับ ที่ไดก็ได้ ถ้าเป็น Windows  ไปที่ C:\ ยังได้เลยครับ แล้วพิมพ์

แค่นี้เองครับ Source code ของเว็ปนี้ทั้งเว็ปจะไปอยู่ที่เครื่องของท่าน  ท่านจะพบว่ามี directory HPCThai โผล่ขึ้นมา ภายในก็มี source code ทั้งหมดของเว็บนี้ (ซึ่งมีเพิ่มเติมและปรับปรุงเรื่อยๆ) ท่านลองเอา mouse ไปทาบบนส่วนหนึ่งส่วนใดของ source code ข้างบนครับ ท่านจะเห็นบรรทัดบนสุดที่บ่งชี้ตำแหน่งของแฟ้มในโครงสร้าง directory และชื่อแฟ้มครับ อย่างในกรณีนี้คือ HPCThai/introHPC/05/b_recur_v1.py   ท่านไล่ไปตามลำดับชั้นนี้ท่านจะพบแฟ้มที่ต้องการครับ

ตราบใดที่เว็บนี้ยังมีชีวิตอยู่ repository นี้จะได้รับการปรับปรุงเรื่อยๆ ครับ ถ้าท่านพิมพ์

พิมพ์ที่ใดก็ได้ในโครงสร้าง directory ที่สร้างขึ้นครับนี้เพื่อเป็นการ “refresh” repository ให้ทันสมัย ตรงกับรุ่นล่าสุดใน Github ทำได้บ่อยครั้งเท่าที่ต้องการ แต่แนะนำว่า ควรทำเมื่อเห็น source code ใหม่ๆ ในเว็ปที่ท่านไม่เคยเห็นมาก่อน แต่ก็ขอบอกไว้นะครับว่าใน repository นี้จะมักจะมี source code มากกว่าที่ท่านเห็นในเว็บนะครับ ไม่ต้องแปลกใจ เพราะผมเขียนเว็ปไว้ล่วงหน้า ดังนั้น repository จึงมีเท่าที่ผมเขียนไม่ใช่เท่าที่ผมแสดงให้ท่านเห็นครับ

ทดสอบความถูกต้องของโปรแกรมโดยใช้ ipython

เปลี่ยน directory ปัจจุบันไปยัง 05 ครับ เราจะทดสอบ b_recur_v1.py กันแล้ว จากนั้นให้พิมพ์ว่า

ถ้าไม่ชอบแบบ shell ก็สามารถใช้

จะเปิดมาเป็น windows mode เลย หรือถ้าชอบ web ก็สามารถใช้

ipython เป็น REPL ที่ใช้ทดแทน ‘python’ ที่เป็น REPL หลัก  แต่ถ้าเป็น *buntu ให้ใช้ ipython3 แทนนะครับ เพราะ ถ้าใช้ ipython จะเป็น python2

ipython มีลูกเล่นแพรวพราวกว่า python เยอะครับ วันนี้เราจะเรียนรู้ลูกเล่นเหล่านั้นบางส่วน เมื่อเข้าสู่ ipython ได้แล้วเราก็ load module ของงานครับ ไม่ต้องงงนะครับแฟ้มที่ .py ทั่วไปนั่นแหละครับคือ module ที่ว่า สำหรับของเราในนั้นมีฟังก์ชันเดียวครับคือ B การ load module เขียนดังนี้ครับ

บอกเผื่อไว้ครับ ถ้าต้องการออกจากโปรแกรม ipython ให้กดปุ่ม Ctrl-d ครับ เอาหละครับ ตอนนี้ ipython รู้จัก B แล้ว เรียกใช้งานได้เลยเช่น
B8เราใช้ REPL ของ ipython เพื่อหา ฺBernoulli ลำดับที่ 8 หรือ B(8) ซึ่งได้คำตอบเป็น Fraction ครับ ซึ่ง Fraction ก็คือเลขเศษส่วนนั่นเอง คำตอบคือ -1/30 ครับ ลองไปที่  http://en.wikipedia.org/wiki/Bernoulli_number ด้านขวาจะมีรายการ Bernoulli Number 20 ลำดับแรก ลองทดสอบดูได้ครับว่าผลลัพธ์ถูกต้องหรือไม่

การประเมินประสิทธิภาพ

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

คำสั่งขึ้นด้วย % เป็นเอกสิทธิ์ของ ipython ครับ ใน ‘python’ ไม่มี คำสั่งแรกที่เราจะเรียนรู้กันก็คือ %time เป็นการจับเวลานั่นเอง ซึ่งผลลัพธ์ที่ได้ คือ 48.6s หรือถ้าต้องการให้แม่นยำขึ้นก็ทำหลายรอบ แล้วมาเฉลี่ยค่ากัน ซึ่ง ipython ก็เตรียม %timeit ใช้แทน %time ครับ คำสั่งนี้จะตัดสินใจเองครับว่าจะทำกี่รอบที่ไม่นานจนเกินไปนัก ซึ่งว่ากันจริงๆ แล้ว %timeit จะเหมาะกับการทดสอบ module ที่กินเวลาน้อย ระดับ milliseconds ไปถึงไม่กี่วินาทีครับ การทดสอบข้างต้นนี้ผมเลยเลือก %time แทน  คราวนี้ลองเปลี่ยนมาเป็นบอร์ด banana pi  ที่ใช้ CPU A20 ดูบ้างจะได้

B20_bp_v1

เห็นชัดครับเมื่อเทียบกับกัน core เดียว A20 ช้า Core i7 3630 (2.4GHz) อยู่ประมาณ 6 เท่า

การเปรียบเทียบเวลาที่ใช้หา B(m) เมื่อไล่ค่า m

ผมรู้สึกตะหงิดๆ ว่าค่า m ยิ่งเยอะยิ่งเวลาที่ใช้ในการหาจะเพิ่มเป็นทวีคูณ เพื่อพิสูจน์ทฤษฏีของผม ผมจะเขียนโปรแกมทดสอบ โดยที่โปรแกรมทดสอบนี้ไปเรียกใช้ B อีกที  เพื่อให้โปรแกรมนี้สามารถใช้งานกับอัลกอริทึมที่มีการปรับปรุงในโอกาสต่อๆ  ผมจึงเขียนโปรแกรมในส่วนนี้ให้ยืดหยุ่นหน่อย ดังนี้ครับ

โปรแกรมนี้ใช้ทั้ง numpy และ matplotlib และเนื่องจากต้องเอาเวลาไปแสดงผลเป็นกราฟเองดังนั้นจึงไม่สามารถใช้ %time ได้ โปรแกรมต้องมีกลไกในการเก็บสะสมเวลาเอง ซึ่งการจับเวลาจะอยู่บรรทัด 12 และ 14 ครับ โปรแกรมนี้ใช้ matplotlib ในการสร้างกราฟ อยู่ที่บรรทัด 23-30 ใครเคยใช้ MATLAB จะรู้สึกเลยครับว่าคล้ายกันมาก ยิ่งถ้าผมใช้ pylab ก็จะยิ่งคล้ายกว่านี้ครับ และผมเลือกที่ไม่แสดงกราฟที่เป็นผลลัพธ์ออกมาหน้าจอเพราะไม่ใช่ทุกเครื่องมีหน้าจอนั่นเอง (ฺบอร์ด banana pi ที่ผมใช้ไม่ได้ต่อหน้าจอครับ) ผมจึงเก็บเป็นแฟ้ม .png แทน รายละเอียดที่เหลือแกะจาก source code ดูได้ครับ ลองทดลองรันโปรแกรมเลยดีกว่าครับ

ผลลัพธ์ที่ได้ดังนี้

r1 คืออัลกอริทึมแบบ Recursive version 1 ครับ ตอนนี้มีแค่ version เดียว เดี๋ยวต่อไปจะมี r2  ส่วน 20 ก็คือค่าทำไปจนถึง B(20) นั่นเอง เมื่อรันโปรแกรมเสร็จ ระบบจะสร้างแฟ้ม b_recur_v1.png มาให้ด้วยครับเป็นกราฟ และ b_recur_v1.out เป็นข้อมูลดิบของกราฟจะนำไปใช้ในขั้นตอนต่อไป ใส่ส่วนของกราฟนั้นเปิดดูได้ครับ ถ้าเป็น windows ไม่ต้องไปหาแฟ้มเปิดที่ไหนครับ พิมพ์

Windows จะไปหาโปรแกรมที่เหมาะสมมาเปิดให้เอง แต่ถ้าเป็น *buntu ให้ใช้

ผลลัพธ์ที่ได้จะประมาณนี้ครับ

b_v1_graph

จะเห็นได้ว่ากราฟชันขึ้นเรื่อยๆ นั่นหมายความว่า ค่า m ใน B(m) ยิ่งสูงขึ้นเท่าไหร่ก็ยิ่งใช้เวลามากขึ้นเป็นทวีคูณครับ

การทำ Curve Fitting

การทำ Curve Fitting คือการหาสมการมาสมการหนึ่งซึ่งเมื่อใส่ค่าในแกน x เข้าไป จะสามารถให้ผลลัพธ์ใกล้เคียงกับค่าในแกน Y ที่เก็บข้อมูลมา ในกรณีนี้คือแกนเวลา พูดง่ายๆ ก็คือหาสมการที่สามารถหาค่าให้ทับทุกค่าที่เก็บมาหรือใกล้เคียงที่สุด เอา source code นี้ไปใช้เลยครับ ดังนี้ครับ

 โปรแกรมนี้จะให้เราป้อนสองค่าครับ ค่าแรกคือชื่อแฟ้มข้อมูลดิบที่ได้จากโปรแกรม bern.py นั่นเอง ส่วนค่าที่สองเป็นตัวเลขประมาณการสูงสุด ซึ่งคืออะไรช่างมันก่อนครับ ในที่นี้ใส่ค่า 30 ไปก่อน ดังนั้นเวลาใช้งานก็จะเป็น

จะได้ผลลัพธ์

curvefit_v1

เรามาพิจารณาแฟ้ม b_recur_v1_fit.png ดูครับ จะได้ภาพนี้

b_recur_v1_fit

(แก้ไข: ที่หัวกราฟสะกดผิดครับที่ถูกต้องเป็น polynomial)

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

เกิดมีใครอยากหา B(30) โดยใช้สมการนี้ ท่านจะว่าอย่างไรครับ ปล่อยให้หาไปอย่างนั้นหรือ จะเสร็จเมื่อใด 10 นาทีเสร็จไหม หรือ 5 ชั่วโมงกันแน่ สิ่งที่ท่านทำได้ก็คือลองดูก่อน สู้กันครับว่าค่าจะออกมาก่อนหรือท่านจะหมดความอดทนก่อน มันคงจะดีกว่าถ้าเรารู้ล่วงหน้า เช่นถ้าเรารู้ล่วงหน้าว่ามันจะใช้เวลารัน 3 วัน เราก็จะได้ไม่สั่งให้มันรัน รอไม่ไหวอยู่แล้ว เพราะเหตุนี้เองจึงเป็นที่มาของการหาสมการตัวนี้ ทำไมหรือครับ ก็สมการตัวนี้ไม่ได้แค่ลากเส้นสีน้ำเงินให้เท่าที่ท่านเห็นในกราฟนะครับ แต่หากมันสามารถลากต่อออกนอกตารางได้ นั่นคือการพยากรณ์นั่นเอง เป็นอีกศาสตร์หนึ่งทางสถิติว่าด้วยเรื่องของการวิเคราะห์แนวโน้ม ผมทำมาให้ท่านเสร็จเรียบร้อยแล้วครับ จำเลข 30 ที่ใส่เป็นพารามิเตอร์ไหมครับ เป็นการบอกให้มันสร้างกราฟแนวโน้มให้ถึง B(30) นั่นเอง คำตอบมันอยู่ในเครื่องของท่านแล้วครับ แฟ้ม b_recur_v1_est.png โปรแกรมสร้างให้พร้อมกับกราฟแรกนั่นเอง ผลลัพธ์จะได้ประมาณนี้ครับ

b_recur_v1_est

จากกราฟจะเห็นค่าพยากรณ์ว่าประมาณ 17,000 วินาทีใช่ไหมครับ หรือ 4.7 ชั่วโมงนั่นเองรอไหวไหมครับ เป็นผมผมไม่รอครับ ผมถามเล่นๆ ดูจากกราฟพยากรณ์นี้ ท่านว่าถ้าหาค่า B(100) จะใช้เวลาประมาณเท่าไหร่ครับ ให้ทายได้คนละครั้ง เชื่อว่าส่วนมากทายน้อยกว่าที่ทำนาย คำตอบคือ 1158 ปีครับ บน Core i7 แล้วนะครับ เกิดใหม่ร้อยชาติก็ยังไม่เสร็จเลยครับ แกน Bm เปลี่ยนเล็กน้อยแต่แกนเวลาวิ่งกันเป็นทวีคูณ ท่านก็อาจจะสงสัยอีกว่า ตัวเลขที่ได้นี้เชื่อถือได้แค่ไหน เห็นผลลัพธ์ที่ได้จากโปรแกรมล่าสุดไหมครับ ทำนาย B(22)  จะใช้เวลาประมาณ 190 วินาที ผมลองใน ipython ได้ผลลัพธ์

B22_nb_v1

3:14m แปลงเป็นวินาทีได้ 194s เมื่อเทียบกับ 190s ที่พยากรณ์ไว้ถือได้ว่าใกล้เคียงกันมากครับ ผิดพลาดไปเพียงแค่ 2% ถือว่ายอมรับได้ครับ

ย้อนกลับมา Banana Pi

ที่ผ่านมาผมแสดงให้เห็นถึงการทำงานบน notebook และ banana pi โดยเน้นการทำงานที่ notebook เป็นหลัก เพื่อเปรียบเทียบให้เห็นความแตกต่างของ CPU ทั้งสองเครื่อง ก็น่าจะเห็นภาพเปรียบเทียบได้คร่าวๆ แล้วนะครับ แต่นี้ต่อไปผมคงตัด notebook ออกไปสักพัก เพื่อลดขนาดของบทความไม่ให้ใหญ่จนเกินไป เพราะเป้าหมายที่แท้จริงของไม่ได้อยู่ที่การเปรียบเทียบความสามารถของ CPU แต่หากอยู่ที่การตั้งฐานอ้างอิงในการรันรอบแรกและพยายามหาวิธีทำให้ได้ผลลัพธ์ที่เร็วขึ้นต่างหาก ดังนั้นผมขอใช้ banapa pi เป็นหลัก สาเหตุหลักก็คือมันนิ่งกว่าครับ เวลาจับเวลาตัวเลขจะน่าเชื่อถือกว่า เพราะไม่มี process มากมายอยู่ในเครื่องเหมือนใน Windows แม้แต่ GUI ก็ยังไม่มีเลยครับ ผมเลยขอตามรอยการทำงานข้างต้นแต่ใช้ banana pi ซึ่งได้ผลลัพธ์ดังนี้ครับ

B20_bp_v1

B_bp_v1

B_curvefit_bp_v1

b_recur_v1_bp_fit

b_recur_v1_bp_est1

B22_bp_v1

สรุปได้ว่า B(22) ประเมินว่า ใช้ 1049s ของจริงใช้ 1545s (25:35m) ซึ่งคลาดเคลื่อนประมาณ 30% เลขไม่ค่อยสวยแต่ไม่เป็นไรครับ ยังยอมรับได้ และประเมินว่า B(100) จะใช้เวลา 1,319,055,732,707s หรือ 41,826 ปี  นานมากครับ เวลานั้นมนุษย์โลกคงสูญพันธุ์ไปหมดแล้ว

การปรับปรุง v2

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

ท่านเคยฉุกคิดหรือไม่ครับว่าทำไมผมจึงหา Bernoulli Number แต่เลขคู่ไม่หาเลขคี่ จริงๆ แล้วนอกจาก B(1) ที่มีค่า 1/2 แล้ว B(เลขคื่) อื่นๆ มีค่าเป็น 0 ทั้งหมดครับ เดาได้ว่ามันต้องมี 0 ไปคูณไม่ขั้นตอนใดก็ขั้นตอนหนึ่ง ซึ่งเราไม่จำเป็นต้องสนใจ เราสนแต่เพียงว่าการหาค่า B เลขคี่ใช้เวลาใกล้เคียงกับ B เลขคู่ที่อยู่ติดกันแม้ว่าเรารู้คำตอบล่วงหน้าแล้วว่าเป็น 0 เสียเวลาโดยใช่เหตุครับ  ดังนั้นถ้าเราปรับอัลกอริทึม เมื่อส่งค่าเป็นเลขคี่ยกเว้น 1 มาให้ตอบกลับไปว่า 0 เลยไม่ต้องคำนวณ โปรแกรมเป็นดังนี้ครับ

ผมเพิ่ม if ง่ายๆ 2 ตัว ตั้งแต่บรรทัด 10 ถึง 13 ท่านคิดว่าโปรแกรมเร็วขึ้นไหม ถ้าเร็วจะเร็วขึ้นเท่าใด บางท่านคิดไวตอบได้ทันทีว่า 2 เท่า ก็เพราะเลขคี่ไม่ต้องทำแล้ว งานก็หายไปครึ่งหนึ่ง มาลองดูกันนะครับว่ามันจะจริงไหม B_bp_v2 จาก B(22) เดิมที่เคยใช้เวลา 25:35m เหลือเพียง 1.17s  และ B(30) ที่เคยประเมินว่า 41,826 ปี เหลือเพียง 17.9s ไม่น่าเชื่อใช่ไหมครับ คำตอบที่ได้ก็ถูกต้องไม่ได้ผิดแต่อย่างใด ท่านสามารถไปพิสูจน์โดยพิมพ์bernoulli(30) ใน sagemath.org จะได้เลขเดียวกันนี้ ทำไมถึงเป็นเช่นนี้ แค่ใส่ if ไปแค่สองตัว ทำไมผลลัพธ์มันต่างกันราวฟ้ากับดินขนาดนั้น ลองนึกภาพดูนะครับ ถ้าทำ B(20) อยู่ ก็ต้องมีการอ้างอิงถึง B(19) ถ้า B(19) ไม่ต้องทำเลยก็จะประหยัดยกยวงที่เป็นหางว่าวของ B(19) มาทำ B(18) ซึ่ง B(18) ก็ไม่ต้องทำ B(17) ประหยัดไปอีกยวง แบบนี้มันทวีหารชัดๆ (ผมนิยามคำเองครับ) เลยลดฮวบของฮวบเลย สมมติว่าท่านมี CPU หลายตัว และสามารถให้มันช่วยกันทำได้โดยไม่เสียเวลา overhead ใดๆ ทั้งสิ้น และอัลกอริทึม v1 นี้สามารถกระจายงานให้ทุก CPU ทำได้เท่าๆ กันอย่างเต็มที่  ท่านต้องหา CPU มา 41,826 ตัวและเปิดเครื่องทุกเครื่องทำงานไป 1 ปี จึงจะสามารถหาคำตอบของ B(30) ได้ แต่แค่เติม if ไปสองตัว ใช้ CPU ตัวเดียว หาคำตอบได้ภายใน 17.9 วินาที เรื่องนี้บอกอะไรเรา ใช่แล้วครับอัลกอริทึมสำคัญกว่า hardware ครับ ถ้าอัลกอริทึมแย่แต่ hardware ดีแค่ไหน ก็แทบไม่มีประโยชน์ครับ ดังนั้นสิ่งที่ต้องทำในงาน HPC คือต้องปรับแต่งโปรแกรมของเราให้ทำงานได้เร็วที่สุดก่อนค่อยคิดขยายไปทำงานหลาย CPU ครับ ทำย้อนศรไม่ได้ครับ คราวนี้เรามาลองหาถึง B(36) บ้างครับ

B_list_v2 จากรูปจะเห็นว่าเราสามารถหา B(36) ได้โดยใช้เวลาเพียง 2 นาทีกว่าเท่านั้น นับว่าเร็วขึ้นมาก ถ้าลองพิจารณาชุดตัวเลข จะเห็นบางอย่างน่าสนใจนะครับ คือ B แต่ละตัวที่เพิ่มขึ้นดูเหมือนว่าจะใช้เวลาเพิ่มขึ้นเป็น 2 เท่า ถ้าแบบนี้คิดสูตรเองได้เลยครับ นั่นคือ T_{B_{m}}=2^{\frac{m-22}{2}} อันนี้แบบหยาบๆ นะครับ ถ้าจะให้แม่นกว่านี้ก็ต้องใช้ curve fitting จะได้ผลลัพธ์ดังนี้ครับB_curvefit_bp_v2 เห็นไหมครับว่าโปรแกรมพยากรณ์ว่า B(38) นั้นจะใช้เวลา 262s เป็นสองเท่าของ B(36) ตามที่คาดไว้เลยครับ แต่ครั้นจะหา B(100) ก็ยังใช้เวลาเกือบ 57 ปีอยู่ดีครับ

การปรับปรุง v3

ยังมีอะไรให้ปรับปรุงอีกไหมครับ ยังครับยังไม่หมด ยังไปได้ต่อ ผมอยากให้คุณพิจารณาสูตรในการหาค่า Bernoulli ให้ดีนะครับ ผมขอยกมาอีกครั้งหนึ่งเพื่อให้เห็นได้ชัดเจน ดังนี้ครับ B_{0}=1 B_{m}=1-\sum_{k=0}^{m-1}\binom{m}{k}\frac{B_{k}}{m-k+1} มองสูตรแล้วอยากให้จินตการตามนะครับ  ถ้าเราต้องการหา B(20) ซึ่ง B(20) ต้องการ B(19) ซึ่งเป็น 0 ไม่สนใจ ต่อมาก็ต้องการ B(18) ตรงนี้เกิดการขยายกิ่งแล้วครับ แยกไปหา B(18) สายหนึ่ง ส่วนตัว B(20) เองก็ต้องไปหา B(16) ถามว่ากิ่งที่แยกไปของ B(18) นั้นต้องการ B(16) หรือไม่ แน่นอนครับ กลายเป็น B(16) เกิดการหา 2 หน และใน B(16) ก็ใช่ย่อยนะครับ ยังหา B เล็กๆ อีกจำนวนมาก เกิดทวีคูณอีกแล้วครับ B(18) 1 ครั้ง B(16) 2 ครั้ง B(14) 4 ครั้ง B(12) 8 ครั้ง B(10) 16 ครั้ง B(8) 32 ครั้ง B(4) 64 ครั้ง B(2) 128 ครั้ง จึงไม่น่าแปลกใจครับว่าทำไมมันกินเวลามากขนาดนั้น มันเหมือนเราเล่นเกมส์คอมพิวเตอร์ พอผ่านด่าน 4 ก็ต้องมาเล่นด่าน 1, 2, 3, 4 ใหม่ ก่อนที่จะไปด่าน 5 ได้ แบบนี้คงไม่มีใครเล่นแน่ครับ ใครทำเกมส์แบบนี้ออกมาขายก็ต้องเจ๊งแน่ครับ ดังนั้นโปรแกรมของเรา ถ้าเราทำให้มันทำงานแค่เพียงครั้งเดียวในแต่ละ B ได้ มันจะเร็วขึ้นขนาดไหนนี่  ผมปรับโปรแกรมเป็น v3 ดังนี้ครับ

หลักการนี้เรียกว่า caching นั่นเองครับ หลักการง่ายๆ ครับคือผมสร้าง dictionary ขึ้นมาตัวหนึ่ง (หรือบางภาษาเรียกว่า hash table)  เช่นถ้าส่ง m เข้ามาใน B ก็ตรวจสอบดูว่า m นั้นมีใน dictionary แล้วรึยัง ถ้ามีแล้วก็ไม่ต้องคำนวณครับ เอาค่าจากใน dictionary ไปใช้เลย แต่ถ้ายังไม่มี ก็เข้าสู่กระบวนการคำนวณตามปกติ เมื่อได้ผลลัพธ์แล้วก่อนส่งผลลัพธ์นั้นให้แก่ผู้ร้องขอ โดยที่แอบเก็บเพิ่มไว้ใน Dictionary ก่อน เท่านี้เองครับ ในส่วนของ binomial ก็เช่นกันครับ ใช้ dictionary อีกตัวเพื่อเก็บ ท่านคิดว่าจะเร็วขึ้นแค่ไหนครับ มาลองดูกัน

b_list36_v3

โอ้ว…พระเจ้าช่วยกล้วยทอดมันๆ อะไรมันจะติดปีกขนาดนั้น นี่ทำงานบน Banana Pi นะครับ ไม่ใช่ Supercomputer แต่อย่างใด อย่ากระนั้นเลย ลองถึง B(1000) เลยดีกว่า

b_list1000_v3

โอ้โหขนาด B(1000) ยังหาได้ภายใน 8s เองหรือนี่ และดูตัวเลขที่เพิมขึ้นสิครับ การเพิ่มเกือบๆ เป็นเส้นตรงเลย ดูดีมากเลยใช่ไหมครับ แต่เดี๋ยวก่อนครับ เรื่องนี้มีเงื่อนงำ มันไม่ได้ดีอย่างที่เห็นหรอกครับ B(1000) ไม่ได้สามารถหาได้ภายใน 8 วินาทีอย่างที่แสดง เรื่องนี้ต้องระวังให้มากครับ มีงานวิจัยขนาดใหญ่หลายงานเลยที่ตกม้าตายตอนคำนวณลักษณะนี้ครับ เรากำลังเปรียบเทียบโดยใช้ฐานคนละฐานครับ B(1000) ในวิธีก่อนหน้านั้นเป็นการหาโดยอิสระครับ หมายความว่าคุณเข้า REPL ของ ipython แล้วพิมพ์หา B(1000)  เวลาที่ได้ก็ตามนั้นครับ แต่วิธีนี้ B(998) นั้นเขาหา B ตัวที่ถัดลงไปเรียบร้อยหมดทุกค่าแล้วครับ  Cache เก็บไว้แล้ว เมื่อร้องขอค่าข้างล่าง ก็ไม่ต้องไปหาอีก หาแต่ตัวมันเองแล้วครับ ดังนั้นเมื่อเข้า REPL แล้วหา B(1000) รับรองได้ว่าวิธีที่ 3 นี้ไม่สามารถหาได้ภายใน 8.83s แน่นอนครับ แล้วจะใช้เวลาเท่าไหร่ บอกได้เลยครับต้องเอาเวลาทั้งหมดก่อนหน้านั้นบวกรวมกันทั้งหมดครับ ถึงจะเป็นของจริง รูปข้างล่างนี้คงพอจะอธิบายภาพให้เห็นได้ชัดครับ ลองพิจารณาดู

b_v3_solve

B(300) ใช้เวลา 17.6s ซึ่งก็นับว่าเร็วมากนะครับเมื่อเทียบกับสองวิธีที่ผ่านมา ซึ่ง v2 นั้นใช้เวลาหาแค่ B(100) ก็เกือบ 57 ปีแล้ว แต่ลองดูครับ ทำไม B(302) ซึ่งใช้เวลาแค่กระพริบตาเท่านั้น ก็อย่างที่อธิบายไว้ครับว่า B(300) มัน cache ข้อมูลมาให้หมดแล้ว B(302) เลยแทบไม่ต้องทำอะไร ก็ได้คำตอบแล้วครับ ถ้าเริ่มต้นมายังไม่ทำอะไรเลย แล้วหา B(302) สิครับ ลองดูครับ

b_v3_solve2

คงไม่ต้องอธิบายอะไรเพิ่มแล้วนะครับ ปิดคดี ผมขอเก็บ B(1000) เป็นตัวอ้างอิงใหม่ เอาใช้กับบทความหน้านะครับ แต่ผมไม่แสดงค่าคำตอบเพราะยาวเป็นหน้าเลย (แก้ไข: ไม่รู้ว่าตอนเขียนบทความตอนนั้นคิดยังไงถึงไปตัดเอาคำตอบออก ผมไม่อยากแก้แล้ว ใครอยากรู้ว่า B(1000) ได้ผลลัพธ์เท่าไร ดูในส่วน C++ ข้างล่างครับ)

b_1000_v3

 C++ กับหลักฐานที่น่าสนใจ

ผ่านไปแล้วนะครับ Python คราวนี้มาถึงตา C++ บ้าง เนื่องจากบทความนี้เริ่มจะยาวแล้ว ผมเลยขอเขียนโดยสังเขปครับ ไปดูกันที่ v3 เลย ในงานจริงก็จะอารมณ์ประมาณนี้เหมือนกันครับ ใช้ Python เป็น prototype เมื่อเขียนเรียบร้อยเป็นที่น่าพอใจแล้วก็มาเขียนด้วย C++ แทน เมื่อสำรวจดูแล้ว C++ มาตรฐานไม่สามารถรองรับตัวเลขจำนวนเต็มที่ยาวได้เรื่อยๆ หรือที่เรียกว่า multiple precision number ได้ จำเป็นต้องหา library เสริมมาใช้ครับ ที่นิยมใช้คือ GNU Multiple Precision Arithmetic Library หรือ GMP ซึ่งเรานำมาลง *buntu ได้ดังนี้ครับ

หรือถ้านเป็น MinGW บน Windows ใช้ใช้

เมื่อติดตั้งเรียบร้อยแล้วก็เขียนภาษา C++ ดังนี้ครับ

เมื่อเขียนเสร็จคอมไพล์และรันโดยใช้ใช้คำสั่ง

จะได้ผลลัพธ์บน Banana Pi ดังนี้ครับ

b_bp_cpp

ค่าทุกค่าถูกต้องทั้งหมดครับ รวมทั้ง B(1000) ด้วย ผมตรวจสอบกับ sagemath.org เรียบร้อยแล้ว ขอให้จำในส่วนหัวคือ -182 กับตอนท้ายคือ 78901 เอาไว้นะครับ นี่คือตัวเลขที่ถูกต้อง พอจับเวลาบรรทัดสุดท้ายเฉพาะ B(1000) ได้เรื่องครับ กินเวลา 1300s หรือ 21:40m ซึ่งช้ากว่า Python ที่กินเวลาเพียง 19:31m  เกิดอะไรขึ้นครับ ทำไมภาษา C++ ถึงช้ากว่า Python ได้ ทั้งๆ ที่ C++ เป็นภาษาแบบคอมไพล์น่าจะเร็วกว่า อัลกอริทึมก็ตัวเดียวกัน และจริงๆแล้วโปรแกรมของ C++ ก็มีแอบโกงเล็กๆ อยู่ครับคือ มันแอบ cache ถึง B(20) ไว้ก่อนหน้าที่คำนวณ B(1000) แต่ก็นับได้ว่าเล็กน้อยมากครับแทบไม่มีผล เพื่อความสะดวกผมเลยรวบเป็นโปรแกรมเดียวกัน แต่ขนาดต่อให้ขนาดนี้แล้วยังช้ากว่า น่าสนใจมากครับ ฝากไว้ให้คิดเป็นการบ้านดีกว่าครับ บทความยาวเกินแล้วครับ ขอยกยอดไปเฉลยกันในบทความต่อไปครับ

ทิ้งท้าย

บทความนี้นำเสนออัลกอริทึมอย่างง่ายและการปรับแต่งประสิทธิภาพเบื้องต้นเพื่อการหาค่า Bernoulli Number เราอาจปรับโปรแกรมเล็กน้อย ไม่ต้องใช้ recursive ก็ได้ เอาผลลัพธ์ B ที่หาได้ไปเก็บไว้ใน array เพื่อเป็นฐานให้หาค่า B ที่ใหญ่ขึ้น ผมเชื่อว่าอัลกอริทึมที่ผมเพิ่งว่าไปนี้น่าจะเป็นอัลกอริทึมที่ Ada ใช้ในการสร้างโปรแกรมแรกของโลก แต่นั่นก็เป็นแค่เรื่องคาดเดาครับ ผมยังไม่มีเวลาศึกษาเอกสารฉบับนั้นจริงๆ

พบกันใหม่ในบทความต่อไปครับ

[Total: 2    Average: 2/5]

You may also like...

Leave a Reply

Your email address will not be published. Required fields are marked *