Комментарии 47
В случае диагональных матриц матричная экспонента состоит из скалярных экспоненциальных функций. В формуле (8) коэффициенты r будут неявно определять часть их разложения в степенной ряд. Когда значения показателей степени в таких функциях достаточно велики, то наш алгоритм проиграет — слишком много членов разложения в ряд придется брать, здесь лучше оптимизировать вычисление скалярных экспонент, выделив целую и дробную части в показателе.
Порядок матрицы изменялся от 2 до 100.
А на графике порядок по оси абсцисс продолжается до 140, а точки следуют, судя по всему, до 130-го порядка. Это опечатка или результат экстраполяции?
В подписи к рисунку написано, что классический алгоритм — это верхняя кривая, то есть с крестиками, а на самом рисунке — что это кривая со штрихами, то есть нижняя. В тексте статьи об этом не написано ничего. Так что же показало сравнение быстродействия алгоритмов?
Одно из основных следствий из закона Мэрфи гласит, что, если что-то может быть не понято, то оно будет непонято.
Когда я ещё был молодым учёным и специалистом, умные дяди меня учили, что именно поэтому крестики — очень плохой выбор. А главное — обсуждение каждого графика и вывод по каждому графику должны быть в тексте.
Когда я ещё был молодым учёным и специалистом, умные дяди меня учили, что именно поэтому крестики — очень плохой выбор. А главное — обсуждение каждого графика и вывод по каждому графику должны быть в тексте.
Честно говоря, неубедительный прирост скорости, с учётом возросших на порядок требований к памяти.
Хм, выходит никто не считает e^At с помощью жордановой формы?
Мне на почту пришло сообщение о том, что в пакете MATLAB вычисление матричной экспоненты для матрицы 140-ого порядка, заполненной случайными числами в интервале от (0;1), занимает 0.005371 секунд (пакет выдал такое время). Используемые команды
В описании команды expm указано, что пакет использует модификацию аппроксимации Паде, рассмотренную в работе Nicholas J. Higham «The Scaling and Squaring Method for the Matrix Exponential Revisited». Наверное, такая скорость вполне может возможна. Я попробовал вычислить матричную экспоненту в WolframAlpha. Например, для матрицы 60x60
вычисления идут несколько секунд (прямая ссылка). Для матрицы 140-ого порядка система не выдала результата.
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-ого порядка система не выдала результата.
Простенькая формула Exp(A)=Exp(A/(2^N))^(2^N) позволила ускорить поиск экспоненты матрицы 500*500 в 3.4 раза (с 16.2 до 4.8 сек на C# безо всяких оптимизаций). Относительная погрешность результата — лучше 10^-13.
Совсем не впечатлило ускорение, вот например в scipy функция expm3 (считает экспоненту через разложение в ряд) вычисляет для матрицы 150х150 всего за 13 мс (и это на достаточно старом ноуте). Хотя я сам реализацию её не смотрел, но в документации написано именно про наивный алгоритм.
А правильный результат он дает на таких размерах?
Каюсь, не проверил — действительно, по-умолчанию даёт весьма неточный результат. Однако, можно выбирать любое число членов ряда (странно, что автоматически результат плохой — видимо, из-за того, что это не основной метод в библиотеке для этого), и в данном случае 200 членов с большим запасом хватает — относительная точность порядка 1e-15. Считается же теперь дольше, но намного быстрее вашего варианта — 150 мс.
Для сравнения, обычный expm оттуда же считает около 15 мс и точно — он как раз использует ту же аппроксимацию, что и вышеупомянутый матлаб.
Для сравнения, обычный expm оттуда же считает около 15 мс и точно — он как раз использует ту же аппроксимацию, что и вышеупомянутый матлаб.
+ наверное, формула, указанная здесь habrahabr.ru/post/233587/#comment_8036405 Точность берется по норме общего члена ряда?
Судя по исходникам scipy, там просто наивный алгоритм суммирования степеней, и всё (да и вроде по времени сходится — матричное умножение A.dot(A) для такой матрицы у меня занимает около 0.6 мс, что при повторении 200 раз будет 120 мс). Про точность вопрос не понял, я сравнивал с expm и expm2 — они используют различные методы, но все три функции выдают практически совпадающие результаты.
Про точность см. неравенство в пункте 3 этого топика. Похоже, что при суммировании там берется фиксированное количество членов ряда.
А, ну да, я же говорю — указывается число членов вручную, автоматический выбор места остановки там не реализован. Просто странно, что у вас наивный метод работает на два порядка дольше.
По-моему, ничего странного, если алгоритм не детерминирован. Я не зря сделал сравнение с временем у системы WolframAlpha.
Всмысле, не детерминирован? Вроде как все шаги вполне однозначны, значит алгоритм детерминирован. Да и как ни крути, в данном случае отличия на два порядка ну никак не должно быть. Сколько у вас членов суммируется в итоге?
С вольфрамальфа как раз зря сравнение сделали — кто знает, с какими приоритетами и на каких мощностях там запускаются такие вычисления, тем более бесплатно… Он же совсем не для этого заточен. Нужно сравнивать с локальными программами/библиотеками.
С вольфрамальфа как раз зря сравнение сделали — кто знает, с какими приоритетами и на каких мощностях там запускаются такие вычисления, тем более бесплатно… Он же совсем не для этого заточен. Нужно сравнивать с локальными программами/библиотеками.
Заканчивает свое исполнение при выполнении неравенства, при этом мы заранее не знаем, сколько итераций будет сделано.
Не понял, к чему именно этот ответ? Если к вопросу
Сколько у вас членов суммируется в итоге?, то несмотря на неизвестность числа членов заранее, постфактум это всё же можно подсчитать.Да, оценив заранее норму матрицы.
Ну так, сколько же у вас шагов выполняется?
Просто как-то складывается впечатление, что взяли заведомо неоптимальную реализацию наивного алгоритма, и показываете, что ваш работает быстрее. Кстати, сравнения точности в статье не увидел.
Просто как-то складывается впечатление, что взяли заведомо неоптимальную реализацию наивного алгоритма, и показываете, что ваш работает быстрее. Кстати, сравнения точности в статье не увидел.
Точность вычисления матричной экспоненты была взята равной 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. Но лучше делить на ближайшую степень двойки по избытку, чтобы потом результат последовательно возвести в квадрат. Ту же самую оптимизацию можно сделать и для рассматриваемого в топике метода.
Кстати, в последнем случае есть маленький недостаток. При последовательном возведении в степень результата матричной экспоненты собирается ошибка, при этом если значения элементов матрицы по модулю достаточно велики, то ошибка может быстро перейти к десятым, а то и хуже.
А собирается эта ошибка потому, что матричную экспоненту мы вычисляем с некоторой точностью.
Ошибка возникает из-за неточного представления исходных данных. Или, если мы их считаем точными двоично-рациональными числами, то из-за ошибки после первого же умножения.
Относительная ошибка произведения двух матриц ведёт себя так же, как и для чисел — она приблизительно разна сумме относительных ошибок сомножителей. Поэтому, относительная ошибка суммы степенного ряда (а она определяется ошибкой наибольшего члена, т.е. 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) — почти такой же, как у классического способа. Возможно, что-то мы ещё потеряли, но немного.
Относительная ошибка произведения двух матриц ведёт себя так же, как и для чисел — она приблизительно разна сумме относительных ошибок сомножителей. Поэтому, относительная ошибка суммы степенного ряда (а она определяется ошибкой наибольшего члена, т.е. 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) — почти такой же, как у классического способа. Возможно, что-то мы ещё потеряли, но немного.
Точность вы какую брали — относительную, или абсолютную? Если абсолютную, как я понял из ваших рассуждений здесь, то вы что-то напутали — достигнуть её с double невозможно (или используется длинная арифметика?) Дело в том, что для матрицы 100х100 со всеми элементами равными 0.5 в матричной экспоненте все элементы будут порядка 1e20 — следовательно, абсолютная точность в 0.00001 никак не достижима с обычными числами с плавающей точкой.
Про число итераций — то, что у вас считается около 149 членов ряда (как я понял, хотя конкретно этого вы не написали) и это занимает столько времени, говорит об огромной неоптимальности реализации. Вы бы ещё на каком-нибудь питоне написали обычными циклами матричное умножение, и сравнивали скорость — выигрыш предлагаемого алгоритма был бы значительно больше :)
Про число итераций — то, что у вас считается около 149 членов ряда (как я понял, хотя конкретно этого вы не написали) и это занимает столько времени, говорит об огромной неоптимальности реализации. Вы бы ещё на каком-нибудь питоне написали обычными циклами матричное умножение, и сравнивали скорость — выигрыш предлагаемого алгоритма был бы значительно больше :)
Точность вычисления матричной экспоненты была взята равной 0.00001. <...> Тогда норма такой матрицы равна 100*0.5 = 50.
В таком случае, норма экспоненты будет exp(50)=5*10^21. Того же порядка (примерно в 100 раз меньше) будет максимальный элемент матрицы. Если вычисления идут в типе double, то хватило бы точности 1000 (у нас только 15 знаков), и на это хватает 119 шагов.
Это без загрузки результата с сервера WolframAlpha.
Хотя у m0nhawk матрица 500x500 считается полсекунды на домашней системе.
Прошу прощения, но так можно писать статьи с каждой лекции Уравнений Мат.Физики. Преподаватель наш, к примеру, максимально подробно объяснял задачу, методы решения, способ аппроксимации и сложность алгоритма. Если попросить — даже ссылки давал на пакеты, которыми можно продемонстрировать решения разными способами и сравнить.
Поэтому спасибо за проделанную работу и статью, вспомнил универ и добро пожаловать в ряды математических факультетов всех, кто заинтересовался, там и нагляднее, и глубже и подробнее. Ну и лет на 6 дольше)
Поэтому спасибо за проделанную работу и статью, вспомнил универ и добро пожаловать в ряды математических факультетов всех, кто заинтересовался, там и нагляднее, и глубже и подробнее. Ну и лет на 6 дольше)
Кстати, для матрицы со случайными элементами от -1 до 1 экспонента с помощью ряда считается гораздо быстрее, чем для матрицы с элементами от 0 до 1. Судя по всему, собственные значения там гораздо меньше (порядка sqrt(N) вместо N/2). Для матрицы 500*500 с элементами из [0,1] потребовалось просуммировать 390 членов, а для матрицы с элементами из [-1,1] — только 55. Поэтому для времени вычислений в других пакетах интересно — какие матрицы им предлагали.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
О вычислении матричной экспоненты