ฉันเกลียด NumPy

• ผู้เขียนอธิบายปัญหาของ NumPy พร้อมตัวอย่างหลายกรณี
• การคำนวณอาร์เรย์แบบง่ายทำได้สะดวกด้วย NumPy แต่เมื่อจำนวนมิติเพิ่มขึ้น ความซับซ้อนและความสับสนก็เพิ่มขึ้นอย่างรวดเร็ว
• การออกแบบของ NumPy เช่น broadcasting และ advanced indexing ยังขาดความชัดเจนและการทำ abstraction ที่ดี
• การเขียนโค้ดจึงต้องพึ่งการเดาและการลองผิดลองถูก แทนที่จะระบุแกนอย่างชัดเจน
• ผู้เขียนเสนอแนวคิดเกี่ยวกับภาษาอาร์เรย์ที่ปรับปรุงแล้ว และจะนำเสนอทางเลือกที่เป็นรูปธรรมในบทความถัดไป

---

บทนำ: ทั้งรักทั้งเกลียด NumPy

• ผู้เขียนระบุว่าใช้ NumPy มานาน แต่ก็ผิดหวังกับข้อจำกัดของมันอย่างมาก
• NumPy เป็นไลบรารีสำคัญและทรงอิทธิพลสำหรับการคำนวณอาร์เรย์ใน Python
• แม้แต่ ไลบรารีแมชชีนเลิร์นนิง สมัยใหม่อย่าง PyTorch ก็ยังมีปัญหาคล้าย NumPy อยู่ด้วย

สิ่งที่ง่ายและสิ่งที่ยากใน NumPy

• การคำนวณง่าย ๆ เช่นการแก้สมการเชิงเส้นพื้นฐาน ทำได้ด้วย ไวยากรณ์ที่ชัดเจนและสวยงาม
• แต่เมื่อมิติของอาร์เรย์สูงขึ้นหรือการคำนวณซับซ้อนขึ้น ก็จำเป็นต้องประมวลผลแบบเป็นชุดโดยไม่ใช้ for loop
• ในสภาพแวดล้อมที่ใช้ลูปไม่ได้ เช่นการคำนวณบน GPU จำเป็นต้องใช้ ไวยากรณ์ vectorization แบบเฉพาะทาง หรือรูปแบบการเรียกฟังก์ชันพิเศษ
• แต่ วิธีใช้งานที่ถูกต้อง ของฟังก์ชันเหล่านี้กลับคลุมเครือ และยากจะเข้าใจให้ชัดเจนจากเอกสารเพียงอย่างเดียว
• ตัวอย่างเช่น ในกรณีอาร์เรย์หลายมิติ ก็แทบไม่มีใครมั่นใจนักว่าควรใช้ฟังก์ชัน linalg.solve ของ numpy อย่างไรจึงจะถูกต้อง

ปัญหาของ NumPy

• NumPy ขาดทฤษฎีที่สอดคล้องกันสำหรับการนำการคำนวณไปใช้กับบางส่วนหรือบางแกนของ อาร์เรย์หลายมิติ
• เมื่อมิติของอาร์เรย์ไม่เกิน 2 ทุกอย่างยังค่อนข้างชัดเจน แต่พอเป็น 3 มิติขึ้นไป การระบุว่าแต่ละอาร์เรย์จะคำนวณกับแกนใดกลับไม่ชัดเจน
• จึงบังคับให้ใช้วิธีซับซ้อนอย่าง None, broadcasting, np.tensordot เพื่อให้มิติตรงกันอย่างชัดแจ้ง
• วิธีเหล่านี้ทำให้เกิดความผิดพลาดได้ง่าย ลดความอ่านง่ายของโค้ด และเพิ่มโอกาสเกิดบั๊ก

ลูปกับความชัดเจน

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

np.einsum: ฟังก์ชันที่ดีเป็นพิเศษ

• np.einsum ทรงพลังเพราะให้ภาษาเฉพาะโดเมนที่ยืดหยุ่นสำหรับการตั้งชื่อแกน
• einsum ทำให้เจตนาของการคำนวณชัดเจนและยัง generalized ได้ดี จึงสามารถเขียนการคำนวณแกนที่ซับซ้อนได้อย่าง explicit
• แต่แนวทางแบบ einsum กลับรองรับเพียง การคำนวณบางประเภท เท่านั้น และใช้กับ linalg.solve ไม่ได้ เป็นต้น

ปัญหาของ broadcasting

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

ความไม่ชัดเจนของการทำ indexing

• advanced indexing ของ NumPy ทำให้การคาดเดา shape ของอาร์เรย์เป็นเรื่องยากและไม่ชัดเจนมาก
• shape ของอาร์เรย์ผลลัพธ์จะเปลี่ยนไปตามการผสมรูปแบบ indexing ที่ใช้ จึงคาดเดาได้ยากหากไม่มีประสบการณ์ตรง
• แม้แต่เอกสารอธิบายกฎของ indexing ก็ยังยาวและซับซ้อน ทำให้เสียเวลาศึกษาอย่างมาก
• ต่อให้ตั้งใจจะใช้เพียง indexing แบบง่าย ในบางการคำนวณก็หลีกเลี่ยงการใช้ advanced indexing ไม่ได้

ข้อจำกัดของการออกแบบฟังก์ชัน NumPy

• ฟังก์ชันจำนวนมากของ NumPy ถูกออกแบบมาให้เหมาะกับ shape ของอาร์เรย์บางแบบเท่านั้น
• สำหรับอาร์เรย์หลายมิติ จำเป็นต้องใช้อาร์กิวเมนต์ axes เพิ่มเติม ใช้ชื่อฟังก์ชันแยก หรืออาศัยขนบเฉพาะ ซึ่งแต่ละฟังก์ชันก็ไม่สอดคล้องกัน
• โครงสร้างแบบนี้สวนทางกับหลักการเขียนโปรแกรมที่ให้ความสำคัญกับ abstraction และการนำกลับมาใช้ซ้ำ
• แม้จะมีฟังก์ชันแก้ปัญหาเฉพาะหน้าได้ แต่หากต้องนำกลับไปใช้กับอาร์เรย์และแกนแบบอื่น ก็ต้องเขียนโค้ดใหม่แทบทั้งหมด

ตัวอย่างจริง: การทำ self-attention

• เมื่อต้องเขียน self-attention ด้วย NumPy หากใช้ลูปโค้ดจะชัดเจน แต่ถ้าถูกบังคับให้ทำ vectorization โค้ดจะซับซ้อนขึ้นมาก
• โดยเฉพาะเมื่อเป็น multi-head attention ที่ต้องใช้การคำนวณหลายมิติ ก็ต้องผสมทั้ง einsum และการแปลงแกนเข้าด้วยกัน ทำให้โค้ดเข้าใจยาก

บทสรุปและทางเลือก

• ผู้เขียนระบุว่า NumPy เป็น "ตัวเลือกเดียวที่ทั้งมีจุดด้อยมากกว่า array language อื่น ๆ แต่ก็กลายเป็นสิ่งสำคัญในตลาดมากที่สุด"
• เพื่อแก้ปัญหาต่าง ๆ ของ NumPy เช่น broadcasting, ความไม่ชัดเจนของ indexing, และความไม่สอดคล้องของฟังก์ชัน ผู้เขียนบอกเป็นนัยว่าได้สร้างต้นแบบของ ภาษาอาร์เรย์ที่ปรับปรุงแล้ว ขึ้นมา
• แนวทางปรับปรุงที่เป็นรูปธรรม (API ของภาษาอาร์เรย์แบบใหม่) มีแผนจะนำเสนอในบทความแยกภายหลัง

[youn17]

ฟังดูเหมือนเป็นเรื่องราวว่าทำไม Julia ถึงถือกำเนิดขึ้นนะครับ แม้จะต้องศึกษาไลบรารีต่าง ๆ แต่ในแง่ที่มันช่วยแก้ปัญหาหลายอย่างของ NumPy ได้ ก็ดูเป็นตัวเลือกที่น่าสนใจมากจริง ๆ ครับ

[ahwjdekf]

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

[domino]

ผมไม่ชอบ NumPy

ดูเหมือนไลบรารี Python ที่ค่อนข้างเก่า ๆ จะมีปัญหาคล้าย ๆ กันหมด

[GN⁺]

ความคิดเห็นจาก Hacker News

• ในตัวอย่างแรก ถ้าดูเอกสารจากแค่ชนิดของ b จะอ่านเข้าใจยาก แต่เพราะมีคำอธิบายเกี่ยวกับ shape ที่คืนกลับมาด้วย จึงต้องตรวจสอบว่าเวกเตอร์ b นั้นมีรูปเป็นเมทริกซ์จริงหรือไม่ (โดยเฉพาะกรณี K=1)
• ถ้าอาร์เรย์มีมิติมากกว่า 2 มิติ แนะนำให้ใช้ Xarray ซึ่งเพิ่มชื่อมิติให้กับอาร์เรย์ NumPy เพราะการ broadcast/จัดแนวจะทำให้อัตโนมัติโดยไม่ต้องคอยจัดมิติหรือ transpose เอง จึงช่วยแก้ปัญหาพวกนี้ได้เกือบทั้งหมด แม้ Xarray จะอ่อนกว่า NumPy ในแง่ linear algebra แต่ก็กลับไปใช้ NumPy ได้ง่ายและแค่ทำฟังก์ชันช่วยไว้เล็กน้อยเท่านั้น พอใช้ Xarray แล้วประสิทธิภาพการทำงานกับข้อมูล 3 มิติขึ้นไปจะดีขึ้นมาก
  • Xarray ให้ความรู้สึกเหมือนรวมข้อดีของ Pandas และ NumPy เข้าด้วยกัน การ indexing แบบ da.sel(x=some_x).isel(t=-1).mean(["y", "z"]) ทำได้ง่าย และการ broadcast ก็ชัดเจนเพราะเคารพชื่อมิติด้วย เด่นมากสำหรับงานข้อมูลภูมิสารสนเทศที่มีหลาย CRS และยังใช้งานร่วมกับ Arviz ได้ดีมาก ทำให้จัดการมิติเพิ่มเติมในงานวิเคราะห์แบบเบย์ได้ง่าย อีกทั้งยังรวมหลายอาร์เรย์ไว้ใน dataset เดียวและใช้พิกัดร่วมกันได้ จึงสั่งงานกับทุกอาร์เรย์ที่มีแกนเวลาได้ง่าย เช่น ds.isel(t=-1)
  • Xarray ทำให้การใช้ NumPy แบบพื้นฐานลดลงไปมากและทำงานได้มีประสิทธิภาพขึ้นมาก
  • อยากรู้ว่าในเฟรมเวิร์กอย่าง Tensorflow, Keras, Pytorch มีอะไรคล้ายกันไหม เพราะจำได้ว่าเคยดีบักเรื่องที่พูดถึงก่อนหน้านี้อย่างยากลำบาก
  • ขอบคุณสำหรับคำแนะนำและตั้งใจจะลองใช้แน่นอน เคยนึกว่ามีแค่ตัวเองที่รู้สึกอึดอัดกับไวยากรณ์อย่าง array[:, :, None] พอเห็นว่ามีคนคิดเหมือนกันก็ดีใจ
  • ในสาย biosignal มี NeuroPype ที่อยู่บน NumPy และรองรับ named axes สำหรับ n-dimensional tensor พร้อมเก็บข้อมูลต่อ element ในแต่ละแกนได้ด้วย เช่น ชื่อช่องสัญญาณ ตำแหน่ง เป็นต้น
  • ทำให้นึกถึงช่วงที่ NumPy แตกสายมาจากไลบรารี Numeric กับ Numarray ลองจินตนาการว่าฝั่ง Numarray เถียงเรื่องนี้ต่อมา 20 ปี จนได้เงินทุนแล้วเปลี่ยนชื่อเป็น Xarray และในที่สุดชนะ NumPy ได้ (แน่นอนว่าส่วนใหญ่เป็นเรื่องแต่ง)
• เหตุผลหนึ่งที่เริ่มใช้ Julia คือไวยากรณ์ของ NumPy ยากเกินไป พอย้ายจาก MATLAB มา NumPy กลับรู้สึกว่าเขียนโปรแกรมแย่ลง และใช้เวลาไปกับการเรียนรู้เทคนิครีดประสิทธิภาพมากกว่าคณิตศาสตร์ ใน Julia ทั้ง vectorization และ loop ต่างก็ทำงานได้ดี จึงโฟกัสกับความอ่านง่ายของโค้ดได้อย่างเดียว รู้สึกได้ถึงประสบการณ์และอารมณ์นี้จากบทความเต็ม ๆ มองว่าแนวทางแบบ “กล่องดำ” ที่คิดว่าอะไรอย่าง np.linalg.solve เร็วที่สุดเสมอแล้วบังคับให้ทุกอย่างต้องเข้ารูปนั้นไม่ถูกต้อง เพราะยังมีอีกหลายเหตุผลที่การเขียน problem-specific kernel เองจะดีกว่า
  • สาเหตุมาจาก Julia เป็นภาษาที่ออกแบบมาสำหรับ scientific computing โดยตรง ส่วน NumPy เป็นไลบรารีที่พยายามยัดลงบนภาษาที่ไม่ได้ออกแบบมาเพื่อสิ่งนี้ หวังว่าสักวัน Julia จะชนะ และคนที่ใช้ Python เพราะแรงเฉื่อยจาก network effect จะได้เป็นอิสระ
  • MATLAB เองถ้าไม่ vectorize แล้วใช้ loop ก็ช้าไม่ต่างจาก Python ปัญหาใหญ่สุดคือความช้าของ Python และแม้ Julia จะมีข้อดีชัดเจน แต่ในทางปฏิบัติก็ยังใช้ได้ในงานที่ค่อนข้างจำกัด Python เองก็มี JIT hack โผล่มาบ้างแต่ก็ยังไม่สมบูรณ์ จึงต้องการทางเลือกแทน Python อย่างมาก
  • MATLAB ต่างจริงหรือ? loop ก็ยังช้าเหมือนเดิม และสิ่งที่เร็วที่สุดก็คือกล่องดำที่ optimize มาเต็มที่อย่างตัวดำเนินการ \\
  • Fortran รุ่นใหม่ก็เหมือน Julia ตรงที่ทั้ง vectorization และ loop ทำงานได้เร็ว จึงสนใจแต่ความอ่านง่ายได้เลย
• ถ้าสรุปความไม่พอใจต่อ numpy เมื่อเทียบกับ Matlab และ Julia ก็คือ แต่ละฟังก์ชันมีอาร์กิวเมนต์เกี่ยวกับแกน ชื่อเรียก และวิธีรองรับ vectorization ไม่เหมือนกันเลย และถ้าจะใช้ฟังก์ชันกับอีกแกนหนึ่งก็มักต้องเขียนโค้ดใหม่แทบทั้งหมด ทั้งที่พื้นฐานของการเขียนโปรแกรมคือ abstraction แต่ NumPy กลับทำให้เรื่องนี้ยาก ใน Matlab โค้ดแบบ vectorized มักใช้ต่อได้เกือบตรง ๆ หรืออย่างน้อยก็แก้แบบชัดเจน แต่ใน NumPy ต้องคอยเปิดเอกสาร ดู transpose/reshape และการจัดรูปชนิดข้อมูลก็ไม่สม่ำเสมอจนคลุมเครือ
  • การรองรับอาร์เรย์ 3 มิติขึ้นไปของ Matlab อ่อนมาก จึงกลับทำให้ปัญหาแบบที่บทความพูดถึงเกิดได้ยาก
  • ปัญหาข้อที่สองอาจลองใช้ jax กับ vmap ดูได้
  • ถ้าอยากเขียนฟังก์ชันสำหรับอาร์เรย์ 2x2 แล้วนำไปใช้กับบางส่วนของอาร์เรย์ 3x2x2 ก็ทำได้ด้วย slice และ squeeze เป็นต้น ตัวปัญหาเองคลุมเครือจนแทบไม่เข้าใจว่ากังวลอะไร
  • จัดการด้วย reshape ได้
• สิ่งที่สับสนที่สุดใน numpy คือไม่ชัดเลยว่าการดำเนินการไหนจะทำงานแบบ vectorized และก็ไม่มีไวยากรณ์ dot แบบ Julia ให้ระบุอย่างชัดเจน อีกทั้งยังมีหลุมพรางเรื่องชนิดค่าที่คืนกลับมาอีกมาก เช่น ถ้ามีอ็อบเจ็กต์ poly1d ชื่อ P แล้วคูณจากด้านขวาด้วย z0 จะได้ poly1d กลับมา แต่ถ้าคูณจากด้านซ้ายเป็น z0*P จะได้เพียงอาร์เรย์ ทำให้เกิดการแปลงชนิดแบบเงียบ ๆ นอกจากนี้ค่าสัมประสิทธิ์นำของ quadratic ก็เข้าถึงได้ทั้ง P.coef[0] และ P[2] ซึ่งชวนให้สับสน ทางการบอกว่า poly1d เป็น API แบบ “เก่า” และแนะนำให้โค้ดใหม่ใช้คลาส Polynomial แต่ในความเป็นจริงก็ไม่มีแม้แต่คำเตือน deprecated กับดักอย่างการแปลงชนิดและความไม่สอดคล้องของ datatype แบบนี้มีอยู่ทั่วทั้งไลบรารีจนกลายเป็นฝันร้ายเวลา debug
• เห็นด้วยกับประเด็นที่ผู้เขียนชี้ไว้ ตอนย้ายจาก Matlab ไป Numpy ก็เจอความอึดอัดหลายอย่าง และรู้สึกว่าการ slice ข้อมูลของ Numpy ก็ใช้งานลำบากกว่า Matlab/Julia ด้วย แต่ถ้าคิดถึงค่าไลเซนส์ของ toolbox ใน Matlab แล้ว ข้อเสียของ Numpy ก็พอรับได้ ปัญหาในบทความส่วนใหญ่เกิดกับ tensor ที่เกิน 2 มิติ และ Numpy เองก็มีรากฐานจากเมทริกซ์ (2D) จึงเป็นข้อจำกัดที่พอเข้าใจได้ ไลบรารีเฉพาะทางอย่าง Torch อาจดีกว่าแต่ก็ไม่ได้ง่ายนัก สุดท้ายความจริงคงเป็นว่า "NumPy แย่กว่าภาษา/ไลบรารีอาร์เรย์อื่น ๆ อยู่นิดหน่อย แต่ก็แทบไม่มีอย่างอื่นให้ใช้"
  • NumPy ตั้งเป้าหมายเป็นอาร์เรย์ n มิติมาตั้งแต่แรกในฐานะภาคต่อของ numarray ดังนั้นไม่ได้หยุดอยู่แค่ 2D
• ปัญหาใหญ่ที่สุดของ ecosystem data science บน Python คือทุกอย่างไม่เป็นมาตรฐาน ไลบรารีเป็นสิบตัวทำงานต่างกันราวกับเป็นคนละ 4 ภาษา และที่พอเหมือนกันก็มีแค่ประมาณ to_numpy() สุดท้ายจึงเสียเวลากับการแปลง format ข้อมูลมากกว่าการแก้ปัญหาจริง Julia เองก็ไม่ได้มีแต่ข้อดี แต่การเชื่อมต่อกันระหว่างไลบรารีหลายตัว เช่น เรื่องหน่วยและความไม่แน่นอน ทำได้ดีกว่า ส่วน Python มักต้องเขียน boilerplate เยอะเสมอ
  • โครงการ array-api กำลังพยายามทำ API มาตรฐานสำหรับการจัดการอาร์เรย์ให้ครอบคลุม ecosystem Python ทั้งหมด
  • ฝั่ง R กลับซับซ้อนยิ่งกว่าเพราะมีระบบคลาสถึง 4 แบบ
• สงสัยว่าทำไมผู้คนถึงใช้ numpy แทน sage
• บางปัญหาแก้ได้ด้วยการใช้ numpysane กับ gnuplotlib หลังจากมีชุดนี้ก็เริ่มใช้ numpy กับทุกงานอย่างจริงจัง ถ้าไม่มีคงแทบใช้ต่อไม่ไหว
  • numpysane สุดท้ายก็ยังเป็น Python loop ไม่ใช่ vectorization จริง
  • ขอบคุณสำหรับคำแนะนำ ปกติก็บ่นเรื่องนี้อยู่บ่อย ๆ เลยไม่เคยนึกว่าจะมีไลบรารีชั้นบนที่เรียบง่ายแบบนี้
• เคยลองเอา matrix multiplication ทั้งหมดไปใส่ใน einsum เพื่อทำ vectorized multi-head attention และใช้ optimize="optimal" ให้เลือกอัลกอริทึม matrix chain multiplication เพื่อเร่งความเร็ว ผลคือเร็วขึ้นจาก implementation แบบ vectorized ทั่วไปราว 2 เท่าจริง แต่ที่น่าประหลาดคือ implementation แบบตรงไปตรงมาที่ใช้ loop กลับเร็วกว่า ใครอยากรู้เหตุผลลองดูโค้ดได้ คาดว่ายังมีพื้นที่ให้ปรับปรุง cache coherency ภายใน einsum ได้อีก