Исследование способов численного нахождения собственных значений и векторов матрицы.

Задача лабораторной работы заключалась в написании и исследовании работы трёх методов нахождения собственных значений и векторов матрицы: степенной метод, метод Данилевского и QR-алгоритм. При написании кода я опирался на конспект лекций по вычислительным методам алгебры Фалейчика Бориса Викторовича.

Степенной метод

Вариации метода

Я реализовал три вариации степенного метода для нахождения собственных чисел и векторов матрицы: простой степенной метод для нахождения одного наибольшего собственного числа и соответствующего ему вектора, более сложную реализацию для случая, когда у матрицы два наибольших по модулю собственных значения совпадают по модулю и более простая реализация не сходится, а также еще более изощренный метод для случая, если наибольшими по модулю собственными числами является пара комплексно-сопряженных чисел. Асимптотики методов: для всех алгоритмов ( - количество итераций). Приведу краткое обоснование:

1. В первой модификации на каждой итерации мы умножаем матрицу на вектор, что и дает асимптотику на итерацию с учетом, что остальные операции над векторами производятся за линейное время.

2. Во второй модификации достаточно просто в два раза больше операций, чем в первой. Асимптотика та же.

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

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

Критерий остановки

Критерием остановки итераций я выбрал , так как в степенном методе сходится именно последовательность . В некоторых местах также используется невязка .

Выбор лучшей модификации

Хоть модификация метода для комплексных значений и покрывает все случаи, скорость ее работы меньше скорости двух первых модификаций и, соответственно, можно взять лучшее от трех модификаций если выбрать наиболее подходящий для данной матрицы метод. Для этого я пробовал разные способы определения алгоритма:

• Для этого я вынес итерации методов в отдельные функции и написал алгоритмы, проверяющие схождение методов. Алгоритм заключается в следующем: фиксируем количество проверочных итераций, пусть это будет ; далее запускаем итераций проверяемого метода, после чего следующие итераций проверяем, чтобы модуль разности уменьшался каждую итерацию. Если это так, то запускаем проверяемую на сходимость модификацию метода уже до соответствия условию невязки, используя значения, полученные за проверочных итераций, чтобы не пересчитывать уже полученные данные еще раз. Проверку методов начинаю с самого простого с точки зрения вычислительной сложности, то есть с самой простой модификации, далее вторая модификация, а если и она не сходится, то мы имеем дело с комплексными значениями, и запускается третий метод.

• Но оказалось, что такой метод часто дает ложноотрицательное отсутствие сходимости, поэтому я попробовал другой: снова зафиксировал и теперь еще и если , то метод сходится. Это тоже не дало хороших результатов. Тогда я решил для начала потестировать алгоритмы.

Зависимость скорости работы от размерности матрицы

Я провёл несколько разных тестов. Результаты меня удивили, так как судя по асимптотике (с учетом большой константы третьего метода) было вполне понятно, каких результатов ждать: первый алгоритм самый быстрый, за ним второй и третий.
Для начала методика тестирования. Так как степенной метод является итерационным и не на всех матрицах сходится за разумное время, встал вопрос о методике тестирования: ведь нет большого смысла тестировать алгоритм и засекать время его работы для матриц, на которых он имеет крайне низкую скорость сходимости. Я решил тестировать алгоритмы только на матрицах, на которых они сходятся за разумное количество итераций.
Первым тестом был тест времени работы алгоритмов и скорости их сходимости (кол-во итераций). Для этого я фиксировал размер матрицы, генерировал случайные матрицы (вся генерация в лабораторной работе осуществлялась с помощью и ) и если за разумное (фиксированное) кол-во итераций метод (какой?) сходился, то я добавлял эту матрицу в список пригодных для тестирования (я генерировал по 5-10 матриц для каждой размерности для более точного определения времени работы). Процесс это небыстрый (дальше будут тесты сходимости алгоритмов), поэтому я решил использовать все возможности своего компьютера и написал процесс перебора матриц многопоточно, что конечно дало отличный множитель скорости. На этом моменте я чуть не допустил фатальную ошибку: передавал всем потокам один для генератора. Осталось разобраться только с тем, какой алгоритм брать за эталонный (по скорости сходимости которого будем судить брать ли матрицу в список пригодных). Очевидно, если матрица сходится для первого метода, то она сойдется для всех остальных (возможно, с бОльшим кол-вом итераций, но всё же за разумное время). Как и если матрица сошлась для второго алгоритма, то она сойдется и для второго. Я провел два теста: в одном алгоритмом для выбора матрицы послужил первый, во втором - второй.
На графиках ниже результаты первого и второго теста.

Результаты, ожидаемо, идентичные, но без второго теста тестирование не было бы полным.
Третий алгоритм ожидаемо работает сравнительно долго, а вот первый и второй идут бок о бок. Я это объясняю тем, что хоть у второго алгоритма выше вычислительная сложность, он делает на почти в два раза меньше итераций, чем первый и почти в три меньше, чем третий (графики решил не вставлять, так как из-за предвыбора матриц график тут рисовать неуместно. Тенденция с разницей в количестве итераций прослеживается на всех размерностях). Теперь становится понятно, как выбирать алгоритм для того самого режима.
Стоит отметить, что границы генерации чисел в этом тесте не влияют на тенденции. Прослеживается другая зависимость, но это уже следующее исследование.

Общее исследование сходимости алгоритма

Для проверки сходимости алгоритма (в том числе на разных размерностях матриц и границах генерации чисел) я провёл следующее исследование. Я построил диаграммы распределения числа сделанных итераций для трех алгоритмов для 22000 матриц


Тут слишком долго тестировалось, пришлось снизить количество матриц до 2200.


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

Выбор лучшей модификации 2

Теперь ясно, что стоит выбирать один из двух алгоритмов: быстрый или медленный , который сходится для комплексных чисел. Первый алгоритм работает примерно за то же время, что и второй, обладая худшей сходимостью, так что даже не стоит тратить время на проверку сойдется ли матрица для него. Остается проверить, сходится ли второй алгоритм.
Я придумал просто запускать алгоритм с намного меньшей точностью и если он сошелся для этой, то предполагается, что сойдется и для исходной. Если же не сошелся для одной из них, то запускается третий метод. Это, кстати, тоже не сработало.
Сработал следующий алгоритм выбора: сделаем какое-то фиксированное количество итераций второго алгоритма, а затем проверим невязки . Если невязка достаточно мала, предполагается, что алгоритм сходится.
Вектор (им описывается всё состояние алгоритма) можно переиспользовать, чтобы итоговый алгоритм сделал меньше итераций, я это также реализовал.
Приведу графики скорости работы гибридного алгоритма на матрицах, сходящихся для второго алгоритма, а также на матрицах, сходящихся для третьего (время выполнения на разных графиках сравнивать некорректно, так как использовались разные наборы для тестирования).

Оптимизации

Теперь, когда собрано много информации о работе алгоритмов, можно подумать об оптимизации. Так как в третьем алгоритме много времени занимает вычисление собственных значений из вектора , а итераций алгоритм делает очень много, то есть смысл делать какое-то количество итераций без вычисления собственных значений и векторов. Хоть это и не асимптотическое улучшение, оно приносит заметные результаты. Я подобрал количество итераций для предварительного расчета, равное . График 9 показывает успех оптимизации:

Метод Данилевсвого и нахождение корней многочлена

Метод Данилевского. Общие сведения

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

Преобразование к форме Фробениуса

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

Ожидаемый .

Метод Данилевского

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

Извиняюсь за обрезанный текст, почему-то решил, что тут он именно так его поместит.

Нахождение корней многочлена. Метод Ньютона и бисекции

Итак, у нас есть многочлен и нам нужно найти его корни. Для начала я реализовал метод, который основывается на том, что можно точно посчитать производную многочлена. В такому случае можно найти экстремумы многочлена как корни производной (ее корни находим так же - рекурсивно), а искомые корни многочлена будут лежать между экстремумами, будем находить их бисекцией. Нужно не забыть сокращать многочлен на порядок производной чтобы коэффициенты не становились слишком большими. Звучит отлично, но на практике я не был впечатлен результатами. Я придумал другой метод: будем запускать метод Ньютона от какой-то очень дальней точки, он будет сходится к какому-то корню, который мы будем уточнять методом бисекции. Производную для метода Ньютона посчитаем аналитически, так как это будет точнее, чем считать ее численно (еще и быстрее). Затем поделим многочлен на и будем повторять это пока не дойдем до многочлена первой степени. Достаточно большую границу будем брать по формуле: \max(\frac{|a_n|}{b + |a_n|}, 1 + \frac{c}{|a_0|}) b = \max(|a_0|, |a_1|, ..., |a_{n - 1}|) c = \max(|a_1|, |a_2|, ..., |a_n|)

Нахождение собственных векторов

Во-первых, зная собственное число, можно рассмотреть слау . Так как - собственное число матрицы, то матрица будет иметь нулевой определитель, следовательно однородная слау будет иметь ненулевое решение - искомый вектор. Такое решение будет находить каждый собственный вектор за , что не очень-то хорошо с учётом того, что таких векторов тоже может быть . Я думал над каким-либо треугольным разложением для такой матрицы, но всё портит.
Также я нашел способ быстрее. Можно заметить, что для ( - матрица формы Фробениуса) собственный вектор для значения считается как степени начиная с правой стороны (легко показать, если последовательно решить СЛАУ, то есть занулить всю первую строку, а затем положить последний элемент искомого вектора равным единице). Тогда остается только привести его к вектору изначальной матрицы. Для этого запомним преобразования столбцов, которые мы делали для приведения матрицы к форме Фробениуса и выполним их для этого вектора в обратном порядке (с оговоркой, что при прибавлении столбца к будем наоборот прибавлять к ). Количество преобразований, очевидно, ограничено , тогда асимптотика алгоритма будет ( - количество собственных значений).

QR Алгоритм

Снова-таки алгоритм написан ровно по Фалейчику: приведение к форме Хессенберга, а затем итерации вращения. Асимптотика алгоритма: . Критерий остановки такой: делаем итерации пока на диагонали есть пересекающиеся квадраты, либо количество нулей под главной диагональю . Если оба этих критерия перестали выполняться, то итерации можно останавливать. Далее извлекаем собственные значения, тут совсем просто. Можно найти собственные векторы способом, описанным выше, за . Или же привести к форме Фробениуса и найти за .
Из проблем в ходе реализации можно перечислить все деления в алгоритме: и в отражениях, и в поворотах. Во всех местах, где есть деление, стоит отдельно рассматривать случай нулевого знаменателя (обычно это означает, что нужный нам элемент(ы) уже нулевой).
График времени работы алгоритма (методика тестирования такая же как в задаче 1).

Сходимость

Снова-таки проведем такие же тесты как для задачи 1.


График красивый, а ситуация страшная:

Можно сделать интересный вывод о том, что QR-алгоритм работает лучше степенного метода на матрицах с большими по модулю числами, но хуже для матриц больших размерностей.