Comments 47
UFO just landed and posted this here
UFO just landed and posted this here
В случае диагональных матриц матричная экспонента состоит из скалярных экспоненциальных функций. В формуле (8) коэффициенты r будут неявно определять часть их разложения в степенной ряд. Когда значения показателей степени в таких функциях достаточно велики, то наш алгоритм проиграет — слишком много членов разложения в ряд придется брать, здесь лучше оптимизировать вычисление скалярных экспонент, выделив целую и дробную части в показателе.
Порядок матрицы изменялся от 2 до 100.
А на графике порядок по оси абсцисс продолжается до 140, а точки следуют, судя по всему, до 130-го порядка. Это опечатка или результат экстраполяции?
В подписи к рисунку написано, что классический алгоритм — это верхняя кривая, то есть с крестиками, а на самом рисунке — что это кривая со штрихами, то есть нижняя. В тексте статьи об этом не написано ничего. Так что же показало сравнение быстродействия алгоритмов?
UFO just landed and posted this here
Одно из основных следствий из закона Мэрфи гласит, что, если что-то может быть не понято, то оно будет непонято.
Когда я ещё был молодым учёным и специалистом, умные дяди меня учили, что именно поэтому крестики — очень плохой выбор. А главное — обсуждение каждого графика и вывод по каждому графику должны быть в тексте.
Когда я ещё был молодым учёным и специалистом, умные дяди меня учили, что именно поэтому крестики — очень плохой выбор. А главное — обсуждение каждого графика и вывод по каждому графику должны быть в тексте.
Честно говоря, неубедительный прирост скорости, с учётом возросших на порядок требований к памяти.
Хм, выходит никто не считает 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. Поэтому для времени вычислений в других пакетах интересно — какие матрицы им предлагали.
Sign up to leave a comment.
О вычислении матричной экспоненты