1 คะแนน โดย GN⁺ 2024-07-02 | 1 ความคิดเห็น | แชร์ทาง WhatsApp
  • เมื่อนำพหุนามดีกรีสูงมาขยายด้วยวิธีแบบมัธยมปลาย ต้องคูณทุกคู่ของพจน์ ทำให้ต้นทุน O(n²) กลายเป็นคอขวดอย่างรวดเร็ว
  • การคูณเวกเตอร์สัมประสิทธิ์ของพหุนามเหมือนกับ คอนโวลูชัน ของสัญญาณไม่ต่อเนื่อง และผลลัพธ์ของ [2, 3, 4] กับ [5, 6, 7] คือ [10, 27, 52, 45, 28]
  • DFT ย้ายสัญญาณไม่ต่อเนื่องไปยังโดเมนความถี่ ส่วน FFT คำนวณการแปลงเดียวกันนี้ได้ใน O(n log n) ทำให้เกิดความแตกต่างเมื่ออินพุตมีขนาดใหญ่
  • คอนโวลูชันในโดเมนเวลาจะเปลี่ยนเป็นการคูณแบบรายองค์ประกอบในโดเมนความถี่ ดังนั้นหากแปลงด้วย FFT แล้วคูณ จากนั้นแปลงกลับด้วย IFFT ก็จะจัดการการคูณพหุนามได้เร็วขึ้น
  • สำหรับดีกรีเล็ก ๆ ต้นทุนการไป-กลับของ FFT/IFFT อาจหักล้างประโยชน์ที่ได้ แต่เมื่อดีกรีใหญ่ขึ้น วิธี FFT จะมีประสิทธิภาพมากกว่า

ทำไมการคูณพหุนามจึงช้าลง

  • พหุนาม P(x) แสดงได้ในรูปผลรวมของสัมประสิทธิ์ a_k กับพจน์ยกกำลังของตัวแปร x
    • ตัวอย่าง P(x)=5x²+2x+9 เป็นพหุนามดีกรี 2
    • เวกเตอร์สัมประสิทธิ์อาจแสดงเป็น [5, 2, 9] หรือ [9, 2, 5] ได้ ขึ้นอยู่กับรูปแบบการเขียน
  • การบวกและการลบค่อนข้างตรงไปตรงมา เพราะแค่บวกหรือลบพจน์ที่มีดีกรีเดียวกัน
    • ใน Python สามารถใช้ zip(p, q) เพื่อวนผ่านสัมประสิทธิ์แต่ละตัวแล้วคำนวณ a + b หรือ a - b
    • หากดีกรีต่างกัน สามารถใช้ zip_longest ได้
  • การคูณต้องคูณแต่ละพจน์เข้าด้วยกัน แล้วรวมพจน์ที่มีดีกรีเดียวกันอีกครั้ง ทำให้ปริมาณการคำนวณเพิ่มขึ้น
    • ผลลัพธ์ของ (2x²+3x+4) × (5x²+6x+7) คือ 10x⁴+27x³+52x²+45x+28
    • ความซับซ้อนของวิธีนี้คือ O(n²) และยิ่งดีกรีสูง จำนวนการคูณที่ต้องใช้ก็ยิ่งเพิ่มขึ้น

เวกเตอร์สัมประสิทธิ์และคอนโวลูชัน

  • ในโดเมนไม่ต่อเนื่อง คอนโวลูชันของสัญญาณสองตัว p และ q นิยามเป็น y[n]=Σ p[k]·q[n-k]
  • การคำนวณทำโดยกลับลำดับ q แล้วเลื่อนจากซ้ายไปขวาบน p พร้อมบวกผลคูณขององค์ประกอบที่ซ้อนทับกัน
  • สัญญาณตัวอย่างมีดังนี้
    • p = [2, 3, 4]
    • q = [5, 6, 7]
  • เมื่อกลับลำดับ q แล้วเลื่อน สัมประสิทธิ์เอาต์พุตแต่ละตัวจะถูกสร้างตามลำดับดังนี้
    • 2×5 = 10
    • 2×6 + 3×5 = 27
    • 2×7 + 3×6 + 4×5 = 52
    • 3×7 + 4×6 = 45
    • 4×7 = 28
  • ผลลัพธ์ของคอนโวลูชันคือ y = [10, 27, 52, 45, 28]
    • ซึ่งเหมือนกับสัมประสิทธิ์ของ 10x⁴+27x³+52x²+45x+28 ที่ได้จากการคูณพหุนาม
    • ดังนั้น การคูณพหุนาม จึงมองได้ว่าเป็นคอนโวลูชันของเวกเตอร์สัมประสิทธิ์

Fourier Transform และ FFT

  • Fourier Transform แปลงสัญญาณจาก โดเมนเวลา ไปยัง โดเมนความถี่
    • ในมุมมองเวลา เรามองสัญญาณเป็นค่าของแต่ละช่วงเวลา
    • ในมุมมองความถี่ เราตีความสัญญาณเป็นผลรวมของความถี่การสั่นที่แตกต่างกัน
  • ความถี่การสั่นแสดงด้วยไซน์และโคไซน์ โดยแต่ละตัวมีสัมประสิทธิ์และเฟส
  • เมื่อนำ FFT ไปใช้กับคลื่นไซน์บริสุทธิ์ 5Hz ในโดเมนความถี่จะปรากฏเหมือนเดลตาที่ตำแหน่ง 5Hz
    • สิ่งนี้แสดงให้เห็นว่าคลื่นไซน์ในโดเมนเวลาสามารถแทนได้ด้วยไซน์ 5Hz เพียงตัวเดียว
  • คำศัพท์ที่เกี่ยวข้องแยกได้ดังนี้
    • Fourier Transform(FT): การแปลงฟูเรียร์ที่นิยามบนโดเมนต่อเนื่อง
    • Discrete Fourier Transform(DFT): การแปลงฟูเรียร์ที่นิยามสำหรับสัญญาณไม่ต่อเนื่อง
    • Fast Fourier Transform(FFT): อัลกอริทึมที่คำนวณ DFT ใน O(n log n) แทน O(n²)
  • DFT แปลงสัญญาณเวลาไม่ต่อเนื่อง x[n] เป็น X[k] ในโดเมนความถี่
    • X[k] แต่ละตัวคำนวณโดยนำตัวอย่างอินพุตไปคูณกับจำนวนเชิงซ้อนที่แทนความถี่เฉพาะ แล้วนำมาบวกกัน

เปลี่ยนเป็นการคูณในโดเมนความถี่

  • ข้อดีหลักของ DFT และโดเมนความถี่คือ สามารถเปลี่ยนคอนโวลูชันให้เป็นการคูณแบบรายองค์ประกอบได้
    • การคอนโวลูชันสัญญาณสองตัวในโดเมนเวลา เทียบเท่ากับการคูณสัญญาณสองตัวในโดเมนความถี่
    • การคูณสามารถคำนวณได้เร็วกว่าคอนโวลูชัน
  • ขั้นตอนในการคูณพหุนามอย่างรวดเร็วมีดังนี้
    • แปลงพหุนามไปยังโดเมนความถี่ด้วย FFT: O(n log n)
    • คูณแบบรายองค์ประกอบในโดเมนความถี่: O(n)
    • แปลงผลลัพธ์กลับมายังโดเมนเวลาด้วย IFFT: O(n log n)
  • โดยรวมแล้ว การใช้ FFT ทำให้คูณพหุนามได้ด้วยความซับซ้อน O(n log n)
  • สำหรับพหุนามขนาดใหญ่ จะเร็วกว่าการคูณแบบมัธยมปลาย O(n²)

การติดตั้งใช้งานด้วย Python และเบนช์มาร์ก

  • multiply_naive ใช้ลูปซ้อนเพื่อคูณทุกคู่ของสัมประสิทธิ์ แล้วบวกเข้าที่ตำแหน่งผลลัพธ์ i + j
    • ความยาวของผลลัพธ์คือ len(p) + len(q) - 1
    • ความซับซ้อนคือ O(n²)
  • multiply_fft ทำการคูณสัมประสิทธิ์บนพื้นฐาน FFT/IFFT
    • คำนวณความยาวที่เป็นกำลังของ 2 ซึ่งไม่น้อยกว่า len(p) + len(q) - 1 เพื่อให้รองรับความยาวผลลัพธ์ได้
    • แพดอินพุตทั้งสองด้วย np.pad
    • คูณค่าที่แปลงด้วย np.fft.fft แบบรายองค์ประกอบ
    • แปลงกลับด้วย np.fft.ifft แล้วปัดเศษส่วนจริงเพื่อแปลงเป็นสัมประสิทธิ์จำนวนเต็ม
  • สำหรับอินพุตตัวอย่าง p = [2, 3, 4], q = [5, 6, 7] ทั้งสองวิธีคืนค่า [10, 27, 52, 45, 28]
  • ในเบนช์มาร์ก เปรียบเทียบวิธี FFT กับ multiply_convolve ซึ่งใช้ np.convolve แทน multiply_naive
    • เพราะ multiply_naive ใช้ลูป Python ที่ช้า จึงเปรียบเทียบโดยตรงกับวิธี FFT ที่ใช้ np.fft.fft ได้ยาก
    • np.convolve ทำการคำนวณเดียวกันด้วยโค้ด C ระดับต่ำ
  • เพิ่มดีกรีในช่วง range(1, 30000, 1000) และสร้างพหุนามสองตัวในแต่ละดีกรีด้วยสัมประสิทธิ์สุ่มระหว่าง 1 ถึง 999999
    • แต่ละวิธีวัดเวลาเฉลี่ยด้วย n_runs = 5
    • ที่ดีกรีต่ำ วิธี FFT อาจไม่ได้เปรียบเพราะต้นทุนการแปลงไป-กลับของ FFT/IFFT
    • เมื่อดีกรีเพิ่มขึ้น วิธี FFT จะแสดงผลลัพธ์ที่มีประสิทธิภาพกว่ามาก

1 ความคิดเห็น

 
GN⁺ 2024-07-02
ความคิดเห็นจาก Hacker News
  • สิ่งที่น่าขัดใจเสมอในคำอธิบายแบบนี้คือมักลืมเรื่อง ข้อผิดพลาดเชิงตัวเลข
    คุณไม่สามารถนามธรรมการคูณสัมประสิทธิ์ให้เป็นแค่ “เวลาคงที่” ได้ ถ้าจะทำแบบนั้น จะนามธรรมการคูณทั้งหมดไปเลยก็ไม่ต่างกัน เมื่อคำนึงถึงความแม่นยำเชิงตัวเลขแล้ว มันใกล้กับ O(n (log n)^3) มากกว่า [1]
    [1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....

    • ขอบเขตความคลาดเคลื่อน ที่อ้างในบทความนั้นมองในแง่ร้ายเกินไป ใน Knuth ฉบับล่าสุดมีขอบเขตที่ถูกต้องแล้ว เพราะฉันเป็นคนบอกเขาเอง
    • คงจะดีถ้าสามารถใช้ การคำนวณบนฐานควอเทอร์เนียน ที่กล่าวถึงในโพสต์ต้นฉบับเพื่อลดหรือแม้แต่กำจัดความผิดพลาดจากการคูณได้ [1],[2],[3]
      [1] One-Dimensional Quaternion Discrete Fourier Transform and an Approach to Its Fast Computation:
      https://www.mdpi.com/2079-9292/12/24/4974
      [2] Convolution Theorems for Quaternion Fourier Transform: Properties and Applications:
      https://onlinelibrary.wiley.com/doi/10.1155/2013/162769
      [3] On the Matrix Form of the Quaternion Fourier Transform and Quaternion Convolution:
      https://arxiv.org/abs/2307.01836
    • ถ้าสัมประสิทธิ์เป็นจำนวนเต็ม ก็สามารถใช้ NTT กับมอดูลัสที่ใหญ่พอเพื่อให้ได้ผลลัพธ์ที่แม่นยำ และโดยเฉพาะในฮาร์ดแวร์ เวลาคูณอาจเร็วขึ้นด้วย
    • นี่จึงเป็นเหตุผลที่เราแยก วิทยาการคอมพิวเตอร์ ออกจาก วิศวกรรมซอฟต์แวร์ :)
  • วิธีนี้ใช้คูณเลขยาว ๆ เข้าด้วยกันได้ ประเด็นสำคัญคือ การคูณพหุนาม เทียบเท่ากับการคูณเลขยาวแบบทั่วไปที่ยังไม่ทำการทด (carry)
    ตัวอย่างเช่น ถ้ามีตัวเลข 1000 หลัก ก็เอาแต่ละหลักมาเป็นสัมประสิทธิ์ของพหุนามที่มี 1000 พจน์ จากนั้นก็ใช้วิธี FFT ที่บทความอธิบายไว้เพื่อคูณพหุนามเหล่านี้ได้ เมื่อต้องแปลงผลลัพธ์กลับเป็นตัวเลข ก็ต้องจัดการเรื่องการทด ถ้าพจน์ใดมีค่ามากกว่า 10 ก็ส่งส่วนที่เกินไปยังหลักถัดไป แล้วแปลงสัมประสิทธิ์กลับเป็นตัวเลข
    นี่คือแนวคิดพื้นฐาน แต่จะมีรายละเอียดที่ละเอียดอ่อนในเรื่องความแม่นยำที่ต้องใช้สำหรับการทด และการรับประกันว่าการปัดผล FFT ไปเป็นจำนวนเต็มที่ใกล้ที่สุดยังคงถูกต้อง วิธีนี้คือแนวทางที่ GMP ซึ่งเป็นไลบรารีหลักของวงการ ใช้ในการคูณจำนวนขนาดใหญ่

    • ตามที่บอก เลขฐาน 10 สามารถแทนเป็นพหุนามที่ x=10 ได้ จึงเข้าใจได้ เช่น 983 = 9x^2 + 8x + 3 หรือ [9, 8, 3]
      สิ่งที่อยากรู้คือในทางปฏิบัติตัวเลขต้องใหญ่แค่ไหนถึงจะคุ้ม และมันถูกนำไปใช้ที่ไหน
  • ถ้ายังไม่เคยดู วิดีโอนี้น่าดูมาก
    https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
    มันอธิบายการไล่อนุมานอัลกอริทึม FFT จาก การคูณพหุนาม ได้ยอดเยี่ยมมาก ฉันกลับมาดูใหม่ประมาณทุก 6 เดือน

  • คุณสมบัติ “คอนโวลูชันคือการคูณแบบจุดต่อจุด” ของ FFT ใช้ได้กับ กลุ่มการคูณเชิงวัฏจักร ใด ๆ ด้วย สำหรับการไล่อนุมานเชิงพีชคณิตมากขึ้น ดู https://www.sciencedirect.com/science/article/pii/S002200007...
    บางครั้งสิ่งนี้ถูกเรียกว่า “harmonic FFT” และยังมี non-harmonic FFT ด้วย: [LCH14] “additive NTT” บน GF(2^n), [HLP24] circle FFT บนวงกลมหน่วย X^2+Y^2=1 ของฟิลด์จำกัด, และ [BCKL21] ecfft บนสายไอโซจีนีของเส้นโค้งวงรี
    [LCH14]: https://arxiv.org/abs/1404.3458
    [HLP24]: https://eprint.iacr.org/2024/278
    [BCKL21]: https://arxiv.org/pdf/2107.08473

  • ใครกันแน่ที่เป็นคนแรกที่เสนอให้ใช้ FFT เพื่อการคูณพหุนามที่เร็วขึ้น?
    พอดีช่วงนี้สงสัยเลยไปค้นดู แม้จะไล่ตามการอ้างอิงได้ไม่เก่งนัก แต่ก็ย้อนกลับไปได้ถึงบทความปี 1995 ของ David Eppstein [0] ตรงนี้เขาใช้มันเพื่อแก้ปัญหา subset sum แบบอัปเดตเพิ่มทีละน้อยได้อย่างมีประสิทธิภาพ
    รู้สึกว่าใน TAOCP ของ Knuth น่าจะมีอะไรที่เก่ากว่านี้แน่
    ที่น่าตกใจไม่น้อยคือ FFT polynomial multiplication สามารถใช้แก้ปัญหา subset sum แบบ exact ที่อนุญาตให้ทำซ้ำได้ในเวลาแบบกึ่งเอ็กซ์โปเนนเชียล [1] ประเด็นสำคัญคืออัลกอริทึมนี้เป็น O(N log N) โดยที่ N ไม่ใช่ขนาดของเซต แต่เป็นค่าสูงสุดของสมาชิก ดังนั้นจึงไม่ใช่อะไรทำนอง counterexample ต่อ P ≠ NP
    [0] https://escholarship.org/content/qt6sd695gn/qt6sd695gn.pdf
    [1] https://x.com/festivitymn/status/1788362552998580473?s=46&t=...

    • ดูเหมือนว่า Pollard [1], Nicholson [2], และ Schönhage-Strassen [3] จะคิดเรื่องนี้ขึ้นมาอย่างอิสระจากกันในช่วงเวลาใกล้เคียงกันแต่ใช้คนละแนวทาง
      มีการกล่าวกันว่า Strassen ค้นพบแนวทางแบบ Pollard ตั้งแต่ปี 1968 แต่ไม่มีหลักฐานเอกสารรองรับ และแม้จะไม่ใช่จุดกำเนิดของ FFT เองก็ตาม ก็ควรคำนึงด้วยว่าบทความปี 1965 ของ Cooley-Tukey [4] เป็นตัวจุดกระแสการวิจัย FFT และการประยุกต์ใช้อย่างจริงจัง เรื่องนี้เกิดขึ้นหลังจากนั้นไม่กี่ปี
      [1] https://doi.org/10.1090/S0025-5718-1971-0301966-0
      [2] https://doi.org/10.1016/S0022-0000(71)80014-4
      [3] https://doi.org/10.1007/BF02242355
      [4] https://doi.org/10.1090/S0025-5718-1965-0178586-1
    • ที่เก่าสุดอาจเป็น Gentleman กับ Sande ในปี 1966 และชื่อก็เท่มากสำหรับมาตรฐานปี 1966: “Fast Fourier Transforms: for fun and profit”
      https://www.cis.rit.edu/class/simg716/FFT_Fun_Profit.pdf
    • อัลกอริทึม Schönhage–Strassen ในปี 1971 แท้จริงแล้วก็คือการคูณพหุนามโดยใช้ FFT นั่นเอง: https://en.m.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Stras...
  • ผมคิดว่าแมชชีนเลิร์นนิงทั้งหมดก็คือการแก้สมการคอนโวลูชัน
    บทความนี้พูดถึงในบริบทของ reinforcement learning https://arxiv.org/abs/1712.06115 แต่แนวทางส่วนใหญ่ก็เข้ากับกระบวนทัศน์นี้

    • โดยพื้นฐานแล้วก็หมายถึงkernel method ใช่ไหม?
  • เมื่อกี้เพิ่ง implement อัลกอริทึม (matrix profile) ที่ใช้ FFT เพื่อคำนวณinner product สำหรับเซตขนาดใหญ่ของ subsequence ใน time series โดยความยาวของ time series n อาจไปถึงระดับหลายร้อยล้านได้
    การคำนวณคอนโวลูชันแบบเร็วด้วย FFT ทำให้เวลาในการคำนวณลดจาก O(n) เหลือ O(log n) และในสเกลนี้ความเร็วที่เพิ่มขึ้นถือว่ามหาศาล ยิ่งใช้ GPU ด้วยก็เร็วขึ้นไปอีก เช่น ประมวลผลข้อมูล 10 ล้านจุดบนโน้ตบุ๊กได้ใน 0.1 วินาที

  • “ลูกเล่น” สำคัญของการดำเนินการนี้ดูเหมือนจะเป็นความเข้าใจข้อนี้:

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

    • ไม่แน่ใจว่าจะช่วยได้ไหม แต่การคูณพหุนามสองตัวที่เรียนกันในโรงเรียน จริง ๆ แล้วก็คือคอนโวลูชัน
      มันมีสมบัติว่า “การทำคอนโวลูชันของสัญญาณสองตัวในโดเมนเวลา เทียบเท่ากับการคูณสัญญาณสองตัวนั้นในโดเมนความถี่” และ FFT ทำให้เราแปลงจากโดเมนเวลาไปเป็นโดเมนความถี่ได้ ดังนั้นจึงย้ายพหุนามไปยังโดเมนความถี่ด้วย FFT แล้วในโดเมนนั้นก็แค่คูณกัน ซึ่งเร็วกว่า convolution หวังว่าขั้นที่ขาดไปจะชัดขึ้นจากตรงนี้ และถ้ายังมีส่วนที่ขาดอยู่ก็อาจอัปเดตบทความได้
  • ถ้าอย่างนั้น การแยกตัวประกอบจำนวนเต็ม ถือเป็น deconvolution แบบไม่ต่อเนื่องหรือเปล่า? ถ้าวางการแทนแบบ FFT หรือก็คือการคูณแบบ pointwise คู่กับ tableax หรือก็คือการคูณยาว/การบวกแบบทดทั่วไป แล้วความสมมาตรเกิดแตกขึ้นมา ผมสงสัยว่าจะได้ข้อมูลเพียงพอสำหรับอัลกอริทึมที่เร็วขึ้นไหม

  • แน่นอนว่าการคูณพหุนามแบบตรงไปตรงมานั้นช้าเมื่อเทียบกับดีกรีของพหุนาม แต่ในทางปฏิบัติจะมีตอนไหนที่ต้องจัดการกับพหุนามดีกรี 100สองตัวจริง ๆ กัน?
    ด้วยเหตุนี้จึงมักมีภาพจำว่าระบบคอมพิวเตอร์พีชคณิตไม่ได้ใช้วิธีแบบนี้

    • ระบบคอมพิวเตอร์พีชคณิต เช่น chebfun ของ Matlab จะแปลงฟังก์ชันตามอำเภอใจให้เป็นพหุนามดีกรีมากกว่า 100เพื่อให้หาค่า root, ค่าที่เหมาะที่สุด ฯลฯ ได้ง่ายขึ้น
    • พบได้บ่อยมากในงานการแก้ไขข้อผิดพลาดและการประมวลผลสัญญาณ
      https://www.youtube.com/watch?v=CcZf_7Fb4Us
      https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor... เป็นตัวอย่างหนึ่ง
    • ผมเคยสร้างโปรแกรม[1] เพราะอยากทำ reverse engineer กับพารามิเตอร์ CRC checksumของไฟล์ขนาดใหญ่ โดยแปลงไฟล์ให้เป็นพหุนาม GF(2) ดีกรีระดับหลายล้าน แล้วคำนวณตัวหารร่วมมาก หากไม่มีการคูณแบบอาศัย FFT ก็แทบเป็นไปไม่ได้เลยที่จะทำให้เสร็จภายในเวลาที่สมเหตุสมผล
      [1]: https://github.com/8051enthusiast/delsum
    • มุมมองแบบคอนโวลูชันเช่นนี้ และ GPU kernel แบบเร็วสำหรับ FFT ถูกใช้ในโมเดล state space บางตัวก่อนยุค Mamba สำหรับการทำ long-sequence modeling โดยในที่นี้พหุนามก็คือลำดับอินพุต
      มีข้อมูลเกี่ยวกับแนวทางนี้อยู่มากในบล็อกโพสต์ของ Hazy Research ช่วงปี 2020~2023
    • ดู https://news.ycombinator.com/item?id=40306339
      "(...) เคยจัดการกับนิพจน์ที่มีขนาดเกือบ 1 เทราไบต์ และมีพจน์มากกว่า 100 ล้านพจน์ในการวิจัยฟิสิกส์"