Convolution, Fast Fourier Transform और Polynomials (2022)

 (alvarorevuelta.com)
1 पॉइंट द्वारा GN⁺ 2024-07-02 | 1 टिप्पणियां | WhatsApp पर शेयर करें
  • बड़े degree वाले polynomials को स्कूल-स्तर की विधि से expand करने पर हर term pair को multiply करना पड़ता है, इसलिए O(n²) लागत जल्दी bottleneck बन जाती है
  • polynomial coefficient vectors का multiplication, discrete signals के convolution के बराबर होता है, और [2, 3, 4] तथा [5, 6, 7] का परिणाम [10, 27, 52, 45, 28] होता है
  • DFT discrete signal को frequency domain में ले जाता है, और FFT उसी transform को O(n log n) में compute करता है, जिससे बड़े inputs पर बड़ा अंतर पड़ता है
  • time domain का convolution, frequency domain में element-wise multiplication में बदल जाता है, इसलिए FFT से transform करके, multiply करके, फिर IFFT से वापस लाने पर polynomial multiplication को तेज़ किया जा सकता है
  • छोटे degree पर FFT/IFFT round-trip की लागत लाभ को कम कर सकती है, लेकिन degree बढ़ने पर FFT तरीका अधिक efficient हो जाता है

polynomial multiplication धीमा क्यों हो जाता है

  • polynomial P(x) को coefficients a_k और variable x की powers वाले terms के योग के रूप में व्यक्त किया जाता है
    • उदाहरण P(x)=5x²+2x+9 degree 2 का polynomial है
    • coefficient vector को notation के अनुसार [5, 2, 9] या [9, 2, 5] की तरह लिखा जा सकता है
  • addition और subtraction अपेक्षाकृत सरल हैं, क्योंकि समान degree वाले terms को जोड़ना या घटाना होता है
    • Python में zip(p, q) से हर coefficient पर iterate करके a + b या a - b निकाला जा सकता है
    • degree अलग होने पर zip_longest का उपयोग किया जा सकता है
  • multiplication में हर term को दूसरे हर term से multiply करना पड़ता है, फिर समान degree वाले terms को दोबारा जोड़ना पड़ता है, इसलिए computation बढ़ जाता है
    • (2x²+3x+4) × (5x²+6x+7) का परिणाम 10x⁴+27x³+52x²+45x+28 है
    • इस विधि की complexity O(n²) है, और degree बढ़ने पर required multiplications की संख्या बढ़ती जाती है

coefficient vectors और convolution

  • discrete domain में दो signals p और q का convolution y[n]=Σ p[k]·q[n-k] से define किया जाता है
  • computation में q को reverse करके p के ऊपर बाएँ से दाएँ slide किया जाता है, और overlap होने वाले elements के गुणनफल जोड़े जाते हैं
  • उदाहरण signals इस प्रकार हैं
    • p = [2, 3, 4]
    • q = [5, 6, 7]
  • q को reverse करके slide करने पर हर output coefficient इस क्रम में बनता है
    • 2×5 = 10
    • 2×6 + 3×5 = 27
    • 2×7 + 3×6 + 4×5 = 52
    • 3×7 + 4×6 = 45
    • 4×7 = 28
  • convolution का परिणाम y = [10, 27, 52, 45, 28] है
    • यह polynomial multiplication से मिले 10x⁴+27x³+52x²+45x+28 के coefficients के बराबर है
    • इसलिए polynomial multiplication को coefficient vectors के convolution के रूप में देखा जा सकता है

Fourier transform और FFT

  • Fourier transform signal को time domain से frequency domain में बदलता है
    • time के नज़रिए से signal को किसी खास समय पर उसके value के रूप में देखा जाता है
    • frequency के नज़रिए से signal को अलग-अलग oscillation frequencies के योग के रूप में समझा जाता है
  • oscillation frequencies को sine और cosine से व्यक्त किया जाता है, और इनके अपने coefficient तथा phase होते हैं
  • 5Hz की pure sine wave पर FFT लगाने से frequency domain में 5Hz की position पर delta जैसी आकृति दिखाई देती है
    • यह दिखाता है कि time domain की sine wave को एक 5Hz sine से व्यक्त किया जा सकता है
  • संबंधित terms को इस तरह अलग किया जाता है
    • Fourier Transform(FT): continuous domain में define किया गया Fourier transform
    • Discrete Fourier Transform(DFT): discrete signals के लिए define किया गया Fourier transform
    • Fast Fourier Transform(FFT): DFT को O(n²) की जगह O(n log n) में compute करने वाला algorithm
  • DFT discrete-time signal x[n] को frequency domain के X[k] में बदलता है
    • हर X[k] को input samples को किसी खास frequency को दर्शाने वाली complex numbers से multiply करके और फिर जोड़कर निकाला जाता है

frequency domain में multiplication में बदलना

  • DFT और frequency domain का मुख्य लाभ यह है कि convolution को element-wise multiplication में बदला जा सकता है
    • time domain में दो signals का convolution, frequency domain में उन signals के multiplication के बराबर होता है
    • multiplication को convolution की तुलना में अधिक तेज़ी से compute किया जा सकता है
  • polynomial multiplication को तेज़ी से करने की प्रक्रिया इस प्रकार है
    • polynomial को FFT से frequency domain में transform करना: O(n log n)
    • frequency domain में element-wise multiply करना: O(n)
    • परिणाम को IFFT से फिर time domain में transform करना: O(n log n)
  • कुल मिलाकर FFT का उपयोग करके polynomial multiplication को O(n log n) complexity में किया जा सकता है
  • बड़े polynomials के लिए यह स्कूल-स्तर के O(n²) multiplication से तेज़ होता है

Python implementation और benchmark

  • multiply_naive double loop से सभी coefficient pairs को multiply करके परिणाम की position i + j में जोड़ता है
    • परिणाम की length len(p) + len(q) - 1 होती है
    • complexity O(n²) है
  • multiply_fft FFT/IFFT आधारित coefficient multiplication करता है
    • ऐसी power-of-two length निकाली जाती है जो len(p) + len(q) - 1 को समेट सके
    • np.pad से दोनों inputs को pad किया जाता है
    • np.fft.fft से transform किए गए values को element-wise multiply किया जाता है
    • np.fft.ifft से वापस लाकर real part को round करके integer coefficients में बदला जाता है
  • उदाहरण input p = [2, 3, 4], q = [5, 6, 7] पर दोनों तरीके [10, 27, 52, 45, 28] लौटाते हैं
  • benchmark में multiply_naive की जगह np.convolve का उपयोग करने वाले multiply_convolve और FFT तरीके की तुलना की गई है
    • क्योंकि multiply_naive का Python loop धीमा है, इसलिए np.fft.fft का उपयोग करने वाले FFT तरीके से उसकी सीधी तुलना कठिन है
    • np.convolve वही operation low-level C code में करता है
  • degree को range(1, 30000, 1000) दायरे में बढ़ाया जाता है, और हर degree पर 1 से 999999 के बीच random coefficients वाले दो polynomials बनाए जाते हैं
    • हर तरीके के लिए n_runs = 5 रखकर average time मापा जाता है
    • कम degree पर FFT/IFFT round-trip transform cost की वजह से FFT तरीका फ़ायदेमंद नहीं भी हो सकता
    • degree बढ़ने पर FFT तरीका काफ़ी अधिक efficient परिणाम दिखाता है

संबंधित पढ़ाई

1 टिप्पणियां

 
GN⁺ 2024-07-02
Hacker News की राय

  • इस तरह की व्याख्याओं में जो बात हमेशा खटकती है, वह यह है कि आम तौर पर numerical error को भूल जाते हैं
    coefficient multiplication को बस “constant time” मानकर abstract नहीं किया जा सकता। अगर ऐसा करना है, तो शुरुआत से पूरी multiplication को ही abstract कर देना भी वैसा ही होगा। numerical precision को ध्यान में रखें तो यह O(n (log n)^3) के ज्यादा करीब है [1]
    [1]: http://numbers.computation.free.fr/Constants/Algorithms/fft....

    • उस लेख में quoted error bound जरूरत से ज्यादा pessimistic है। Knuth के नए edition में सही bound शामिल है, क्योंकि मैंने उन्हें बताया था
    • अच्छा होगा अगर OP article में बताए गए quaternion-based operations का उपयोग करके multiplication error को घटाया या पूरी तरह खत्म किया जा सके [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
    • अगर coefficients integers हों, तो पर्याप्त बड़े modulus के साथ NTT इस्तेमाल करके exact result मिल सकता है, और खासकर hardware में multiplication time भी तेज हो सकता है
    • इसीलिए computer science और software engineering में फर्क किया जाता है :)

  • इस तरीके से लंबे numbers को आपस में multiply किया जा सकता है। मुख्य बात यह है कि polynomial multiplication, carry किए बिना की गई सामान्य long-number multiplication जैसी ही होती है
    उदाहरण के लिए अगर 1000-digit number है, तो हर digit को 1000 elements वाले polynomial का coefficient मान लें। फिर लेख में बताए गए FFT method से इन polynomials को multiply किया जा सकता है। result को वापस number में बदलने के लिए carry संभालना होगा। अगर कोई element 10 से बड़ा है, तो extra part को अगले digit में carry कर दें, और coefficients को number में बदल दें
    basic idea यही है, और carry के लिए जरूरी precision तथा FFT result को nearest integer पर round करने पर सही result की guarantee देने वाले हिस्से में कुछ subtle points हैं। इसी तरीके से इस क्षेत्र की प्रमुख library GMP big-number multiplication करती है

    • जैसा आपने कहा, decimal number को x=10 वाले polynomial के रूप में represent किया जा सकता है, इसलिए समझ आता है। उदाहरण के लिए 983 = 9x^2 + 8x + 3, यानी [9, 8, 3]
      असल में यह meaningful होने के लिए numbers कितने बड़े होने चाहिए, और इसका इस्तेमाल कहाँ होता है, यह जानना चाहूँगा

  • अगर अभी तक नहीं देखा है, तो यह video देखना अच्छा रहेगा
    https://youtu.be/h7apO7q16V0?si=bmgUEMTQSqU3flIv
    polynomial multiplication से FFT algorithm को derive करता है, और सचमुच शानदार है। मैं लगभग हर 6 महीने में इसे फिर से देखता हूँ

  • FFT की “convolution is pointwise multiplication” property किसी भी cyclic multiplication group में भी valid होती है। ज्यादा algebraic derivation के लिए https://www.sciencedirect.com/science/article/pii/S002200007... देखें
    इसे कभी-कभी “harmonic FFT” भी कहा जाता है, और non-harmonic FFT भी हैं: GF(2^n) पर [LCH14] “additive NTT”, finite field की unit circle X^2+Y^2=1 पर [HLP24] circle FFT, elliptic-curve isogeny sequence पर [BCKL21] ecfft
    [LCH14]: https://arxiv.org/abs/1404.3458
    [HLP24]: https://eprint.iacr.org/2024/278
    [BCKL21]: https://arxiv.org/pdf/2107.08473

  • तेज़ polynomial multiplication के लिए FFT का इस्तेमाल करने का सुझाव सबसे पहले किसने दिया होगा?
    हाल ही में जिज्ञासा हुई तो खोजा; citations को बहुत अच्छे से trace तो नहीं कर पाया, लेकिन David Eppstein के 1995 के paper [0] तक पीछे जा पाया। इसमें incremental updates के बाद partial sum problem को efficiently solve करने के लिए इसका इस्तेमाल किया गया है। यकीनन Knuth के TAOCP में यह इससे पहले रहा होगा।
    यह बात भी काफी चौंकाने वाली थी कि FFT polynomial multiplication से repetition allowed exact partial sum problem को sub-exponential time में solve किया जा सकता है [1]. अहम बात यह है कि यह algorithm O(N log N) है, जहां N set का size नहीं बल्कि maximum element है, इसलिए यह P ≠ NP के लिए कोई counterexample जैसा नहीं है।
    [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] ने लगभग उसी समय अलग-अलग approaches से इसे स्वतंत्र रूप से सोचा था।
      कहा जाता है कि Strassen ने 1968 में Pollard-style approach खोजी थी, लेकिन इसका कोई documented record नहीं है। साथ ही, भले ही यह FFT का जन्म खुद न हो, Cooley-Tukey का 1965 paper [4] FFT और उसके applications पर research को सचमुच शुरू करने वाला रहा—यह भी देखना चाहिए। यह उसके कुछ साल बाद की बात है।
      [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
    • सबसे शुरुआती शायद 1966 में Gentleman और Sande हो सकते हैं, और 1966 के हिसाब से title भी काफी शानदार है: “Fast Fourier Transforms: for fun and profit”
      https://www.cis.rit.edu/class/simg716/FFT_Fun_Profit.pdf
    • 1971 का Schönhage–Strassen algorithm असल में FFT का उपयोग करके किया गया polynomial multiplication ही है: https://en.m.wikipedia.org/wiki/Sch%C3%B6nhage%E2%80%93Stras...

  • मुझे लगता है सारी machine learning convolution equations solve करने का काम है।
    यह paper reinforcement learning के context में इसे discuss करता है https://arxiv.org/abs/1712.06115, लेकिन ज्यादातर approaches इस paradigm में फिट बैठती हैं।

    • मूल रूप से इसका मतलब kernel methods नहीं है क्या?

  • अभी-अभी मैंने FFT का इस्तेमाल करके time-series subsequences के बड़े set के लिए inner products calculate करने वाला algorithm (matrix profile) implement किया है। time series की length n करोड़ों में जा सकती है।
    FFT से fast convolution calculation के कारण computation time O(n) से O(log n) तक घट जाता है, और इस scale पर speedup जबरदस्त है। GPU का इस्तेमाल करें तो यह और तेज हो जाता है—जैसे laptop पर 1 करोड़ data points को 0.1 seconds में process करना।

  • इस operation की core “trick” यह realization लगती है:

    दूसरे शब्दों में, time domain में दो signals का convolution करना frequency domain में दो signals को multiply करने जैसा ही है।
    complex idea को बहुत छोटे steps में तोड़कर, maths में कमजोर मुझे भी किसी तरह समझ आ जाए—ऐसा करने वाला यह अच्छा article है। लेकिन क्या एक intermediate step छूट गया? या उसे reader के लिए ढूंढने वाली exercise के रूप में छोड़ दिया गया? उस point तक मैं पहले ही अपनी math ability पूरी तरह झोंक चुका था, और थोड़ा “अब बाकी का उल्लू बना दो” जैसा लगा। क्या सिर्फ मुझे ऐसा लगा? article खुद वाकई अच्छा था।

    • पता नहीं मदद करेगा या नहीं: school में सीखे जाने वाले दो polynomials का multiplication असल में convolution होता है।
      “time domain में दो signals का convolution करना frequency domain में दो signals को multiply करने जैसा है” यह property है, और FFT हमें time domain से frequency domain में transform करने देता है। इसलिए polynomials को FFT से frequency domain में ले जाने के बाद, उस domain में बस multiplication करना होता है। यह convolution से तेज है। जानना चाहूंगा कि क्या इससे missing step स्पष्ट होता है; अगर कुछ छूटा हो तो article update कर सकता हूं।

  • तो क्या integer factorization discrete deconvolution है? सोच रहा हूं कि FFT representation, यानी pointwise multiplication के inverse operation, और tableax, यानी सामान्य long multiplication/carry addition को साथ-साथ रखने पर symmetry टूटती है और क्या fast algorithm के लिए पर्याप्त जानकारी मिल सकती है।

  • बेशक, naive polynomial multiplication polynomial degree के हिसाब से धीमा होता है। लेकिन असल में 100वीं degree के polynomial दो को handle करने की ज़रूरत कब पड़ती है?
    इसी वजह से मुझे ऐसा impression था कि computer algebra systems में यह तरीका इस्तेमाल नहीं होता

    • computer algebra system, जैसे Matlab का chebfun, arbitrary functions को 100वीं degree या उससे ऊपर के polynomials में बदल देता है, ताकि roots, optimal values वगैरह आसानी से निकाले जा सकें
    • error correction और signal processing में यह बहुत आम है
      https://www.youtube.com/watch?v=CcZf_7Fb4Us
      https://en.wikipedia.org/wiki/Reed%E2%80%93Solomon_error_cor... इसका एक उदाहरण है
    • बड़े files के CRC checksum parameters को reverse-engineer करना चाहता था, इसलिए मैंने एक program[1] बनाया जो file को लाखों degree वाले GF(2) polynomial में बदलकर greatest common divisor calculate करता है। FFT-based multiplication के बिना यह reasonable समय में संभव नहीं है
      [1]: https://github.com/8051enthusiast/delsum
    • यह convolution perspective और FFT के लिए fast GPU kernels, Mamba से पहले के कुछ state-space models में long sequence modeling के लिए इस्तेमाल हुए थे; यहाँ polynomial ही input sequence होता है
      Hazy Research के 2020~2023 के blog posts में इस approach के बारे में काफी जानकारी है
    • https://news.ycombinator.com/item?id=40306339 देखें
      “(...) physics research में लगभग 1 terabyte लंबी और 100 million से ज़्यादा terms वाली expressions से काम पड़ा था”