О вычислении матричной экспоненты

    image

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

    1. Вычисление матричной экспоненты


    Вычисление экспоненты
    image
    где E — единичная матрица, t — время, связано с необходимостью расчета высоких степеней матрица A. Получим формулу, позволяющую определить матричную экспоненту с помощью n степеней матрицы A (n — ее порядок).

    Пусть характеристическое уравнение матрицы A имеет вид
    image
    По теореме Гамильтона-Кэли [2] матрица A удовлетворяет матричному уравнению, аналогичному (2):
    image
    откуда
    image
    Следуя методу Д.К. Фаддеева [2], коэффициенты характеристического уравнения определяются по рекуррентному соотношению
    image
    где image — след матрицы image (сумма элементов, стоящих на главной диагонали), image

    Далее введем обозначение: если m=0, то image; иначе (при натуральном m)
    image

    Умножим обе части соотношения (3) на матрицу A с учетом введенных обозначений. Получим
    image
    Выражение (5) можно переписать как
    image
    Теперь умножим обе части равенства (6) на матрицу A, подставив при этом в полученное соотношение формулу (3):
    image
    Тогда из выражения (7) с помощью последовательного умножения на матрицу A обеих его частей следует, что
    image

    Теперь представим матричную экспоненту как
    image
    Откуда имеем
    image

    2. Описание алгоритма


    Для реализации вычисления матричной экспоненты, согласно (8), был применен следующий алгоритм. Сначала нужно инициировать результат значением нулевой матрицы. Вычислить image для k от 0 до n. Далее выполнить для k от 0 до n–1 следующую последовательность операций:
    1. Вычислить сумму image. В качестве критерия прекращения суммирования использовать условие image, где image — положительное число, характеризующее точность вычисления суммы.
    2. Используя значения, полученные ранее, определить произведение image и прибавить его к текущему значению результата.

    При вычислении матричной экспоненты с помощью данного алгоритма используется рекуррентное соотношение (4). При больших значениях m и k большинство значений q будут рассчитываться повторно много раз. Поскольку q(m,k) является чистой функцией (зависит только от входных аргументов), то будет разумно применить стратегию мемоизации.

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

    3. Сравнение с классическим алгоритмом


    Разработанный алгоритм обеспечивает вычисление матричной экспоненты, используя только n степеней матрицы, в то время как классический алгоритм (1) не детерминирован, и вычисления степеней матрицы продолжаются, пока не выполнится условие останова алгоритма. Обычно таким условием является

    image

    Заметим, что классический алгоритм имеет потребление памяти image. Описанный алгоритм, в виду необходимости хранить n-1 степень матрицы A, имеет потребление памяти image.

    Нами был проведен вычислительный эксперимент с целью сравнить быстродействие алгоритмов. Для этого была разработана KipDblK программа из комплекса [3] на языке C++, реализующая оба алгоритма. С помощью данной программы были произведены расчеты матричной экспоненты для матриц различного размера. Порядок матрицы изменялся от 2 до 132. Матрица инициализировалась случайными числами в диапазоне [0;1]. Экспонента вычислялась для t=1. Результаты сравнительного эксперимента представлены на рис. 1. По оси абсцисс отложен порядок матрицы A, по оси ординат — время счета. Полученные точки соединены сплайнами для наглядности.

    image

    Рис. 1. Сравнение временных характеристик классического алгоритма (верхняя кривая) вычисления матричной экспоненты и алгоритма, описанного в данном топике (нижняя кривая).

    P.S.


    Данный топик был подготовлен по материалам нашей статьи [4].

    Литература


    1. Демидович Б.П. Лекции по математической теории устойчивости. — М.: Наука, 1967.
    2. Гантмахер Ф.Р. Теория матриц. – М.: Наука, 1967.
    3. KipDblK maxima_comm.tar.gz.
    4. Безгин С.В., Пчелинцев А.Н. Организация матричных и символьных вычислений для исследования поведения решений обыкновенных дифференциальных уравнений // Системы управления и информационные технологии, 2012. Т. 47, №1. — С. 4-7.

    Similar posts

    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 47

    • UFO just landed and posted this here
    • UFO just landed and posted this here
        0
        В случае диагональных матриц матричная экспонента состоит из скалярных экспоненциальных функций. В формуле (8) коэффициенты r будут неявно определять часть их разложения в степенной ряд. Когда значения показателей степени в таких функциях достаточно велики, то наш алгоритм проиграет — слишком много членов разложения в ряд придется брать, здесь лучше оптимизировать вычисление скалярных экспонент, выделив целую и дробную части в показателе.
        +1
        Порядок матрицы изменялся от 2 до 100.

        А на графике порядок по оси абсцисс продолжается до 140, а точки следуют, судя по всему, до 130-го порядка. Это опечатка или результат экстраполяции?
          +1
          Спасибо за замечание. Это опечатка. Поправил.
          0
          В подписи к рисунку написано, что классический алгоритм — это верхняя кривая, то есть с крестиками, а на самом рисунке — что это кривая со штрихами, то есть нижняя. В тексте статьи об этом не написано ничего. Так что же показало сравнение быстродействия алгоритмов?
          • UFO just landed and posted this here
              0
              Одно из основных следствий из закона Мэрфи гласит, что, если что-то может быть не понято, то оно будет непонято.
              Когда я ещё был молодым учёным и специалистом, умные дяди меня учили, что именно поэтому крестики — очень плохой выбор. А главное — обсуждение каждого графика и вывод по каждому графику должны быть в тексте.
            0
            Честно говоря, неубедительный прирост скорости, с учётом возросших на порядок требований к памяти.
              0
              Хм, выходит никто не считает e^At с помощью жордановой формы?
                0
                Жорданова форма вообще малопригодна для вычислений, поскольку неустойчива.
                  0
                  Вычисляют, но там, как мне известно, сначала нужно найти собственные числа матрицы A.
                  +1
                  Мне на почту пришло сообщение о том, что в пакете MATLAB вычисление матричной экспоненты для матрицы 140-ого порядка, заполненной случайными числами в интервале от (0;1), занимает 0.005371 секунд (пакет выдал такое время). Используемые команды

                  tic; expm(rand(140)); toc
                  

                  В описании команды expm указано, что пакет использует модификацию аппроксимации Паде, рассмотренную в работе Nicholas J. Higham «The Scaling and Squaring Method for the Matrix Exponential Revisited». Наверное, такая скорость вполне может возможна. Я попробовал вычислить матричную экспоненту в WolframAlpha. Например, для матрицы 60x60

                  MatrixExp(RandomReal[{0,1},{60,60}])
                  

                  вычисления идут несколько секунд (прямая ссылка). Для матрицы 140-ого порядка система не выдала результата.
                  • UFO just landed and posted this here
                    0
                    Простенькая формула Exp(A)=Exp(A/(2^N))^(2^N) позволила ускорить поиск экспоненты матрицы 500*500 в 3.4 раза (с 16.2 до 4.8 сек на C# безо всяких оптимизаций). Относительная погрешность результата — лучше 10^-13.
                      0
                      Такой прием можно и здесь применять для ускорения сходимости рядов.
                      0
                      Совсем не впечатлило ускорение, вот например в scipy функция expm3 (считает экспоненту через разложение в ряд) вычисляет для матрицы 150х150 всего за 13 мс (и это на достаточно старом ноуте). Хотя я сам реализацию её не смотрел, но в документации написано именно про наивный алгоритм.
                        0
                        А правильный результат он дает на таких размерах?
                          0
                          Каюсь, не проверил — действительно, по-умолчанию даёт весьма неточный результат. Однако, можно выбирать любое число членов ряда (странно, что автоматически результат плохой — видимо, из-за того, что это не основной метод в библиотеке для этого), и в данном случае 200 членов с большим запасом хватает — относительная точность порядка 1e-15. Считается же теперь дольше, но намного быстрее вашего варианта — 150 мс.
                          Для сравнения, обычный expm оттуда же считает около 15 мс и точно — он как раз использует ту же аппроксимацию, что и вышеупомянутый матлаб.
                            0
                            + наверное, формула, указанная здесь habrahabr.ru/post/233587/#comment_8036405 Точность берется по норме общего члена ряда?
                              0
                              Судя по исходникам scipy, там просто наивный алгоритм суммирования степеней, и всё (да и вроде по времени сходится — матричное умножение A.dot(A) для такой матрицы у меня занимает около 0.6 мс, что при повторении 200 раз будет 120 мс). Про точность вопрос не понял, я сравнивал с expm и expm2 — они используют различные методы, но все три функции выдают практически совпадающие результаты.
                                0
                                Про точность см. неравенство в пункте 3 этого топика. Похоже, что при суммировании там берется фиксированное количество членов ряда.
                                  0
                                  А, ну да, я же говорю — указывается число членов вручную, автоматический выбор места остановки там не реализован. Просто странно, что у вас наивный метод работает на два порядка дольше.
                                    0
                                    По-моему, ничего странного, если алгоритм не детерминирован. Я не зря сделал сравнение с временем у системы WolframAlpha.
                                      0
                                      Всмысле, не детерминирован? Вроде как все шаги вполне однозначны, значит алгоритм детерминирован. Да и как ни крути, в данном случае отличия на два порядка ну никак не должно быть. Сколько у вас членов суммируется в итоге?
                                      С вольфрамальфа как раз зря сравнение сделали — кто знает, с какими приоритетами и на каких мощностях там запускаются такие вычисления, тем более бесплатно… Он же совсем не для этого заточен. Нужно сравнивать с локальными программами/библиотеками.
                                        0
                                        Заканчивает свое исполнение при выполнении неравенства, при этом мы заранее не знаем, сколько итераций будет сделано.
                                          0
                                          Не понял, к чему именно этот ответ? Если к вопросу Сколько у вас членов суммируется в итоге?, то несмотря на неизвестность числа членов заранее, постфактум это всё же можно подсчитать.
                                            0
                                            Да, оценив заранее норму матрицы.
                                              0
                                              Ну так, сколько же у вас шагов выполняется?
                                              Просто как-то складывается впечатление, что взяли заведомо неоптимальную реализацию наивного алгоритма, и показываете, что ваш работает быстрее. Кстати, сравнения точности в статье не увидел.
                                                0
                                                Точность вычисления матричной экспоненты была взята равной 0.00001. Я не задавался целью определить количество членов суммы, просто запускается цикл суммирования до выполнения неравенства в пункте 3. При заданном числе слагаемых вряд ли можно попасть в заданную точность по норме текущего члена матричного ряда. Если интересует количество, то приведу следующую оценку. Предположим, матрица A имеет размер 100x100 и заполнена одинаковыми числами, равными 0.5 (как указано в топике, я беру из отрезка [0;1]). Тогда норма такой матрицы равна 100*0.5 = 50. Также положим t = 1. Из неравенства Коши-Буняковского следует, что степень из нормы можно вынести. Тогда неравенство в пункте 3 сведется к 50^i/i! < 0.00001. Осталось подобрать такое наибольшее значение числа i, чтобы оно выполнялось. Получается i_max = 149. Вы брали 200 членов — примерно соизмеримый результат. Если честно, то я не знаю — откуда берется такое время вычислений по сравнению с пакетами. Может, при замерах лишнее время где-то собрали (хотя, вряд ли). Оптимизация, которую здесь можно провести, — предварительное деление элементов матрицы на ее норму * t. Тогда неравенство из пункта 3 превратится в оценку 1/i! < eps. Но лучше делить на ближайшую степень двойки по избытку, чтобы потом результат последовательно возвести в квадрат. Ту же самую оптимизацию можно сделать и для рассматриваемого в топике метода.
                                                  0
                                                  Кстати, в последнем случае есть маленький недостаток. При последовательном возведении в степень результата матричной экспоненты собирается ошибка, при этом если значения элементов матрицы по модулю достаточно велики, то ошибка может быстро перейти к десятым, а то и хуже.
                                                    0
                                                    А собирается эта ошибка потому, что матричную экспоненту мы вычисляем с некоторой точностью.
                                                      0
                                                      Ошибка возникает из-за неточного представления исходных данных. Или, если мы их считаем точными двоично-рациональными числами, то из-за ошибки после первого же умножения.
                                                      Относительная ошибка произведения двух матриц ведёт себя так же, как и для чисел — она приблизительно разна сумме относительных ошибок сомножителей. Поэтому, относительная ошибка суммы степенного ряда (а она определяется ошибкой наибольшего члена, т.е. A^50), будет примерно 50*10^(-16).
                                                      Если мы возьмём матрицу B=A/64, и вычислим C=exp(B)-1=B+B^2/2+..., то её относительная ошибка будет примерно такой же, как у матрицы A (допустим, 10^(-16)). После вычисления (1+C)^64 (повторяя итерации C_{n+1}=C_n^2*2*C_n), мы получим, что погрешность увеличилась в 64 раза, и стала 64*10^(-16) — почти такой же, как у классического способа. Возможно, что-то мы ещё потеряли, но немного.
                                                      0
                                                      Точность вы какую брали — относительную, или абсолютную? Если абсолютную, как я понял из ваших рассуждений здесь, то вы что-то напутали — достигнуть её с double невозможно (или используется длинная арифметика?) Дело в том, что для матрицы 100х100 со всеми элементами равными 0.5 в матричной экспоненте все элементы будут порядка 1e20 — следовательно, абсолютная точность в 0.00001 никак не достижима с обычными числами с плавающей точкой.

                                                      Про число итераций — то, что у вас считается около 149 членов ряда (как я понял, хотя конкретно этого вы не написали) и это занимает столько времени, говорит об огромной неоптимальности реализации. Вы бы ещё на каком-нибудь питоне написали обычными циклами матричное умножение, и сравнивали скорость — выигрыш предлагаемого алгоритма был бы значительно больше :)
                                                        0
                                                        0.5 я взял для примера оценки количества членов. Соглашусь про реализацию — там допиливать надо.
                                                        +1
                                                        Точность вычисления матричной экспоненты была взята равной 0.00001. <...> Тогда норма такой матрицы равна 100*0.5 = 50.

                                                        В таком случае, норма экспоненты будет exp(50)=5*10^21. Того же порядка (примерно в 100 раз меньше) будет максимальный элемент матрицы. Если вычисления идут в типе double, то хватило бы точности 1000 (у нас только 15 знаков), и на это хватает 119 шагов.
                                                          0
                                                          Спасибо за интересное замечание, учту.
                                            0
                                            Это без загрузки результата с сервера WolframAlpha.
                                              0
                                              Хотя у m0nhawk матрица 500x500 считается полсекунды на домашней системе.
                                              • UFO just landed and posted this here
                                                  0
                                                  Согласен. Вообще, здесь нужно сравнивать не время работы алгоритма в конкретной системе, а по сути используемые методы.
                                                    0
                                                    Ну, раз заговорили о времени, мне прислали информацию о том, что для матрицы 140x140 экспонента в Octave считается доли секунды.
                                                  • UFO just landed and posted this here
                                      +2
                                      Прошу прощения, но так можно писать статьи с каждой лекции Уравнений Мат.Физики. Преподаватель наш, к примеру, максимально подробно объяснял задачу, методы решения, способ аппроксимации и сложность алгоритма. Если попросить — даже ссылки давал на пакеты, которыми можно продемонстрировать решения разными способами и сравнить.
                                      Поэтому спасибо за проделанную работу и статью, вспомнил универ и добро пожаловать в ряды математических факультетов всех, кто заинтересовался, там и нагляднее, и глубже и подробнее. Ну и лет на 6 дольше)
                                        0
                                        Кстати, для матрицы со случайными элементами от -1 до 1 экспонента с помощью ряда считается гораздо быстрее, чем для матрицы с элементами от 0 до 1. Судя по всему, собственные значения там гораздо меньше (порядка sqrt(N) вместо N/2). Для матрицы 500*500 с элементами из [0,1] потребовалось просуммировать 390 членов, а для матрицы с элементами из [-1,1] — только 55. Поэтому для времени вычислений в других пакетах интересно — какие матрицы им предлагали.
                                          0
                                          У вас отличие по числу членов меньше порядка, а на матрицах 100х100 оно будет ещё меньше. Время из графика в статье же как минимум на 2 порядка больше, чем должно быть.

                                          И я, к примеру, брал числа от 0 до 1, да и число итераций со значительным запасом.

                                        Only users with full accounts can post comments. Log in, please.