Ранг матрицы — это или число её линейно независимых строк, или число её линейно независимых столбцов, или порядок её максимального минора с ненулевым определителем. Любого из этих трёх равнозначных определений ранга вполне достаточно на все случаи жизни. При решении задач мы будем иметь дело только с квадратными матрицами, поэтому о прямоугольных матрицах здесь говорить не будем.

Если дана матрица A размером n×n, то её ранг — это целое число от нуля до n. Рассмотрим крайние варианты. Если A — нулевая матрица, т.е. состоит из одних нулей, то rank A = 0. Если A — матрица с ненулевым определителем, т.е. det A ≠ 0, то rank A = n, поскольку в этом случае максимальный невырожденный минор — это сама матрица A и есть. Во всех остальных случаях 1 ≤ rank A ≤ n–1.

Обычно ранг матрицы находят используя слегка модифицированные приёмы, полезные при нахождении определителя (далее эти приёмы даются для столбцов; во всех формулировках столбцы можно заменить на строки):

1. Все элементы столбца можно умножить на любое число, отличное от нуля. Например, можно сократить все элементы столбца на общий множитель.
2. К одному столбцу можно прибавить другой, умноженный на любое число.
3. Нулевой столбец можно выкинуть (ранг от этого не изменится). Из двух одинаковых или пропорциональных столбцов один можно выкинуть.
4. Столбцы можно переставлять.

Поговорим теперь о ранге матриц вида A – λE, где λ — собственное число матрицы A. Для таких матриц det(A – λE) = 0, поэтому обязательно rank(A – λE) < n.

Мне хочется привести ещё один способ решения этой задачи, основанный на операционном исчислении. Подобно тому, как дифференциальное исчисление сформировалось при исследовании понятия «производная», а интегральное исчисление — при исследовании обратного к производной действия, так и операционное исчисление — это теория, базирующаяся на двух понятиях: «оригинал» и «изображение».

Придуманное английским физиком Оливером Хевисайдом на рубеже XIX и XX веков, операционное исчисление поначалу сильно раздражало ортодоксальных математиков, поскольку Хевисайд получал правильные результаты при решении важных прикладных задач, нимало не беспокоясь насчёт обоснования своих выкладок, которые, надо сказать, выглядели весьма непривычно. Лишь после трудов Бромвича и Карсона из Англии, а также польского математика Яна Микусинского, которые подвели теоретическую базу под идеи Хевисайда, операционное исчисление приобрело статус «законной» области математики.

Определение. Пусть  f (t) — некоторая функция, заданная на положительной полуоси, и интеграл

∞ F(p) =   ∫  e–pt f (t) dt 0

(5)

при некотором значении параметра p=p₀ абсолютно сходится.² Тогда сам интеграл называется преобразованием Лапласа функции  f (t). При этом F(p) называется изображением функции  f (t), а  f (t) — оригиналом функции F(p).

Чтобы не записывать каждый раз формулу (5) для указания на вид связи между функциями  f (t) и F(p), была придумана форма записи  f (t) _·=^· F(p), причём она используется в обе стороны: и  f (t) _·=^· F(p), и F(p) _·=^· f (t). Только аргументы показывают «кто есть кто» в этих соотношениях: если аргумент p, то речь идёт об изображении; если аргумент t, то это оригинал.

Для многих функций  f (t) преобразование Лапласа находится достаточно легко, что позволяет без труда составлять пары «оригинал–изображение».

[· · ·]

Возвращаемся от изображений к оригиналам и получаем ответ:

Напоследок пара советов. Если захотелось узнать об операционном исчислении чуть больше, то не надо кидаться с криком на диван: «Полежу полчаса — всё пройдёт!», а наоборот — надо отправляться в библиотеку за книгами:

1. В. А. Диткин,  А. П. Прудников. «Операционное исчисление» (М., Высшая школа, 1975).
2. Ф. А. Шелковников,  К. Г. Такайшвили. «Сборник упражнений по операционному исчислению» (М., Высшая школа, 1976).

Изучение лучше начинать со второй. Читаем условия задач, затем лезем в ответы и читаем решения (в этой книге они довольно подробные). Через какое-то время начинаем самостоятельно решать задачи, у которых приводятся только ответы. Таким образом можно пройти достаточно большой путь без теоремы обращения и обобщённых функций. Когда захочется разобраться с обобщёнными функциями, то желательно для начала прочитать Добавление в книге: Я. Б. Зельдович. «Высшая математика для начинающих» (М., Наука, 1970). Её автор рассказывает о дельта-функции Дирака очень доходчиво и не стремится напугать читателя многоэтажной математикой.

Матрицы в последних двух примерах были частными случаями n×n-матрицы

Давно понятно, что, несмотря на дифурные формулировки задач, занимаемся мы чистой воды алгеброй. Поэтому при исследовании матрицы A на собственные числа и векторы не будем утруждать себя даже такой малостью, как запись системы уравнений.

Корень из единицы ζ = e^2πi/n — это корень уравнения xⁿ – 1 = 0. Многие свойства этого числа можно получить непосредственно из уравнения. Например, поскольку

и ζ ≠ 1, то 1 + ζ + ζ ² + ... + ζ ^n–1 = 0. Более того, поскольку числа ζ ^k, k = 2, n – 1, наряду с ζ, также являются корнями уравнения xⁿ – 1 = 0, то

	n–1
1 + ζ k + (ζk)2 + ... + (ζ k)n–1 =	∑	(ζ k)m = 0     для всех   k = 1, n – 1.
	m=0

И всё же, если есть возможность выбора, я предпочитаю использовать методы анализа, а не алгебраические средства. В частности, получать сумму степеней корней из единицы через геометрическую прогрессию:

n–1 ∑ m=0 ζ km =  n–1 ∑ m=0 (ζ k)m =  ζ kn – 1 ζ k – 1  =  ì í î n,    k ≡ 0 (mod n), 0,    k ≠ 0 (mod n).

(13)

Здесь k может быть любым целым числом: и положительным, и отрицательным, и нулём.

Простыми следствиями формулы (13) являются равенства

n–1 ∑ m=0	(ζ k)m·(ζ l)m =	ì í î	n,    k + l ≡ 0 (mod n), 0,    k + l ≠ 0 (mod n);
n–1 ∑ m=0	(ζ k)m·(ζ l)m =	ì í î	n,    k – l ≡ 0 (mod n), 0,    k – l ≠ 0 (mod n).

Первое из них полезно при нахождении квадрата матрицы A, второе — при нахождении произведения AA:

Эти формулы показывают, что матрица U = A/√n — более удобный объект, чем первоначальная матрица A. Во-первых, UU = E и, таким образом, U ^–1 = U. Во-вторых, U — унитарная матрица⁵ и можно воспользоваться плодами общей теории. В-третьих, U ² состоит только из нулей и единиц, причём нулей в ней ну очень много. Это позволяет надеяться на достаточно простой вид характеристического многочлена для U ², а там и до собственных чисел матрицы U рукой подать.⁶

Итак, ищем собственные числа матрицы U ²:

[· · ·]

det (U 2 – λE) =

ì í î

(λ – 1)n1 + 1(λ + 1)n1 – 1,  –(λ – 1)n1 + 1(λ + 1)n1,     n = 2n1;    n = 2n1 + 1.

Таким образом, собственные числа n×n-матрицы U ² — это числа +1 и –1, а их кратности зависят от чётности n. Соответственно, для матрицы U собственными числами будут квадратные корни из +1 и –1, т.е. числа ±1 и ±i. Поэтому

Насчёт кратностей a, b, c, d пока можно сказать только одно:

a + b =    кратность числа  +1 для матрицы U 2    =

ì í î

n1 + 1, n1 + 1,   n = 2n1;   n = 2n1 + 1;

c + d =    кратность числа  –1 для матрицы U 2    =

ì í î

n1 – 1, n1,   n = 2n1;   n = 2n1 + 1.

Чтобы найти a, b, c, d нужно выявить ещё какие-нибудь связи между этими числами. Используем следующие соображения. Матрица U унитарная, её жорданова матрица J диагональна, и на этой диагонали располагаются +1 в количестве a штук, –1 в количестве b штук, а также +i и –i по c и d штук соответственно. В сумме a+b+c+d дают, естественно, n. Далее,

поскольку, в силу общей теории, след и определитель матрицы являются инвариантами. Определитель det U вычисляется легко (определитель Вандермонда), но однозначного нахождения чисел a, b, c, d не гарантирует. Зато след tr U идеально подходит для нашей цели:

				n–1
a – b + ci – di = tr U =	1  √n	tr A =	1  √n	∑	ζm².
				m=0

Последняя сумма называется суммой Гаусса. Именно Гаусс впервые доказал, что

n–1  n–1  ∑  ζm² =  ∑  e2πim²/n = m=0 m=0 ì í î (1 + i)√n,  √n,  0,  i√n,    n ≡ 0(mod 4);    n ≡ 1(mod 4);    n ≡ 2(mod 4);    n ≡ 3(mod 4).

(15)

Поэтому

a – b =

ì í î

1, 1, 0, 0,   n ≡ 0(mod 4);   n ≡ 1(mod 4);   n ≡ 2(mod 4);   n ≡ 3(mod 4);

и     c – d =

ì í î

1, 0, 0, 1,   n ≡ 0(mod 4);   n ≡ 1(mod 4);   n ≡ 2(mod 4);   n ≡ 3(mod 4).

что, совместно с найденными выше значениями для a+b и c+d, приводит к равенствам:

a =

1 4

ì í î

n + 4 n + 3 n + 2 n + 1

÷ ÷ ÷

,   b =

1 4

ì í î

n n – 1 n + 2 n + 1

÷ ÷ ÷

,   c =

1 4

ì í î

n n – 1 n – 2 n + 1

÷ ÷ ÷

,   d =

1 4

ì í î

n – 4 n – 1 n – 2 n – 3

÷ ÷ ÷

при

n ≡ 0(mod 4);   n ≡ 1(mod 4);   n ≡ 2(mod 4);   n ≡ 3(mod 4).

Таким образом,

n + 4 4     n + 2 4     n + 1 4     n – 1 4   det(U – λE) = (–1)n (λ – 1)  (λ + 1)  (λ – i)  (λ + i) .

(16)

Аналогичное разложение для det(A – λE) получается заменой скобок (λ ± 1), (λ ± i) на скобки (λ ± √n), (λ ± i√n), при сохранении показателей степеней. Для n = 3 и n = 4, в частности, получаются уже известные формулы

Задача нахождения собственных чисел матрицы A почти завершена: для порядка надо бы доказать формулу (15). Займёмся этим.⁷

Наиболее удачный способ доказательства (удачный в том смысле, что способ не требует никакой специальной изобретательности) предложил в 1835 году Дирихле, который впоследствии использовал сумму Гаусса в своей знаменитой теореме о простых числах в арифметических прогрессиях: если a и b — целые, взаимно простые числа, то последовательность

содержит бесконечно много простых чисел.

Итак, доказательство Дирихле. Пусть функция g(x) задана на отрезке [0, n] явно: g(x) = e^2πix²/n, а за пределы этого отрезка продолжена периодическим образом: g(x+n) = g(x). Получившаяся функция, как легко видеть, непрерывна, но не дифференцируема в точках склейки. Например, g(0) = g(n) = 1, g′(0+) = 0 ≠ g′(0–) = g′(n–) = 4πi.

Разложим g(x) в ряд Фурье⁸ по системе функций {e^2πikx/n, k Î R}:

[· · ·]

Тогда

n–1  n–1 ∑ ζm² = ∑ g(m) = [· · ·] = n(1 + i–n)  ∫  e2πint² dt =  (1 + i–n)(1 + i) 2 √n. m=0 m=0 R

На завершающем этапе было использовано известное значение интеграла Френеля:

∫   eiat² dt =  Ö π  a  eπi/4 =  Ö π 2a  (1 + i)       (a > 0). R

(17)

Другой способ вычисления суммы Гаусса был предложен в «Курсе математического анализа» Гурса (М., Гостехиздат, 1933, т. 2, ч. 1, стр. 120) и основан на контурном интегрировании.⁹

[· · ·]

Можно ещё отметить, что дважды использованный на последней стадии выкладок трюк ∫_a^b + ∫_b^c = ∫_a^c — это не свойство интеграла по вещественным отрезкам, а теорема Коши для голоморфной функции и треугольного контура. Экспонента, которую мы интегрировали, — целая аналитическая функция, поэтому можно брать треугольники с любыми R и n.

Ещё один способ вычисления суммы Гаусса предложил в 1840 году В. А. Лебег¹⁰ на основе формулы Пуассона

∑   e–π(k+t)²z =   1  √ z  ∑   e2πikt e–πk²/z. kÎZ kÎZ

(18)

Здесь tÎR, zÎC, причём Re z > 0 (чтобы ряды сходились), а при вычислении √ z выбирается то значение, для которого Re √ z > 0.

Формула (18) обычно получается как частный результат или при рассмотрении преобразования Фурье обобщённых функций,¹¹ или при рассмотрении дробно-линейных преобразований тэта-функций.¹² Для наших целей, однако, проще проигнорировать эти теории и доказать равенство (18) непосредственно.

[· · ·]

Вернёмся теперь к формуле Пуассона и посмотрим, что с ней произойдёт, если положить z = ε – im/n, где m, n Î N, а затем устремить ε к нулю. Очевидно, что мы не можем просто взять и подставить z = –im/n в равенство (18), поскольку ряды тогда разойдутся. Весь интерес заключается в том, чтобы, опираясь на абсолютную и равномерную сходимость рядов при ε>0, выяснить их асимптотическое поведение. Грубо говоря, если ряд слева ведёт себя при ε → +0 как (A+o(1))/ε, а ряд справа — как (B+o(1))/ε, то A=B.

[· · ·]

Приравнивая главные члены асимптотик (19) и (20), меняя индекс суммирования k₂ на k и считая, что mn + l — чётное число, получаем в итоге плоды наших трудов:

n–1 m–1  1  √ n  exp ( πil 2 4mn ) ∑  exp ( πim n  k2 +  πil n k ) =   1  √ m  exp (  πi 4 ) ∑  exp ( –   πin m  k2 +  πil m k ) . k=0 k=0

(21)

В частности, при l = mn

n–1 m–1  1  √ n ∑  (–1)mk eπimk²/n =  eπi(1 – mn)/4  √ m  ∑  (–1)nk e–πink²/m. k=0 k=0

(22)

Если теперь положить в этом равенстве m = 2, то получим

	n–1
1  √ n	∑	e2πik²/n =	eπi/4  √ 2	e–πin/2(1 + (–1)n e–πin/2) =	(1 + i)(1 + i –n) 2	,
	k=0

и тем самым выражение для суммы Гаусса найдено в третий раз. Как говорил Гекке: «Точное знание поведения аналитической функции в окрестности её особых точек — источник многих арифметических теорем».

Кстати, если обозначить сумму в левой части равенства (22) через g(n, m), то легко видеть, что само равенство представимо в виде:

т.е. функция g(n, m) хоть и не симметрическая, но очень терпимо относится к перестановке своих аргументов. Такого рода формулы, когда некая функция  f (n, m) связана с  f (m, n) простым соотношением, называются теоремами или законами взаимности, поэтому равенство (22) выражает закон взаимности для обобщённой суммы Гаусса.¹³

Все приведённые выше способы доказательства формулы (15) являются аналитическими по методам и получаются с использованием непрерывных функций, несобственных интегралов, рядов Фурье и т.д., хотя исходная сумма была конечным набором дискретных объектов и, казалось бы, должна находиться в рамках манипуляций подобными же объектами. Хассе, занимаясь в своих «Лекциях по теории чисел» (М., ИЛ, 1953) вычислением суммы Гаусса и приводя «почти арифметический» вариант вычисления, пишет, что вклад анализа хоть и можно свести к минимальному, но полностью обойтись без него нельзя. И далее с явным неудовольствием продолжает: «Существуют другие доказательства, в которых, наоборот, арифметике отводится по возможности меньшая роль, а иногда она полностью вытесняется аналитическими рассуждениями». Я же, в свою очередь, при любом удобном случае пытаюсь избавиться от дискретных объектов и «обменять» их на непрерывные, поскольку арсенал методов работы с последними гораздо богаче.

Итак, с собственными числами мы разобрались,¹⁴ займёмся теперь собственными векторами. Их нахождение — более тяжёлая задача. У меня, по крайней мере, простого решения этой задачи нет. Если вы смогли найти такое решение или знаете по литературе, то занесите его мне на кафедру.

Начнём с некоторых замечаний общего плана. Поскольку A = √nU, то собственные векторы матриц A и U совпадают. Уже было отмечено, что U — унитарная матрица, поэтому для каждого её собственного числа λ кратности #(λ) найдётся ровно сколько же ЛНЗ собственных векторов. Вообще, предварительное обсуждение — это удобный момент, чтобы вспомнить кое-какие определения и теоремы из алгебры.¹⁵

• Матрица S называется симметрической, если S = S ^t. Матрица H называется эрмитовой, если H = H ^t (при этом H ^t обозначается H^*). Вещественная матрица O называется ортогональной, если OO ^t = E. Комплексная матрица U называется унитарной, если UU ^* = E. Матрица N называется нормальной, если NN ^* = N ^*N (в частности, эрмитова матрица нормальна, унитарная матрица тоже нормальна).
• Для вещественной симметрической матрицы S собственные числа вещественны, а собственные векторы, соответствующие различным собственным числам, ортогональны в евклидовом пространстве Rⁿ. Если λ — собственное число кратности m, то для него существует ровно m ЛНЗ собственных векторов.¹⁶ Как следствие, из собственных векторов вещественной симметрической n×n-матрицы можно собрать базис пространства Rⁿ. Существует ортогональная матрица O, приводящая S к диагональному виду:
  S = OΛO ^–1
  

  где Λ = diag(λ₁, λ₂,..., λ_n), а λ_k, k = 1, n, — все собственные числа матрицы S с учётом кратности.

Легко заметить соответствие «симметрическая матрица → эрмитова матрица», «ортогональная матрица → унитарная матрица». Поэтому не удивительно, что комплексная переформулировка предыдущего пункта делается простой заменой. И всё же наша исходная матрица A комплексна, поэтому нас, в первую очередь, должен интересовать именно комплексный вариант теории, а о вещественных матрицах я тут всего наговорил просто так, для преемственности перехода от Rⁿ к Cⁿ.

Сначала обсудим нормальные матрицы, как наиболее общие среди определённых выше.

• Любая нормальная матрица N приводится к диагональному виду с помощью некоторой унитарной матрицы: N = UΛU ^–1. Если Nh = λh (т.е. λ — собственное число матрицы N, а h — соответствующий собственный вектор), то N ^*h = λh. (Как следствие, все собственные числа эрмитовой матрицы вещественны.) Для нормальных матриц собственные векторы, соответствующие различным собственным числам, ортогональны в унитарном пространстве Cⁿ. Если λ — собственное число матрицы n кратности m, то найдётся m ЛНЗ собственных векторов, соответствующих этому числу. Найденные собственные векторы можно ортогонализовать, затем нормировать и получить в результате ортонормированный базис пространства Cⁿ. ¹⁷

Теперь можно было бы сказать несколько слов об унитарных матрицах. Например, что унитарные матрицы в пространстве Cⁿ осуществляют перевод ортонормированного базиса в ортонормированный, что собственные числа унитарных матриц по модулю равны единице, т.е. λ = e^ic, где c Î R, но, пожалуй, пора переходить к нашей конкретной матрице U = A/√n. Собственные числа матрицы U равны {±1, ±i}, кратности тоже известны, а про собственные векторы мы теперь знаем, что из них можно собрать базис в Cⁿ. Кроме того, для разных собственных чисел собственные векторы ортогональны. Поэтому пространство Cⁿ можно разложить в прямую сумму ортогональных подпространств, каждое из которых является линейной оболочкой собственных векторов, соответствующих одному из собственных чисел.

[· · ·]

Взглянем теперь на проблему поиска собственных векторов матрицы A с другой точки зрения. Если x, y Î Rⁿ, то y = Ax и x = Ay/n называют соответственно прямым и обратным дискретным преобразованием Фурье (или просто ДПФ). В координатной форме записи эти преобразования выглядят так:

	n–1			n–1
yk =	∑	e2πikm/n xm     (k = 0, n–1)     и     xm =	1  n	∑	e–2πikm/n yk     (m = 0, n–1),
	m=0			k=0

Если не применять специальных ухищрений, то для нахождения вектора y по известному вектору x требуется n² комплексных умножений. Можно слегка выгадать, если не делать умножений на 1, но порядок это не уменьшит. На практике, чтобы его уменьшить, применяют алгоритмы быстрого преобразования Фурье (или просто БПФ), которым для нахождения вектора y достаточно порядка n ln n комплексных умножений.

По сути все алгоритмы БПФ заключаются в том, чтобы выбрать составное число n, разложить его на множители, n = n₁n₂, и перейти от задачи умножения n×n-матрицы A на n-мерный вектор x к подзадачам размерности n₁ и n₂. При этом структура корней из единицы ζ^k позволяет реализовать такой переход весьма эффективно с точки зрения времязатрат.

Пусть, например, n — чётное число. Тогда

[· · ·]

Тем самым n-точечная задача вычисления ДПФ свелась к двум n/2-точечным задачам. При этом вектор x тоже разделился на две составляющие — чётную и нечётную. На этом этапе уже ощущается потребность введения более детальной информации в обозначениях. Раньше для корней из единицы запись ζ = e^2πi/n казалась вполне достаточной. Теперь, чтобы не путаться в формулах, в которых присутствуют корни из единицы разных степеней, я иногда буду использовать нижние индексы: ζ_n = e^2πi/n, A_n = ||ζ_n^kl||, U_n = A_n /√n. Таким образом, формулы (24) осуществляют редукцию от n-точечного ДПФ y = A_nx к двум n/2-точечным с матрицей A_n/2.

Обозначим через M(n) число комплексных умножений при вычислении ДПФ y = A_nx. Тогда формулы (24) позволяют составить рекуррентное соотношение M(n) = 2M(n/2) + n/2. Если n является степенью двойки, то процесс редукции можно довести до логического конца, n → n/2 → n/4 → ... → 1, не прибегая ни к чему иному, кроме полученных формул. При этом

т.е. число необходимых комплексных умножений существенно меньше, чем при прямом перемножении матрицы A_n на вектор x.

Этот простой и наглядный пример реализации БПФ, который был опубликован Кули и Тьюки в 1965 году, не следует рассматривать как единственно возможный. Блейхут в своей книге «Быстрые алгоритмы цифровой обработки сигналов» (М., Мир, 1989) называет «злополучным мифом» широко распространившееся мнение о том, что ДПФ становится быстрым лишь при длине блока n, равной степени двойки. На самом же деле хорошие БПФ-алгоритмы существуют практически для произвольной длины блока.

Другим алгоритмом БПФ, сильно отличающимся от алгоритма Кули–Тьюки и гораздо менее известным, является алгоритм Гуда–Томаса. Он был опубликован в 1960 году, но публикация эта прошла почти незамеченной. Для нашей задачи поиска собственных векторов матрицы A_n именно этот алгоритм представляет интерес.

Начнём издалека — с линейного диофантова уравнения c двумя неизвестными и построения его решения на основе алгоритма Евклида. Пусть n₁, n₂ — положительные взаимно простые числа.¹⁸ Тогда существуют такие целые m₁, m₂, что n₁m₁ + n₂m₂ = 1. На практике эти числа находят, раскладывая n₂/n₁ в цепную дробь, затем отбрасывая последнюю дробь в этой цепи, вновь всё сворачивая и получая в итоге отношение чисел m₁ и m₂ с точностью до знака.

Чтобы лучше «вспомнить всё» рассмотрим конкретный пример. Найдём целочисленные решения уравнения 14m₁ + 45m₂ = 1. Раскладываем дробь 45/14 в цепную дробь:

Выкидываем «хвост», ½, и сворачиваем всё обратно:

Ясно, что |m₁| > |m₂|, поскольку 45 > 14. Поэтому |m₁| = 16, a |m₂| = 5. Знаки чисел m₁, m₂ противоположны, и способ их определения я всегда забываю. В данном случае все рассматриваемые числа невелики, так что знаки можно определить прямым умножением. Можно также заметить, что 14·16 оканчивается цифрой 4, а 45·5 — цифрой 5, и отсюда сделать вывод, где должен стоять знак минус. В общем, m₁ = –16, m₂ = 5. После нахождения одной пары решений легко выписать и общее решение уравнения 14m₁ + 45m₂ = 1, но нам это не нужно.

Ещё один пример, чуть сложнее. Пусть n = n₁n₂ и числа n₁, n₂ взаимно просты. Тогда весь набор корней из единицы n-й степени получается перемножением наборов корней из единицы n₁-й и n₂-й степеней:

Для доказательства достаточно заметить, что оба множества содержат по n элементов и что каждый корень из единицы n-й степени можно представить в виде произведения корней из единицы меньших степеней:

С другой стороны,

а значит, если k₁ пробегает набор чисел {0, 1, ..., n₁–1}, а k₂ — набор чисел {0, 1, ..., n₂–1}, то k₁n₂ + k₂n₁ пробегает полную систему вычетов по модулю n, а k ≡ k₁n₂ + k₂n₁ (mod n) — набор чисел {0, 1, ..., n–1} (т.е. приведённую систему вычетов). Как следствие, если определить функцию

	n–1
G(n, m, l) =	∑	exp	(	πim n	k2 +	πil n	k	)	,
	k=0

где mn+l — чётное число, то, используя разложение n в произведение взаимно простых чисел n₁ и n₂, получаем

Пример в примере:

274
∑	exp	(	2πi 275	(14k2 + 45k)	)	= [· · ·] = 5√11 exp	(	–	5πi 22	)	,
k=0

и мы переходим к ещё одному вспомогательному вопросу, который нам понадобится в дальнейшем — китайской теореме об остатках.

Пусть n = n₁n₂ ... n_k и все n_k попарно взаимно просты. Тогда n_k и n/n_k взаимно просты, и можно найти целые числа m_k, M_k, удовлетворяющие равенству n_km_k + (n/n_k)M_k = 1. Китайская теорема об остатках гласит, что у системы сравнений

существует единственное решение x, лежащее в 0, n–1, и оно удовлетворяет сравнению

	K
x ≡	∑	ak	n  nk	Mk  (mod n).
	k=1

Простая иллюстрация. Пусть K = 2, т.е. n = n₁n₂. Тогда решение системы

имеет вид x ≡ a₁n₂m₂ + a₂n₁m₁ (mod n).

В частности, если n = 2·5·7·9 = 630 = n₁n₂, где n₁ = 2·7 = 14, n₂ = 5·9 = 45, то из предыдущих наработок нам уже известно, что m₁ = –16, m₂ = 5. Поэтому x ≡ a₁·45·5 – a₂·14·16 (mod 630). Для a₁ = 8, a₂ = 22, например, получаем x ≡ –3128 ≡ 22 (mod 630) и, значит, x = 22 — это единственное решение в диапазоне 0 ≤ x ≤ 629.

Теперь мы во всеоружии и готовы обсудить алгоритм Гуда–Томаса, который сводит процедуру вычисления ДПФ y = A_nx к перемножению прямоугольных матриц.

Пусть n = n₁n₂, а числа n₁, n₂ взаимно просты. Возьмём n-мерный вектор x и преобразуем его в n₁×n₂-матрицу X следующим образом:

Компоненты вектора x и элементы матрицы X специально обозначены одной и той же буквой x, чтобы подчеркнуть, что мы не делаем никаких действий над вектором x, кроме перегруппировки. Компонента x_j вектора x ставится в ту ячейку x_j1, j2 матрицы X, для которой выполняются условия  j₁ ≡ j (mod n₁),  j₂ ≡ j (mod n₂). Согласно китайской теореме об остатках, матрицу X можно наоборот преобразовать в вектор x, если элемент x_j1, j2 располагать на месте x_j, где  j ≡ j₁n₂m₂ + j₂n₁m₁ (mod n). Здесь m₁, m₂ — это те самые числа, которые удовлетворяют равенству n₁m₁ + n₂m₂ = 1.

Далее, опять же только перегруппировкой компонент вектора y, организуем матрицу Y:

В данном случае взаимно-однозначное соответствие компонент y_k ↔ y_k1,k2 строится на основе формул (25). Если задан индекс k, то пара индексов (k₁, k₂) определяется по формулам k₁ ≡ km₂ (mod n₁), k₂ ≡ km₁ (mod n₂). В свою очередь, по индексам k₁, k₂ восстанавливается индекс k:

Теперь, переходя от векторов x, y к матрицам X, Y, мы можем записать ДПФ y = A_nx в чисто матричном виде

Y = [· · ·] = A		X A		.
	n2		n1

Вот такая интерпретация ДПФ была предложена Гудом и Томасом. Мы не будем далее рассматривать вычислительные аспекты этого алгоритма (например, что он является быстрым), а перейдём, наконец, к его использованию для задачи поиска собственных векторов.

[· · ·]

Рассмотрим пример. Пусть n = 6 = n₁n₂, где n₁ = 2 и n₂ = 3. В этом случае, как мы уже знаем, spec U₂ = {+1, –1}, spec U₃ = {+1, –1, +i} и

[· · ·]

Ещё один пример, почти аналогичный. Пусть n = 12 = n₁n₂, где n₁ = 4 и n₂ = 3. В этом случае spec U₄ = {+1, –1, +i}, #(+1) = 2, #(–1) = #(+i) = 1, и

[· · ·]

И тем не менее, несмотря на всю изящность построения собственных векторов матрицы U_n по собственным векторам матриц меньшего порядка, при таком подходе к решению нам пришлось бы столкнуться с некоторыми весьма серьёзными проблемами. Во-первых, пришлось бы доказывать, что равенство Y = ±X справедливо всегда. Во-вторых, нет никакой уверенности, что мы всегда будем получать линейно независимые векторы. А вдруг при n = 28, построив семь собственных векторов для λ = 1 и достраивая последний восьмой, мы получим в итоге линейную комбинацию предыдущей семёрки? В-третьих, даже зная разложение n в произведение простых чисел, n = ∏ p_k^αk, нам не обойтись без построения собственных векторов матриц типа U_pα, причём непонятно из каких соображений.¹⁹ Поэтому хорошо было бы предложить более основательный подход, который не нуждается в необходимости предварительной «разборки» n на множители p^α и последующей «сборки» векторов h⁽ⁿ⁾ из компонент h^(pα).

Такой подход действительно существует. Аналитический в своей основе, он был опубликован Владимиром Матвеевым,²⁰ и к его изложению мы сейчас и приступаем.

Начнём, как всегда, издалека. Пусть

— интегральное преобразование Фурье. Спрашивается, как построить функцию  f (x), инвариантную к этому преобразованию?

Для начала заметим, что для гладких быстро убывающих функций  f (x)

[· · ·]

Конечно, продемонстрированный подход к решению задачи может показаться разочаровывающе скучным, поскольку он ничего не говорит о способе получения таких красивых фурье-инвариантов, как

и всё же у него есть одна важная деталь: он связывает возможность построения собственных функций оператора F, соответствующих собственному числу λ = 1, со степенями оператора применительно к произвольной функции  f (x).

Можно предложить другой, не такой абстрактный способ построения функций, инвариантных к преобразованию Фурье. Если вспомнить интеграл

∫	e–at² + ibt dt =	Ö	π  a	e–b²/4a,       (Re a > 0,   Re(√a) > 0),
R

который был найден при обсуждении формулы Пуассона, то, полагая a = ½ и b = –x, мы получаем равенство

и самый известный фурье-инвариант — функцию e^–x²/2 — собственную функцию оператора F, соответствующую собственному числу λ = 1. Дифференцируя (33) по x, находим ещё одну собственную функцию:

теперь для λ = –i. На этом лёгкая жизнь заканчивается — ещё одно дифференцирование приводит к равенству

и для нахождения следующей собственной функции мы вынуждены устраивать линейные комбинации с предыдущими наработками:

Получается функция e^–x²/2(x² – ½) — собственная функция оператора Фурье, соответствующая собственному числу λ = –1.

Понятно, что таким способом (дифференцирование + линейные комбинации) можно двигаться и дальше. Ещё понятно, что такое черепашье движение быстро надоест. Чтобы одним махом построить всё семейство собственных функций, проделаем следующее.

Пусть  f (x) — бесконечно дифференцируемая на R функция, такая что

(иначе говоря,  f (x) принадлежит пространству Шварца). Возьмём преобразование Фурье от  f (x) и продифференцируем его по x:

С другой стороны,

Поэтому

Далее,

Таким образом, если F[ f ](x) =  f (x), т.е.  f (x) — собственная функция оператора Фурье, соответствующая собственному числу λ=1, то (x – d/dx)ⁿ f (x) — собственная функция, соответствующая числу λ=(–i)ⁿ. Выбирая  f (x) = e^–x²/2, мы получаем семейство функций

где H_n(x) — некоторый многочлен степени n. В частности, H₀(x) = 1, H₁(x) = 2x, H₂(x) = 4x² – 2. В общем случае,

Многочлены H_n(x) называют многочленами Эрмита, а последнее равенство — формулой Родрига. Здесь дифференциальный оператор более прост, поэтому формула Родрига часто используется как отправная точка для получения различных свойств многочленов Эрмита. Например, производящей функции, значений в нуле и разложения по степеням x.

Самым важным свойством функций ψ_n(x), наряду с инвариантностью к преобразованию Фурье, является тот факт, что множество {ψ_n(x), n = 0, 1, ...} образует базис в пространстве L₂(R), причём ортогональный:²²

(ψm, ψn) =   ∫  ψm(x)ψn(x) dx =   ∫  e–x²Hm(x)Hn(x) dx = √π 2nn! δn,m. R R

(35)

Чтобы доказать равенство (35), достаточно построить для скалярных произведений производящую функцию и сравнить коэффициенты при s^mtⁿ в исходном и конечном рядах.

Н-да... Что-то я чересчур увлёкся, и первоначальная задумка кратенько обсудить оператор Фурье F и его собственные функции, а затем вернуться к собственным векторам матрицы U обернулась длинным разговором. Постараюсь быть менее многословным.

Итак, собственные функции ψ_n(x) оператора Фурье образуют ортогональный базис в L₂(R). Поскольку F [ψ_n](x) = (–i)ⁿψ_n(x), то spec F = {±1, ±i}, а само пространство L₂(R) раскладывается в прямую сумму ортогональных подпространств, соответствующих собственным числам оператора F :

L2(R) = S+1 + S–1 + S+i + S–i =   ∑  Sλ. λ Î spec F

(36)

Для каждого собственного числа λ определим проекционный оператор (или проектор) P_λ: L₂(R) → S_λ. Тогда разложение (36) можно представить в операторном виде:

т.е. тождественный оператор — это сумма проекторов на ортогональные подпространства.²³ Наша цель теперь — выразить эти проекторы через оператор Фурье.

Возьмём произвольную функцию  f (x) Î L₂(R) и представим её в виде суммы четырёх функций из ортогональных подпространств:

Тогда

F f  =	∑	F fλ =	∑	λ fλ =	∑	λPλ f,
	λ Î spec F		λ Î spec F		λ Î spec F

или

Используя очевидные свойства проекторов, P_λ² = P_λ и P_λ1P_λ2 = 0 для λ₁ ≠ λ₂, находим разложения степеней оператора Фурье:

Четвёртая степень оператора Фурье — это тождественный оператор, поэтому достаточно оставить только равенства для степеней k = 1, 2, 3, 4 и рассматривать получившиеся соотношения как систему линейных уравнений для нахождения проекторов:

Вот теперь мы можем взять любую функцию  f (x) Î L₂(R) и построить из неё собственные функции оператора Фурье для каждого λ, просто навесив на  f  проектор P_λ и посчитав всего один интеграл, поскольку F ²[ f ](x) =  f (–x):

Как видно, проекторы не только разбивают пространство L₂(R) на четыре ортогональных подпространства, но проделывают это таким образом, что

Перед тем как покончить с обсуждением оператора Фурье и вернуться из бесконечномерного L₂(R) в конечномерное Cⁿ (где, как вы уже догадываетесь, мы тоже будем строить проекторы на собственные подпространства, но только теперь для матрицы U) рассмотрим один пример. Возьмём в L₂(R) базис {e^–ax²xⁿ; Re a > 0, n Î N} и посмотрим, что из него сделают наши проекторы P_λ.

[· · ·]

Окончательные выражения довольно громоздки и только при a = ½ немного упрощаются. Вытекающие отсюда разнообразные следствия относительно многочленов Эрмита (интегралы, разложения в ряды и конечные суммы) мы не станем обсуждать.

В каком-то смысле, долгое исследование положения дел с оператором Фурье, действующим в пространстве L₂(R), оказалось выгодным. Многие моменты переносятся на матрицу U, действующую в пространстве Cⁿ, без всяких изменений. Причина кроется в унитарности,

и совпадении спектров,

Не нужно даже проверять, что U ⁴ = E, поскольку спектральное разложение оператора U сделает всё за нас. Тем не менее повторим некоторые утверждения (практически, в тех же обозначениях!). Раскладываем Cⁿ в прямую сумму ортогональных подпространств S_λ:

или, что то же самое, раскладываем единичную матрицу в сумму проекторов:

Поскольку все собственные векторы матрицы U вещественны, то каждый проектор P_λ — это n×n-матрица с вещественными, а не комплексными элементами. Раскладываем по проекторам матрицу U и её степени:

По ходу дела получаем, что U ³ = U и U ⁴ = E. Записываем систему и её решение:

Теперь небольшая перегруппировка слагаемых в правой части, учитывающая вещественность матрицы U ², завершающее добавление вектора v и всё — дело, в принципе, сделано! P₊₁v&=\dfrac{1}{2} \biggl(\dfrac{v+U²v}{2}+ \dfrac{Uv+\bar Uv}{2}\biggr),&     P_+iv&=\dfrac{1}{2} \biggl(\dfrac{v-U²v}{2}-i \dfrac{Uv-\bar Uv}{2}\biggr),\\ P_–1v&=\dfrac{1}{2} \biggl(\dfrac{v+U²v}{2}- \dfrac{Uv+\bar Uv}{2}\biggr),& P_–iv&=\dfrac{1}{2} \biggl(\dfrac{v-U²v}{2}+i \dfrac{Uv-\bar Uv}{2}\biggr).\\ \end{array} \eqno(38)

Осталось выбрать какой-нибудь определённый вектор v и проектор P_λ превратит его (если только не занулит) в собственный вектор матрицы U, соответствующий собственному числу λ.

А теперь построение полной системы собственных векторов матрицы U. Берём набор базисных ортов и напускаем на этот стандартный набор проекторы P_λ.

[· · ·]

Осталось показать, что получившийся набор из n векторов действительно является базисом (т.е. набором ЛНЗ векторов), и это последнее задание для вас на сегодня. :)

Рассмотренные выше вопросы — лишь самое начало большой и многогранной теории. После суммы Гаусса ∑ζ^m² можно перейти к суммам корней из единицы в более высоких степенях, например, ∑ζ^m³, заняться теорией характеров и высшими законами взаимности. А можно вернуться к одномерным функциям ψ_n(x) = e^–x²/2H_n(x) и заняться их переносом на N-мерный случай. Например, при N = 2 моды Эрмита–Гаусса

обретают весьма достойных родственников — моды Лагерра–Гаусса

где L_n^m(t) — многочлены Лагерра. Оба семейства функций ведут своё происхождение от одной гауссовой функции: H_0,0(x, y) = L_0,0(x, y) = e^{–(x² + y²)/2}. Оба семейства являются ортогональными базисами в L₂(R²) и инвариантны к двумерному преобразованию Фурье:

И, подобно тому как суммы Гаусса, функции Гаусса и ряды квадратичных экспонент тесно переплетаются и плавно трансформируются друг в друга (заставляя задуматься о почти философских вопросах взаимосвязи дискретного и непрерывного), оба семейства связаны между собой и алгебраическими, и интегральными соотношениями:

где P_n^(μ,ν)(t) — многочлены Якоби.

Гильберт, не будучи очень высокого мнения о способностях среднего студента, считал, что ничего нельзя усвоить, пока не услышишь несколько раз.

— Пять раз, Герман, пять раз! — говорил он Герману Вейлю, когда тот начинал свою педагогическую деятельность. — Вычисления проводи не выше, чем на уровне таблицы умножения, и начинай всё с простейших примеров.

Более точно, для любого отрезка [τ₁, τ₂], вложенного в интервал (t₁(μ₀), t₂(μ₀)), на котором существует решение u(t, μ₀), найдётся такая окрестность U(μ₀) = {μ Î R^m: |μ – μ₀| < δ}, что

1) для любого μ Î U(μ₀) решение u(t, μ) определено при всех t Î [τ₁, τ₂];

2) решение u(t, μ) является непрерывной функцией на [τ₁, τ₂] ×U(μ₀).

Такая замена чем-то сродни выделению полного квадрата для уравнения ax² + bx + c = 0, т.е. шагу, который позволяет решить уравнение 2-й степени. Для уравнений 3-й и 4-й степеней линейная замена тоже полезна, хотя, конечно, её роль уже не столь велика. О решении алгебраических уравнений (и не только) можно прочесть в книге: В. В. Прасолов, Ю. П. Соловьёв. «Эллиптические функции и алгебраические уравнения» (М., Факториал, 1997), в которой, в частности, приведены три способа решения уравнения 4-й степени. Не могу удержаться, чтобы не привести ещё один, вычитанный в заметке В. М. Дубровского «Об уравнениях четвёртой степени» (УМН, 1973, № 4).

Рассмотрим общее уравнение 4-й степени

где e — не обязательно основание натуральных логарифмов, хотя может и быть им. :) Пусть α* = –b/4a, т.е. точка, в которой P'''(x) = 0. Введём обозначения P = P(α*), P' = P'(α*), P'' = P''(α*). Пусть H — какой-либо корень уравнения

отличный от нуля. Тогда корни исходного уравнения имеют вид

x1,2 = α* + √H ±	√	–H –  P'' 4a  –  P' 4a√H	,     x3,4 = α* – √H ±	√	–H –  P'' 4a  +  P' 4a√H	.

И, уж совсем уходя в сторону, несколько слов об уравнениях 5-й степени (см. ещё одну заметку В. М. Дубровского, опубликованную в УМН, 1973, № 5). Как известно (теорема Абеля), уравнение общего вида

в радикалах не решается. Пусть, однако, кубические уравнения

имеют общий корень, который мы обозначим через α*. Тогда четыре корня уравнения P(x) = 0 имеют вид x_1,2 = α* ± √H₁, x_3,4 = α* ± √H₂, где H₁, H₂ — корни квадратного уравнения

Если A — произвольная матрица n×n, то

Замена λ → –λ:

Перемножаем обе формулы:

Ещё одна замена λ² → λ: