Домівка > Математичний аналіз > Многочлени і перетворення Фур’є

Многочлени і перетворення Фур’є

Прямолінійний метод додавання двох многочленів вимагає \Theta(n) часу, а для множення — \Theta(n^2). Оскільки множення многочленів має значну практичну важливість, у цій статті ми покажемо як швидке перетворення Фур’є допомагає скоротити час на множення двох многочленів до \Theta(n \lg{n}).

Многочлени

Многочлен від змінної x над полем F представляє функцію A(x) як суму: \sum_{j=0}^{n-1} a_j x^j.

a_0, \dots, a_{n-1} ми називаємо коефіцієнтами многочлена. Ми кажемо, що A(x) має степінь k якщо найбільший індекс ненульового коефіцієнта k. Будь-яке ціле число більше ніж k є степінь-межею многочлена.

У випадку множення двох многочленів з векторами коефіцієнтів a і b отримуємо многочлен з вектором коефіцієнтів c. c також називають згорткою векторів a, b. І позначається як a\otimes b.

Представлення многочленів

Тут ми наведемо два види представлення многочленів: коефіцієнтне представлення і поточкове представлення.

Коефіцієнтне представлення

Коефіцієнтне представлення многочлена A(x) = \sum_{j=0}^{n-1} a_j x^j степінь-межі n це вектор коефіцієнтів a = (a_0, a_1, \dots, a_{n-1})^T. Таке представлення вигідне для певних дій на многочленами. Наприклад, операція обчислення многочлена A(x) у точці x_0. Використовуючи правило Горнера, ми можемо обчислити многочлен за \Theta(n) часу:
A(x_0) = a_0 + x_0(a_1 + x_0(a_2 + \dots + x_0(a_{n-2} + x_0 a_{n-1}))\dots)).
Так само, додавання в цій формі займає \Theta(n) часу.

Якщо ж розглянути множення двох многочленів A(x), B(x) степінь-межі n кожен, то ця операція має складність \Theta(n^2), оскільки нам необхідно помножити кожен з коефіціентів вектора a на кожен з коефіцієнтів вектора b. Результовний вектор коефіцієнтів c називається згортка.

Поточкове представлення

Поточкове представлення многочлена A(x) степінь-межі n є множиною двійок точка-значення
\{(x_0, y_0), (x_1, y_1), \dots, (x_{n-1}, y_{n-1})\}
такою, що всі x_k різні і y_k = A(x_k), k = 0, 1, \dots, n-1. Многочлен має багато різних поточкових представлень оскільки ми можемо використати будь-який набір точок x_0, x_1, \dots, x_{n-1} як базис для представлення.

Із використання метода Горнера обчислення многочлена в n точках займає \Theta(n^2) часу. Далі ми покажемо, що, якщо вправно вибрати n точок, то ми можемо прискорити обчислення до часу \Theta(n\lg{n}). Оберненою операцією для обчислення поточкового представлення є визначення коефіцієнтів многочлена — інтерполяція. Інтерполяція однозначно означена якщо многочлен має степінь на один меншу ніж кількість двійок точка-значення.

Інтерполяція зводиться до розв’язання системи лінійних рівнянь. Використовуючи LU-розклад ми можемо визначити коефіцієнти за час \O(n^3). Швидший алгоритм використовує формулу Лагранжа:
A(x) := \sum_{k=0}^{n-1} y_j \prod_{\begin{smallmatrix}0\le k\le n-1\\ j\neq k\end{smallmatrix}} \frac{x-x_j}{x_k-x_j}.

Поточкове представлення зручне для множення многочленів. Якщо C(x) = A(x)B(x), тоді C(x_k)=A(x_k)B(x_k) \forall x_k. Тобто, ми можемо поточково перемножити двійки точка-значення. Однак, ми повинні врахувати, що degree(C) = degree(A) + degree(B); якщо A і B степінь-межі n, тоді C степінь-межі 2n. Стандартне поточкове представлення для A, B містить n двійок точка-значення для кожного многочлена, отже ми повинні розширити їх поточкові представлення до 2n двійок точка-значення в кожному. Відтак маємо для A:
\{(x_0, y_0), (x_1, y_1), \dots, (x_{n-1}, y_{n-1})\},
для B:
\{(x_0, y'_0), (x_1, y'_1), \dots, (x_{n-1}, y'_{n-1})\},
для C:
\{(x_0, y_0y'_0), (x_1, y_1y'_1), \dots, (x_{n-1}, y_{n-1}y'_{n-1})\}.
Час необхідний на множення двох многочленів представлених у розширеній поточковій формі становить \Theta(n).

Швидке множення многочленів заданих коефіцієнтами

Ми можемо використати множення многочленів у поточковій формі, яке має лінійну обчислювальну складність для пришвидшення множення многочленів у коефіцієнтному представленні. Для цього ми повинні вміти швидко переводити між двома формами, тобто обчислювати й інтерполювати.

Викрут полягає в тому, що ми маємо право вибирати базис для точкового представлення як нам заманеться. Якщо правильно вибрати базис, то ми можемо змінювати представлення за час \Theta(\lg{n}). Пізніше ми покажемо, що для цього ми повинні використати комплексні корені з одиниці. Далі нам буде потрібне дискретне перетворення Фур’є і його обернене для переходу з представлення до представлення. А швидке перетворення Фур’є допоможе нам зробити це швидше.

Перед обчисленням многочленів A, B ми подвоюємо їхні степені до 2n додаванням нульових коефіцієнтів вищих порядків. Оскільки вектори мають 2n елементів, ми використаємо комплексні 2n корені з одиниці, які ми позначатимемо як w_{2n}.

Це наша ціль:

DFT та FFT

Комплексні корені одиниці

n-й комплексний корінь з одиниці це комплексне число w таке, що w^n = 1.

Існує саме n n-х комплексних коренів одиниці: e^{2\pi ik/n} для k = 0, 1, \dots, n-1. Щоб уявити собі таке число ми скористаємось формулою Ейлера, що пов’язує комплексну експоненту з тригонометричними функціями:
e^{iu}=\cos u + i \sin u.

На зображені видно комплексні корені при n рівному 8. w_n = e^{2\pi i/2} є головним значенням кореня.

n коренів утворюють абелеву групу щодо множення. Обернений елемент для кожного елементу цієї групи дорівнює спряженому елементу. Група коренів з одиниці ізоморфна адитивній групі лишків $latex\Z_n,$ оскільки w_n^n = w_n^0 = 1 тягне за собою, що w_n^j w_n^k = w_n^{j+k} = w_n^{(j+k)\mod n}. Так само, w_n^{-1} = w_n^{n-1}.

Лема скорочення Для будь-якого цілих чисел n\ge 0, k \ge 0, d > 0, w_{dn}^{dk} = w_n^k.

Лема половинення: Якщо n > 0 і парне, тоді квадрати n n-х комплексних коренів одиниці є n/2 n/2 коренів одиниці.

Лема підсумовування: Для будь-якого цілого n\ge 1 і ненульового цілого k неділимого на n, \sum_{j=0}^{n-1}(w_n^k)^j = 0.
Для доведення леми можна використати таку тотожність для \sum_{k=0}^nx^k = \frac{x^{n+1}-1}{x-1}, x \ne 1.

DFT

Повернемось до многочлена A(x) = \sum_{j=0}^{n-1}a_jx^j степінь-межі n, який ми хочемо обчислити в w_n^0, w_n^1, \dots, w_n^{n-1}. Тут ми вже маємо справу з розширеним многочленом, тобто під n ми розуміємо те, що раніше було 2n. Визначимо результовні y_k для k = 0, 1, \dots, n-1,
y_k = A(w_n^k) = \sum_{j=0}^{n-1}a_jw_n^{kj}.

Вектор y = (y_0, y_1, \dots, y_{n-1}) є дискретним перетворенням Фур’є (DFT) вектора коефіцієнтів a = (a_0, a_1, \dots, a_{n-1}). Ми позначатимемо y = \mbox{DFT}_n(a).

FFT

Завдяки використанню швидкого перетворення Фур’є (FFT), яке отримує перевагу завдяки використанню властивостей комплексних коренів одиниці, ми можемо обчислити \mbox{DFT}_n(a) за час \Theta(\lg n), тоді коли прямолінійний метод давав нам \Theta(n^2). Для спрощення, ми припускаємо, що n є точним степенем 2. Якщо це не так, ми можемо додати нульових коефіцієнтів вищих порядків. Хоча підходи для обробки ситуації коли n не є точним степенем двійки відомі.

FFT метод використовує стратегію розділяй і володарюй, порізно використовуючи парно- і непарноіндексовані коефіцієнти A(x) для визначення двох нових многочленів степінь-межі n/2:
A^{[0]}(x) = a_0 + a_2x + a_4x^2 + \dots + a_{n-2}x^{n/2-1},
A^{[1]}(x) = a_1 + a_3x + a_5x^2 + \dots + a_{n-1}x^{n/2-1}.

Маємо, що A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2), отже проблема обчислення A(x) у w_n^0, w_n^1, \dots, w_n^{n-1} зводиться до
1. обчислення многочленів A^{[0]}(x), A^{[1]}(x) степінь-межі n/2 в точках
(w_n^0)^2, (w_n^1)^2, \dots, (w_n^{n-1})^2, (1)
2. поєднання результатів.

Згідно з лемою половинення, набір значень (1) містить не n різних значень, але n/2 n/2 коренів одиниці, кожен корінь зустрічається двічі. Отже, на кожному етапі рекурсії ми отримуємо вдвічі більше задач вдвічі меншого розміру. Така декомпозиція становить базу для наступного рекурсивного алгоритму.
RECURSIVE-FFT(a)
n = a.length // n є степенем 2
якщо n == 1
    повернути a
w_n = e^{2\pi i/n}
w = 1 // встановлення w_n, w необхідно для правильного обчислення w у циклі
a^{[0]} = (a_0, a_2, \dots, a_{n-2})
a^{[1]} = (a_1, a_3, \dots, a_{n-1})
y^{[0]} = RECURSIVE-FFT(a^{[0]})
y^{[1]} = RECURSIVE-FFT(a^{[1]})
для k = 0 \dots n/2-1
    y_k = y_k^{[0]} + wy_k^{[1]}
    y_{k+n/2} = y_k^{[0]} - wy_k^{[1]}
    w = ww_n
повернути y // y – стовпчик-вектор

Два виклики RECURSIVE-FFT обчислюють \mbox{DFT}_{n/2},
y_k^{[0]} = A^{[0]}(w_{n/2}^k),
y_k^{[1]} = A^{[1]}(w_{n/2}^k),
або, відповідно до леми скорочення,
y_k^{[0]} = A^{[0]}(w_n^{2k}),
y_k^{[1]} = A^{[1]}(w_n^{2k}).

У тілі циклу ми комбінуємо результати отримані від двох рекурсивних викликів. Отже, для y_0, y_1, \dots, y_{n/2-1}
y_k = y_k^{[0]} + w_n^k y_k^{[1]}
    = A^{[0]}(w_n^{2k})+w_n^kA^{[1]}(w_n^{2k})
    = A(w_n^k).

Для y_{n/2}, y_{n/2+1}, \dots, y_{n-1} маємо
y_{k+n/2} = y_k^{[0]} - w_n^k y_k^{[1]}
    = y_k^{[0]} + w_n^{k+n/2} y_k^{[1]}
    = A^{[0]}(w_n^{2k}) + w_n^{k+n/2} A^{[1]}(w_n^{2k})
    = A^{[0]}(w_n^{2k+n}) + w_n^{k+n/2} A^{[1]}(w_n^{2k+n})
    = A(w_n^{k+n/2}).
Тут ми використали, що -w_n^k = w_2^1 w_n^k = w_n^{n/2} w_n^k = w_n^{k + n/2}.

Щоб визначити час виконання процедури RECURSIVE_FFT зауважимо, що окрім самих рекурсивних викликів виконується \Theta(n) роботи. Отже, рекурентна формула така:
T(n) = 2 T(n/2) + \Theta(n) = \Theta(n\lg{n}).

Значить, ми можемо обчислити многочлен степінь-межі n в n комплексних коренях одиниці за час \Theta(n\lg{n}).

Інтерполяція в комплексних коренях одиниці

Зараз ми завершимо схему множення многочленів, показавши як інтерполювати комплексні корені одиниці многочленом, що уможливить перетворення з форми точка-значення до форми коефіцієнтів. Запишемо DFT як матричне рівняння y=V_na і знайдемо форму оберненої матриці.
\begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ \vdots \\ y_{n-1} \end{pmatrix} =   \begin{pmatrix}   1 & 1 & 1 & \dots & 1 \\   1 & w_n & w_n^2 & \dots & w_n^{n-1} \\   1 & w_n^2 & w_n^4 & \dots & w_n^{2(n-1)} \\   \vdots & \vdots & \vdots & \ddots & \vdots \\   1 & w_n^{n-1} & w_n^{2(n-1)} & \dots & w_n^{(n-1)(n-1)} \\   \end{pmatrix}  \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ \vdots \\ a_{n-1} \end{pmatrix}

(k, j) елемент матриці V_n є w_n^{kj}, для k = 0,1,\dots,n-1. Для зворотньої операції, яку ми запишемо як a = \mbox{DFT}_n^{-1}(y), ми помножимо y на матрицю V_n^{-1}.

Теорема: Для j,k = 0,1,\dots,n-1, (k, j) елемент з V_n^{-1} є w_n^{-kj}/n.
*Доведення:** Розглянемо (j, j') елемент матриці V_n^{-1}V_n:
[V_n^{-1}V_n]_{jj'} = \sum_{k=0}^{n-1}(w_n^{-kj}/n)(w_n^{kj'})=\sum_{k=0}^{n-1}(w_n^{k(j'-j)}/n.
Сума дорівнює 1 якщо j'=j і 0 інакше завдяки лемі підсумовування.

Отже, маючи матрицю V_n^{-1}, отримуємо, що \mbox{DFT}_n^{-1}(y) задається через
a_j = \frac{1}{n}\sum_{k=0}^{n-1}y_k w_n^{-kj}, j=1,2,\dots,n-1.
Звернімо увагу, що ця формула дуже схожа з формулою y = \mbox{DFT}_n(a). Ми лише поміняли ролями x і a, замінили w_n на w_n^{-1}, і поділили кожний елемент результату на n. з цього видно, що ми можемо обчислити \mbox{DFT}_n^{-1} за \Theta(n\lg{n}).

Теорема про згортку

Для будь-яких двох векторів a, b довжини n, де n є степенем 2, згортка векторів
a\otimes b = \mbox{DFT}_{2n}^{-1}(\mbox{DFT}_{2n}(a)\cdot\mbox{DFT}_{2n}(b)),
де вектори a,b доповнені нулями до довжини 2n і \cdot позначає покомпонентний добуток двох векторів.

Приклади

Обчислення DFT для вектора (0,1,2,3)

a = (0,1,2,3), n = 4.
w_4^0 = 1
w_4^1 = i
w_4^2 = -1
w_4^3 = -i
y_0 = 0*w_4^{0*0} + 1*w_4^{0*1} + 2*w_4^{0*2} + 3*w_4^{0*3} = 6

\dots
y_3 = 0*w_4^{3*0} + 1*w_4^{3*1} + 2*w_4^{3*2} + 3*w_4^{3*3} = 0*w_4^0 + 1*w_4^3 + 2*w_4^2 + 3*w_4^1 = -2+2i
y = (6, -2-2i, -2, -2+2i)
Якщо попросити MATLAB порахувати DFT цього ж вектора, то результат буде такий:

>> fft([0 1 2 3])
ans =
   6.0000 + 0.0000i  -2.0000 + 2.0000i  -2.0000 + 0.0000i  -2.0000 - 2.0000i

Відмінні знаки уявного доданку завдячують тому, що MATLAB як w_n використовує e^{(-2\pi i)/n}, а в цій статті ми використовуємо e^{(2\pi i)/n}.

Дивись також

У статті використані матеріали з книги “Вступ в алгоритми” за авторством Томаса Кормена, Чарльза Лейзерсона, Рональда Рівеста та Кліффорда Штайна.

Advertisements
  1. Коментарів ще немає.
  1. No trackbacks yet.

Залишити відповідь

Заповніть поля нижче або авторизуйтесь клікнувши по іконці

Лого WordPress.com

Ви коментуєте, використовуючи свій обліковий запис WordPress.com. Log Out / Змінити )

Twitter picture

Ви коментуєте, використовуючи свій обліковий запис Twitter. Log Out / Змінити )

Facebook photo

Ви коментуєте, використовуючи свій обліковий запис Facebook. Log Out / Змінити )

Google+ photo

Ви коментуєте, використовуючи свій обліковий запис Google+. Log Out / Змінити )

З’єднання з %s

%d блогерам подобається це: