Псевдоспектры матриц: высокоточное построение на базе Python

Спектральные портреты конечномерных линейных операторов (матриц) представляют весьма удобный инструмент анализа устойчивости их собственных значений к возмущениям. Такой анализ может быть полезен, например, в случаях, когда необходимо решить, стоит ли доверять полученным численно собственным значениям или нет, или решить какую-либо прикладную задачу, связанную с исследованием устойчивости поведения динамической системы. Хорошо известно, что собственные значения матриц могут существенно изменяться под действием даже незначительных возмущений; однако насколько эти возмущения могут быть велики, отлично иллюстрирует следующий пример с матрицей, взятой из галереи тестовых матриц MatLab: gallery(5). Эта матрица 5х5 следующего вида: $$
A=\left(\begin{array}{lllll}
-9& 11& -21& 63& -252\\
70& -69& 141& -421& 1684\\
-575& 575& -1149& 3451& -13801\\
3891& -3891& 7782& -23345& 93365\\
1024& -1024& 2048& -6144& 24572
\end{array}\right)
$$

имеет 5-кратное нулевое собственное значение! 

Тем не менее, стандартная функция MatLab для вычисления собственных значений  eigs этого не обнаруживает, как и не обнаруживает того факта функция библиотеки NumPy eigvals:  np.linalg.eigvals; в обоих случаях – результат – ненулевые комплексные собственные значения.

Для NumPy 1.9.2 результат вычислений собственных значений этой матрицы (с удвоенной точностью) такой:
 

>np.linalg.eigvals(A)
array([-0.03697710+0.02749668j, -0.03697710-0.02749668j,
        0.01470225+0.04265918j,  0.01470225-0.04265918j,  0.04454971+0.j        ])

Интересно и то, что интсрументарий MatLab Eigtool , подтверждает неправильные численные собственные значения:
 

Eigtool Gallery(5) failure?!

Таким образом, построение спектрального портрета для матрицы A  в программе Eigtool подтверждает "правильность" ненулевых собственных значений, полученных при помощи функции eigs: вокруг "точных" собственных значений образовались кольца псевдоспектра... уже и сам, пока писал эти строки, стал сомневаться в том, что все собственные значения матрицы нулевые! Однако, символьные расчеты убеждают в обратном – собственные числа матрицы A нулевые (их можно провести, например, в Maple).

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

Было решено использовать самый простой алгоритм, основанный на определении спектрального портрета и вычислять значения \(\varepsilon\) для каждой точки регулярной сетки, покрывающей фрагмент комплексной плоскости. Первый вариант решения этой задачи был (когда-то) представлен на некогда моем сайте pseudosp.narod.ru (теперь управление им утеряно, и я рад бы удалить его, да не получается:): это был MatLab скрипт, использующий движок символьных вычислений Maple для вычисления сингулярных значений матриц с повышенной точностью. Помню, это был 2007 или 2008 г, вычисления на моем компьютере на сетке в 1 миллион точек (при точности не менее 100 знаков) осуществлялись около суток, и, результат был существенно отличен от результата работы Eigtool. Спектральный портрет был весьма изрезан и не концентрировался вокруг неверных собственных значений: по крайней мере, такой спектральный портрет не убеждал в правильности ненулевых собственных значений, что уже было хорошо.

Со временем я отошел от использования MatLab в сторону Python+SciPy+Matplotlib. В качестве "спортивного интереса" это побудило реализовать функцию расчета спектрального портрета матрицы на Python. Поскольку сеточный алгоритм нахождения спектрального портрета прекрасно подходит для параллельной реализации было решено использовать стандартный модуль multiprocessing,  а для осуществления вычислений с произвольной точностью – модуль mpmath.

Из всего этого вышел следующий модуль: mpseudo, позволяющий осуществлять высокоточные параллельные расчеты псевдоспектров матриц. С его помощью я построил спектральный портрет матрицы A (с использованием 100 значащих цифр при расчете сингулярных значений). Он намного правдоподобней, чем тот, что получен с помощью Eigtool (а все из-за точности вычислений):

Pseudospectrum of gallery(5)

Из последнего спектрального портрета следует, что даже незначительные относительные возмущения матрицы (сравнимые с удвоенной точностью представления вещественных чисел в ЭВМ \(10^{-16}\) приводят к существенным искажениям ее собственных чисел; таким образом, причина ненулевых собственных значений при попытке их вычисления с удвоенной точностью становится вполне ожидаемой: "красное пятно" псевдоспектра заполняет весь комплексный квадрат [-0.05, 0.05], он соответствует значению \(\varepsilon\approx 10^{-16}\), а это означает, что собственные значения при уровне относительных возмущений \(\varepsilon\) могут оказаться всюду в этом квадрате.

blog comments powered by Disqus