התמרת פורייה המהירה

בפוסט הקודם הצגתי את התמרת פורייה הבדידה, והבאתי דוגמה אחת לשימוש בה - כפל מהיר של פולינומים. אמרתי גם שכדי שהכפל באמת יהיה מהיר, עלינו לדעת לבצע את התמרת פורייה מהר; ואמרתי גם שהתמרת פורייה הבדידה מיוחדת בכך שהיא עוסקת בכמות סופית של מידע, ולכן אפשר לתכנת אותה במחשב בצורה פשוטה. כפי שאתם ודאי מנחשים מכל זה, הפוסט הנוכחי יהיה הרבה יותר מוטה לכיוון של מדעי המחשב מאשר הפוסטים הקודמים על פורייה; בפרט, אני הולך לדבר על שאלות של סיבוכיות ולהציג אלגוריתמים.

אבל ראשית, בואו ניזכר מהי התמרת פורייה הבדידה. אני הולך להשתמש בפוסט הזה בסימונים קצת שונים מאלו שהשתמשתי בהם בפוסט הקודם, כי הסימונים הישנים יהיו פחות נוחים לי. על כן אני ממליץ לאלו מכם שקראו את הפוסט הקודם לשכוח מהכל ולהתחיל מחדש. בכל הנוגע לסימונים, זאת אומרת.

אנחנו מסמנים את שורש היחידה הפרמיטיבי ה-\( n \)-י בסימון \( \omega_{n}=e^{\frac{2\pi i}{n}} \), והחזקות שלו \( \left(\omega_{n}^{0},\omega_{n}^{1},\dots,\omega_{n}^{n-1}\right) \) הן כל \( n \) שורשי היחידה מסדר \( n \). הקלט להתמרת פורייה הבדידה היא סדרה \( \left(a_{0},a_{1},\dots,a_{n-1}\right) \) של \( n \) מספרים מרוכבים (האינדקס מתחיל מ-0 כי אנחנו במדעי המחשב, כמובן) והפלט הוא סדרה \( \left(b_{0},b_{1},\dots,b_{n-1}\right) \) של סדרה של מספרים מרוכבים, כך שמתקיים הקשר הבא בין הסדרות, לכל \( 0\le k<n \):

\( b_{k}=\sum_{t=0}^{n-1}a_{t}\omega_{n}^{-kt} \) (משוואת אנליזה)

\( a_{k}=\frac{1}{n}\sum_{t=0}^{n-1}b_{t}\omega_{n}^{kt} \) (משוואת סינתזה)

למעשה, בפוסט הקודם המשוואות נראו טיפה אחרת (ה-\( \frac{1}{n} \) הופיע במשוואת האנליזה דווקא) אבל ההבדל אינו מהותי ומעכשיו המשוואות שיעניינו אותי הן אלו שהצגתי כאן.

במבט ראשון, לא נראה שיש בכלל צורך באלגוריתם כלשהו - המשוואות שמאפשרות לנו לחשב את \( b_{k} \) לכל \( k \) הן פשוטות מאוד. פשוט כופלים שורשי יחידה ב-\( a_{t} \)-ים ומחברים. מה הסיפור? ובכן, פשוט מאוד: יש \( n \) איברים \( b_{t} \) שאנחנו רוצים לחשב, ולכל אחד מהם נצטרך לבצע \( n \) פעולות של כפל ו-\( n \) פעולות של חיבור. סה”כ \( \Theta\left(n^{2}\right) \) פעולות. זה לא המון, אבל זה גם לא מעט. אם \( n=1000 \) זה אומר שנזדקק למיליון פעולות בערך.

מה שאנחנו רוצים הוא אלגוריתם שמבצע את כל המהומה הזו עם \( \Theta\left(n\log n\right) \) פעולות, כלומר כאשר \( n=1000 \) הוא מבצע בערך 7,000 פעולות - ההבדל בין זה ובין מיליון הוא די ברור. להשוואת זמני ריצה בפועל במימושים אמיתיים נגיע בסוף, אבל בינתיים תאמינו לי שזה הבדל גדול.

האלגוריתם שאציג לא יהיה סוף דבר - הוא יהיה שיטה נאיבית יחסית להתמודדות עם הבעיה. זה אומר שיש שיטות מורכבות יותר שאני לא מציג, והוא גם יניח הנחה מוזרה - ש-\( n \) הוא חזקה של 2. עם זאת, הוא ישיג לנו את זמן ה-\( \Theta\left(n\log n\right) \) המובטח ויתן תחושה של “איך עושים את הקסם הזה”, שזה מה שאני מחפש.

לפני שנצלול לפרטים, הנה הרעיון, וזה רעיון פשוט מאוד. ראשית, בואו נשים לב שאת משוואות האנליזה ניתן לתאר בתור הצבה של שורשי יחידה בפולינום. נגדיר את הפולינום \( A\left(x\right)=\sum_{t=0}^{n-1}a_{t}x^{t} \). כעת נשים לב לכך ש-\( b_{k}=A\left(\omega_{n}^{-k}\right) \). אז אפשר לשאול את השאלה הכללית: בהינתן פולינום \( A \) עם \( n \) מקדמים ו-\( n \) נקודות \( x_{0},\dots,x_{n-1} \), האם ניתן למצוא את \( y_{0},\dots,y_{n-1} \) המוגדרים על ידי \( y_{k}=A\left(x_{k}\right) \) ב-\( \Theta\left(n\log n\right) \) פעולות במקום ב-\( \Theta\left(n^{2}\right) \)?

זו שאלה טובה, ולמיטב ידיעתי התשובה היא לא (אבל אני לא מומחה ואין לי הוכחה כרגע). כלומר, אם אני מנסה לחשב את הערכים של \( A \) בנקודות שרירותיות כלשהן, אין לי דרך חכמה לעשות את זה. הקסם של התמרת פורייה המהירה הוא שבמקרה הספציפי של שורשי היחידה, דווקא כן אפשר לעשות קיצורי דרך בחישובים. זה אומר שהאלגוריתם יצטרך להסתמך איכשהו על תכונות מיוחדות שיש לשורשי היחידה ולסתם מספרים שרירותיים אין.

הרעיון הבסיסי באלגוריתם הוא פשוט מאוד: הפרד ומשול. הדוגמה הקלאסית שצריכה להיות לכם בראש היא זו של אלגוריתם מיון-מיזוג: הרעיון באלגוריתם המיון הזה הוא לקחת את הרשימה שממיינים, לפצל אותה לשניים, למיין כל תת-רשימה באופן רקורסיבי, ואז למזג את שתי הרשימות הממויינות שהתקבלו, תוך ניצול העובדה שקל יחסית למזג שתי רשימות ממויינות.

גם כן אנחנו הולכים לקחת את הסדרה שאנחנו רוצים לבצע לה התמרה, לפצל אותה לשתי תת-סדרות, למצוא התמרה לכל אחת בנפרד ואז “למזג” את שתי ההתמרות להתמרה אחת, יחסית ביעילות. ברמה הטכנית זה מתבצע בעזרת הפולינום \( A \) שתיארתי: אנחנו מפרקים את \( A \) לשני תת-פולינומים \( A^{0},A^{1} \), ובמקום לחשב את הערך של \( A \) על שורשי היחידה מסדר \( n \), אנחנו מחשבים את הערכים של שני תתי-הפולינומים על שורשי היחידה מסדר \( \frac{n}{2} \). את הערכים הללו אנחנו יכולים “למזג” אחר כך ולקבל את הערכים של \( A \) על שורשי היחידה מסדר \( n \). זה כל הסיפור; רק צריך להבין איך מתבצע שלב ה”מיזוג” במקרה הנוכחי. רק שימו לב לכך שראינו מדוע הכרחי להניח ש-\( n \) הוא חזקה של 2: אנחנו חייבים ש-\( n \) יתחלק ב-2 כדי הפירוק לשני תתי-פולינומים יעבוד ללא בעיות, ומכיוון שאנו פותרים באופן רקורסיבי גם \( \frac{n}{2} \) יצטרך לקיים את התכונה הזו וכן הלאה (עד למקרה של \( n=1 \) שהוא טריוויאלי - ההתמרה של סדרה מאורך 1 הוא הסדרה עצמה).

כלי הנשק המתמטי הבסיסי שלנו כאן הוא התכונה הבאה של שורשי היחידה מסדר \( n \), כאשר \( n \) הוא מספר זוגי: אם נעלה את כל שורשי היחידה מסדר \( n \) בריבוע, נקבל את כל שורשי היחידה מסדר \( \frac{n}{2} \), כאשר כל שורש יחידה שכזה מתקבל בדיוק פעמיים. הנה דוגמה קלה: שורשי היחידה מסדר 4 הם המספרים \( 1,i,-1,-i \). אם נעלה אותם בריבוע נקבל את המספרים \( 1,-1,1,-1 \), שהם שורשי היחידה מסדר 2, וכל אחד מהם התקבל בדיוק פעמיים. כבר מהדוגמה הזו מתקבל הרושם ש-\( \frac{n}{2} \) שורשי היחידה הראשונים, עד \( \omega_{n}^{\frac{n}{2}} \), כבר מספיקים כדי “לכסות” את כל שורשי היחידה מסדר \( \frac{n}{2} \), וש-\( \frac{n}{2} \) השורשים הבאים נותנים בדיוק את אותם ערכים (באותו סדר) כשמעלים אותם בריבוע.

בואו נוכיח את זה פורמלית. האבחנה הכללית היא ש-\( \omega_{dn}^{dk}=\omega_{n}^{k} \), וזה נובע ישירות מההגדרה: \( \omega_{dn}^{dk}=e^{\frac{2\pi idk}{dn}}=e^{\frac{2\pi ik}{n}}=\omega_{n}^{k} \). כמו כן, \( \omega_{n}^{k+n}=\omega_{n}^{k}\cdot\omega_{n}^{n}=\omega_{n}^{k} \) (כי \( \omega_{n}^{n}=1 \) שהרי \( \omega_{n} \) הוא שורש יחידה מסדר \( n \)). לכן:

מכאן נקבל לכל \( 0\le k<\frac{n}{2} \):

\( \left(\omega_{n}^{k+\frac{n}{2}}\right)^{2}=\omega_{n}^{2k+n}=\omega_{n/2}^{k+n/2}=\omega_{n/2}^{k}=\omega_{n}^{2k}=\left(\omega_{n}^{k}\right)^{2} \)

ראינו שהריבועים של \( \omega_{n}^{k} \) ושל \( \omega_{n}^{k+\frac{n}{2}} \) זהים, אבל למה הם שורש יחידה מסדר \( \frac{n}{2} \)? אה, זה קל: כי \( \left[\left(\omega_{n}^{k}\right)^{2}\right]^{\frac{n}{2}}=\left(\omega_{n}^{k}\right)^{n}=1 \). כעת, האם אנחנו באמת תופסים את כל שורשי היחידה מסדר \( \frac{n}{2} \) כך? שוב, די בבירור כן, כי \( \omega_{n/2}^{k}=\omega_{n}^{2k}=\left(\omega_{n}^{k}\right)^{2} \) (כאשר \( 0\le k<\frac{n}{2} \)).

אז אם לסכם: הריבועים של \( \omega_{n}^{0},\omega_{n}^{1},\omega_{n}^{2},\dots,\omega_{n}^{\frac{n}{2}-1} \) נותנים בדיוק את כל שורשי היחידה מסדר \( \frac{n}{2} \), והם זהים (גם בסדר שלהם) לריבועים של \( \omega_{n}^{\frac{n}{2}},\omega_{n}^{\frac{n}{2}+1},\dots,\omega_{n}^{n-1} \). עכשיו אפשר לעבור לאלגוריתם עצמו.

כזכור, אנחנו מתבוננים בפולינום \( A\left(x\right)=\sum_{t=0}^{n-1}a_{t}x^{t} \) ורוצים לחשב את ההצבה של שורשי היחידה בו. מה שנעשה יהיה לפרק אותו לשני פולינומים - זה של המקדמים במקומות הזוגיים, וזה של המקדמים במקומות האי זוגיים. נגדיר:

\( A^{0}\left(x\right)=\sum_{t=0}^{\frac{n}{2}-1}a_{2t}x^{t}=a_{0}x^{0}+a_{2}x^{1}+a_{4}x^{2}+\dots+a_{n-2}x^{\frac{n}{2}-1} \)

\( A^{1}\left(x\right)=\sum_{t=0}^{\frac{n}{2}-1}a_{2t+1}x^{2t+1}=a_{1}x^{0}+a_{3}x^{1}+a_{5}x^{2}+\dots+a_{n-1}x^{\frac{n}{2}-1} \)

שימו לב - החזקות של ה-\( x \)-ים הן רציפות, מ-1 ועד \( \frac{n}{2}-1 \), ולכן קיבלנו שני פולינומים ממעלה נמוכה בחצי מזו של \( A\left(x\right) \), שפשוט מחלקים ביניהם את המקדמים שלו.

עכשיו, אפשר לתאר את \( A\left(x\right) \) בעזרת שני הפולינומים הללו בצורה הפשוטה הבאה:

\( A\left(x\right)=A^{0}\left(x^{2}\right)+xA^{1}\left(x^{2}\right) \)

אם אתם לא מאמינים, נסו לבצע את החישוב בעצמכם, זה קל.

המסקנה פשוטה:

\( A\left(\omega_{n}^{k}\right)=A^{0}\left(\omega_{n}^{2k}\right)+\omega_{n}^{k}A^{1}\left(\omega_{n}^{2k}\right)=A^{0}\left(\omega_{n/2}^{k}\right)+\omega_{n}^{k}A^{1}\left(\omega_{n/2}^{k}\right) \)

וכפי שראינו, עבור \( k\ge\frac{n}{2} \) אין צורך לחשב את הערכים של \( A^{0},A^{1} \) על \( \omega_{n}^{k} \) במפורש, אפשר להשתמש בערך שלו על \( \omega_{n}^{k-\frac{n}{2}} \). זו הנקודה שבה אנחנו מנצלים את שורשי היחידה כדי לחסוך משהו.

הנה פסאודו-קוד של האלגוריתם. מכיוון שאני חסיד גדול של גישת “לתת פסאודו-קוד שבאמת אפשר להריץ”, הוא כתוב עבור מערכת האלגברה הממוחשבת Sage, שמשתמשת בשפה שהיא כמעט פייתון, עם כמה הרחבות (למשל, הסימן של החזקה). הדברים היחידים שנאלצתי לכתוב כאן ולא הייתי כותב בפסאודו קוד הם I במקום i בשביל השורש של מינוס 1, והמרה מפורשת של \( \frac{n}{2} \) למספר שלם, אחרת Sage לא מוכן לבצע פעולות מודולו \( \frac{n}{2} \).

def fft(a):
  n = len(a)
  if (n == 1):
    return a
  omega = e^(2*pi*I/n)
  a0 = [a[2*k] for k in range(n/2)]
  a1 = [a[2*k+1] for k in range(n/2)]
  b0 = fft(a0)
  b1 = fft(a1)
  b = [b0[k % Integer(n/2)] + omega^k*b1[k % Integer(n/2)] for k in range(n)]
  return b

כמובן, הקוד הזה מסתיר בתוכו מורכבות בסיסית ש-Sage מטפל בה אוטומטית אבל בשפות תכנות פשוטות אין אותה - אריתמטיקה עם שורשי היחידה \( \omega_{n}^{k} \). אנחנו מעלים אותם בחזקות, ומחברים אותם, ובלאגן שלם. מי שירצה לממש את התמרת פורייה יצטרך, כמובן, לממש קודם את החלק שמתעסק באריתמטיקה של המספרים הללו (ולכן של מספרים מרוכבים בכלל). לרוב פשוט מתבססים על ספריות קיימות.

קל לנתח את סיבוכיות האלגוריתם על ידי מבט בקוד. בשורות 6-7 (“הפיצול”) יש לנו לולאות עם זמן \( O\left(n\right) \) לכל אחת, וכך גם בשורה 10 (“המיזוג”). בשורות 8-9 יש לנו קריאה רקורסיבית לאלגוריתם עם קלט שגודלו חצי מהגודל הנוכחי. אז סיבוכיות האלגוריתם נתונה על ידי פתרון המשוואה \( T\left(n\right)=2T\left(\frac{n}{2}\right)+\Theta\left(n\right) \) - אותה משוואה של מיון מיזוג, שפתרונה הוא \( T\left(n\right)=\Theta\left(n\log n\right) \).

מכאן גם קל לממש את האלגוריתם עבור ההתמרה ההפוכה. אפשר לכתוב את כל הקוד מחדש, אבל בשביל מה? הנוסחה, כזכור, היא \( a_{k}=\frac{1}{n}\sum_{t=0}^{n-1}b_{t}\omega_{n}^{kt} \). היא מאוד דומה לנוסחה של התמרה רגילה. אולי אפשר לחשב את ההתמרה ההפוכה בעזרת שימוש בהתמרה הרגילה? יש שני הבדלים בין ההתמרה הרגילה וההפוכה. ראשית, יש את הכפל ב-\( \frac{1}{n} \); ושנית, יש את הסימן של החזקה של שורש היחידה: הנוסחה של ההתמרה הרגילה היא \( b_{k}=\sum_{t=0}^{n-1}a_{t}\omega_{n}^{-kt} \) ושם יש מינוס בחזקה, בעוד שכרגע אין.

התעלול המרכזי שנשתמש בו הוא שאנחנו כבר יודעים ש-\( \left(\omega_{n}^{k}\right)^{0}=\left(\omega_{n}^{k}\right)^{n} \) ולכן \( \left(\omega_{n}^{k}\right)^{-t}=\left(\omega_{n}^{k}\right)^{n-t} \). לכן אפשר לכתוב את נוסחת ההתמרה הרגילה גם כך:

\( b_{k}=\sum_{t=0}^{n-1}a_{t}\omega_{n}^{-kt}=\sum_{t=0}^{n-1}a_{t}\left(\omega_{n}^{k}\right)^{n-t}=\sum_{r=1}^{n}a_{n-r}\left(\omega_{n}^{k}\right)^{r} \)

כאן ביצעתי את החלפת המשתנה \( r=n-t \).

אם נסמן \( a_{n}\triangleq a_{0} \) אפשר לכתוב את הסכום האחרון שלמעלה בתור \( \sum_{r=0}^{n-1}a_{n-r}\left(\omega_{n}^{k}\right)^{r} \), שנראה קצת יותר מוכר.

אם כן, אם נתונה לי הסדרה \( c_{0},c_{1},\dots,c_{n-1} \) ואני רוצה לבצע עליה את ההתמרה ההפוכה, אני בעצם רוצה לחשב התמרת פורייה רגילה על הסדרה \( a_{0},\dots,a_{n-1} \) המוגדרת על ידי המשוואה:

\( \sum_{r=0}^{n-1}a_{n-r}\left(\omega_{n}^{k}\right)^{r}=\frac{1}{n}\sum_{t=0}^{n-1}c_{t}\left(\omega_{n}^{k}\right)^{t} \)

זה נותן לי את השוויונות הבאים:

\( a_{n-t}=\frac{1}{n}c_{t} \) כאשר \( 0\le t<n \)

כלומר, \( a_{t}=\frac{1}{n}c_{n-t} \) כאשר \( 0<t\le n \). עבור \( a_{0} \) כזכור מתקיים \( a_{0}=a_{n} \).

אם לסכם את מה שעשינו כאן, ה-\( a \)-ים שלנו מתקבלים מה-\( c \)-ים בצורה הבאה: את כל ה-\( c \)-ים כופלים ב-\( \frac{1}{n} \); את האיבר הראשון שלהם משאירים ללא שינוי; את הסדר של היתר הופכים. במילים אחרות, אם \( c=\left(1,2,3,4\right) \) אז נקבל \( a=\left(\frac{1}{4},1,\frac{3}{4},\frac{1}{2}\right) \). זה יוצא מאוד פשוט בקוד פייתון, שבו יש קונבנציה שאינדקסים שליליים ברשימה עוברים עליה מהסוף להתחלה:

def inverse_fft(b):
  n = len(b)
  return fft([b[k] / n for k in reversed(range(-(n-1),1))])

ומכאן קצרה הדרך לקוד עבור קונבולוציה:

def convolution(a,b):
  a_pad = a + [0 for i in range(len(a))]
  b_pad = b + [0 for i in range(len(b))]
  A = fft(a_pad)
  B = fft(b_pad)
  return inverse_fft([A[k]*B[k] for k in range(len(A))])[0:len(a)+len(b)-1]

הקוד הזה הוא נאיבי, במובן זה שהוא מניח ששתי הרשימות הן באותו האורך (ושהאורך הזה הוא חזקה של 2). יותר מזה, הוא גם עושה משהו מוזר - “מרפד” כל רשימה עם אפסים באופן שמכפיל את גודלה. למה? ובכן, כי שיקרתי קצת בפוסט הקודם: יצרתי את הרושם שקונבולוציה של שתי רשימות סופיות של מספרים אוטומטית מניחה שכל שאר המקומות הם 0, וזה אכן מה שצפוי שיהיה ומה שאנחנו משתמשים בו כשאנחנו כופלים פולינומים, למשל, אבל זה לא מה שהתמרת פורייה הבדידה עובדת איתו. זכרו שהתמרת פורייה הבדידה מתייחסת לרשימה של מספרים כאילו הם מייצגים פונקציה מחזורית, כלומר מחוץ לגבולות הרשימה הפונקציה מתחילה לחזור על עצמה. זה אומר שהקונבולוציה שאותה התמרת פורייה יודעת להפוך לכפל היא כזו שמניחה שהפונקציה היא מחזורית. הריפוד הוא פשוט דרך להתחמק מזה - אנחנו שמים אפסים בכל המקומות שעשויים להשפיע על החישוב, ואז מתעלמים מהאיברים ה”מיותרים” שנוצרו במהלך חישוב הקונבולוציה. זה מרגיש מלוכלך נורא, אבל זה למעשה לא פתרון כזה גרוע.

סיימנו לבינתיים עם FFT ועם התמרות פורייה למיניהן, אבל כמובן שהסיפור רק מתחיל כאן; פורייה נהיה מעניין כשמשתמשים בו בפועל. אני מקווה להציג בעתיד (אולי הלא רחוק) כמה שימושים מעניינים באמת.