Linear System
LG
Не можете ли вы высказаться по вопросу, в каких случаях лучше
программировать решение линейной системы алгебраических уравнений
самому, а в каких использовать готовые библиотеки. И то и другое имеет
как свои плюсы так и минусы. И методов и библиотек в настоящее время
разработано большое количество. Как во всем этом разобраться?
Evgenii Rudnyi
Здравствуйте, Леонид.
Я должен сказать, что я химик по образованию. Однако использование
численных методов было значительной частью моей работы.
Поскольку я не математик, то вопрос, надо ли писать свою программу
решения линейных уравнений, даже не ставился. Ответ был для меня
очевиден - лучше всего взять что-то готовое, поскольку цель была другая
- решить конкретную задачу, где решение системы линейных уравнений было
только частью. А зачем изобретать велосипед?
К счастью, или к сожалению, это зависит от того, как посмотреть, особого
выбора численных библиотек в наше время и не было. Я начинал работать с
численными методами в 1979, во времена исторического материализма. У нас
был доступ к замечательному советскому суперкомьютеру БЭСМ-6 и
замечательный преподаватель, Михаил Григорьевич Анашкин. Скажем, как раз
где-то в 1981 году появился перевод книги Форсайта на русский
Computer Methods for Mathematical Computations, by Forsythe, Malcolm,
and Moler.
http://www.netlib.org/fmm/index.html
а у нас в группе эта книга появилась еще раньше прямо на английском языке.
Соответственно, подпрограммы decomp и solve и использовались для решения
системы линейных уравнений.
Время шло и на каком-то этапе я сменил подпрограммы из книги Форсайта на
LAPACK. Насколько я понимаю, многие, а может быть и большинство
коммерческих программ используют LAPACK - это говорит само за себя.
Далее на каком-то этапе я услышал слова оптимизированный BLAS. Честно
сказать, вначале не придал этому особого значения, поскольку по простоте
душевной казалось, что ну конечно можно превзойти код, оптимизированный
компилятором Фортрана, но не настолько же, чтобы этим серьезно
заниматься. В свое оправдание скажу, что мне требовалось решать
сравнительно небольшие системы линейных уравнений (скажем размерности до
20-30), и в этом случае оптимизированный BLAS все равно не играет особой
роли.
Следующий этап был связан с тем, что я более или менее регулярно
тестировал возможности новых компьютеров. На БЭСМ-6 в 1980 году
вычисление обратной матрицы размерности 20 занимало, если я правильно
помню, 17 секунд. Для меня это была отправная точка для измерения
развития компьютерных технологий во времени.
Так вот, по-моему в году 2000 я увидел, что с помощью ПК уже можно
решать и системы линейных уравнений размерностью 1000, но при этом мой
код, вызывающий LAPACK работает раз в пять медленнее, чем в MATLAB. Это
было особенно обидно, поскольку как раз где в это время MATLAB только
перешел на использование LAPACK, в то время как я использовал LAPACK и
раньше.
Пришлось вспомнить про волшебные слова оптимизированный BLAS, и тут мой
выбор был ATLAS. Это действительно помогло.
Я надеюсь, что эта небольшая история ответит на ваш вопрос и поможет
сделать вам правильный выбор.
С уважением,
Евгений Рудный
LG
Я предполагал, что вы достаточно молодой человек, которому просто
разобраться в обилии и сложностях работы с системами, библиотеками и
т.д. Но вы имеете опыт БЭСМ6, правда, у меня дополнительно к этому
БЭСМ4 и М220.
Вначале был опыт использования функций их Бате К. и Вилсона Е., затем
из Джорджа А. и Лю Дж., затем Ортега Д. На длительный период этого
было достаточно. Но использование современных программных пакетов
показало, что имеются большие наработки в решении линейных систем.
Сколько труда вложено в Metis, а ведь это только часть.
Поэтому стал вопрос об эффективном использовании готовых пакетов. Но
это сопряжено с определенными сложностями, обсуждение которых хотелось
бы увидеть. Мне кожется, что при сравнении боблиотек нужно было бы
включать пункт о сложность применения, хотя этот критерий
субъективный.
А как вы относитесь к итерационным методам?
KSergP
я считаю, что стоит писать свой софт, "заточенный" под свои задачи, а
не использовать универсальные алгоритмы, поскольку, существует
возможность его модернизации в зависимости от ваших идей. Да, есть и
открытый софт, но я так думаю найдется немного лудей готовых "рыться"
в нем.
хотя стоит использовать, например, BLAS как элементарный шаг более
"мощного" алгоритма.
необходимость перехода от прямых методов к итерационным для меня
очевидна. За последнее время в этой области сделано очень много
наработок. Выпускается очень много библиотек с различными реализациями
итерационных методов и предобусловливания, ориетированных на
разреженные СЛАУ. Но, от себя, отмечу, что мне не известен софт
позволяющий без дополнительной модификации использовать эффективные
итерационные методы (крыловского типа) решения СЛАУ с плотной
матрицей. Так, например, в матлабе говорится....для решения плотных
СЛАУ используется LAPACK.
Всего доброго.
Evgenii Rudnyi
> я считаю, что стоит писать свой софт, "заточенный" под свои задачи, а
> не использовать универсальные алгоритмы, поскольку, существует
> возможность его модернизации в зависимости от ваших идей.
Мне кажется стоит разделить это на отдельных утверждения.
1) Универсальный алгоритм или алгоритм для конкретного случая.
Наверное, универсального решения здесь нет. Однако здесь не надо
забывать про время разработки. Можно улучшить производительность, если
найти нечто особенное в своей задаче, и это потом использовать. Однако,
это требует время и вопрос здесь в том, будет ли вложенное время
оправдано. Скажем, если эту задачу требуется решать не так уж и часто,
то может быть проще применить универсальный алгоритм.
Я недавно случал интересный доклад на семинаре NAFEMS
http://www.nafems.de/Veranstaltungen/CFD07/cfd07_agenda.html
Application Specific CFD Codes versus Application Specific User
Interfaces – Where´s the Difference?
R. D. Löffler (Fluent Deutschland GmbH)
Здесь как раз обсуждался похожий вопрос. Так вот, один из выводов был,
что один универсальный код обходится дешевле, чем много кодов на каждый
специальный случай.
2) Свой код или чужой.
Моя позиция здесь сводится к тому, что изобретать велосипед не самое
лучшее занятие. Если уже что-то есть сделанное и это можно использовать,
то мне представляется это крайне разумнум решением.
>Да, есть и
> открытый софт, но я так думаю найдется немного лудей готовых "рыться"
> в нем.
> хотя стоит использовать, например, BLAS как элементарный шаг более
> "мощного" алгоритма.
Я всю жизнь использую библиотеки. Скажем, мой последний код
http://www.imtek.uni-freiburg.de/simulation/mor4ansys/
без использования библиотек было бы невозможно написать за то время,
которое было в наличии.
> необходимость перехода от прямых методов к итерационным для меня
> очевидна.
В моем приложении мне требовались именно прямые методы, поскольку надо
было много раз решать систему уравнений с одной матрицей но с разными
RHS (right hand side).
Также прямые методы гораздо более устойчивы. Если хотите, их поэтому
иногда называют black box preconditioner. Также современные
многофронтальные методы с сочетании с тем, что память дешевая, приводят
к тому, что по времени они не так уж и сильно проигрывают, в особенности
для реальных трехмерных задач.
> Так, например, в матлабе говорится....для решения плотных
> СЛАУ используется LAPACK.
А что, скажем если взять матрицы до размерности несколько тысяч, вы
считаете, что есть итерационные методы, которые побьют LAPACK?
С уважением,
Евгений
KSergP
LG
Я знаком с работами Сергея. Есть такие случаи.
Пусть плотная матрица имеет диагональные элементы, которые значительно
превосходят
по модулю остальных. Простая итерация будет самым эффективным методом.
Ну, а более серьезно, если формально матрица плотная, а по сути, если
пренебречь
малыми элементами-разреженная. Для нее может быть предпочтительный
итерационный метод.
Интереснее проблема использования методов для задачи о собственных
значениях.
Большие задачи, скорее всего нужно решать итерационно. Но нужен
хороший предобуславливатель,
который минимизирует число итераций. Вопрос, где его взять?
Всего доброго. Леонид.
KSergP
>
>> В моем приложении мне требовались именно прямые методы, поскольку надо
>> было много раз решать систему уравнений с одной матрицей но с разными
>> RHS (right hand side).
>
>
>> Также прямые методы гораздо более устойчивы.
>
примеры посмотреть...
>
>> Если хотите, их поэтому
>> иногда называют black box preconditioner.
>
решением СЛАУ прямыми методами и предобусловливанием.... если можно по
подробней...
>
>
>> А что, скажем если взять матрицы до размерности несколько тысяч, вы
>> считаете, что есть итерационные методы, которые побьют LAPACK?
>
другой вопрос...в том, что сравнением скорее всего может оказаться не
совсем корректным....
поскольку для реализации надо использовать одни и теже алгоритмы
"ускорения", что и используются в LAPACKе...
С уважением, Сергей.
Evgenii Rudnyi
>>> В моем приложении мне требовались именно прямые методы, поскольку надо
>>> было много раз решать систему уравнений с одной матрицей но с разными
>>> RHS (right hand side).
> хм... но в принципе существуют блок версии итерационных методов.
Буду признателен, если сообщите ссылки про эти методы. Хотелось бы
понять, на чем они основаны.
А как там растет время вычислений в зависимости от числа векторов в
правой части?
>>> Также прямые методы гораздо более устойчивы.
> интересное замечание, что прямые методы более устойчивы....хотелось бы
> примеры посмотреть...
Я имел в виду слово robustness. При использовании прямых методов я
гарантированно получаю решение (или сообщение о том, что матрица
вырождена) за предсказуемое число операций и это не зависит от матрицы.
Вы согласны с этим утверждением?
В случае итерационных методов все зависит от preconditioner, то есть для
произвольной матрицы решение за заданное число операций не
гарантируется. Также проблема здесь определить, что матрица вырождена.
>>> Если хотите, их поэтому
>>> иногда называют black box preconditioner.
> что-то я не совсем понял, упомянутую вами связь между устойчивым
> решением СЛАУ прямыми методами и предобусловливанием.... если можно по
> подробней...
Соответственно, универсальный preconditioner для всех матриц обычно
называют black box preconditioner и это мечта всех разработчиков
итерационных методов. Когда они обсуждают этот вопрос, то они говорят,
что один black box preconditioner - это прямой метод, поскольку здесь
для всех матриц гарантированно получается решение. Это конечно
утрировано, поскольку в этом случае итерационной метод и не требуется,
но если посмотреть формально, то лучший preconditioner - это обратная
матрица.
Это я не сам выдумал, это из докладов разработчиков итерационных
методов, которые я слышал.
>>> А что, скажем если взять матрицы до размерности несколько тысяч, вы
>>> считаете, что есть итерационные методы, которые побьют LAPACK?
> если чесно, то абсолютно уверен...
> другой вопрос...в том, что сравнением скорее всего может оказаться не
> совсем корректным....
> поскольку для реализации надо использовать одни и теже алгоритмы
> "ускорения", что и используются в LAPACKе...
Мне кажется, это очень просто попробовать. Но я боюсь, что для случайной
матрицы это не выйдет, поскольку, насколько я понимаю, black box
preconditioner, за исключением прямых методов, не существует.
С уважением,
Евгений Рудный
P.S. Мне кажется, что здесь затронуто несколько разных тем, и если вы
захотите продолжить обсуждение, я бы предложил сделать новые темы,
скажем, блочные итерационные методы, black box preconditioner, побить
LAPACK.
KSergP
Здравстуйте, Евгений.
> Буду признателен, если сообщите ссылки про эти методы. Хотелось бы
> понять, на чем они основаны.
например, вот
http://www-sbras.nsc.ru/EMIS/journals/ETNA/vol.16.2003/pp129-142.dir/pp129-142.pdf
> А как там растет время вычислений в зависимости от числа векторов в
> правой части?
>
я, "по молодости ))", уже как в 2005г. экспериментировал.....
но, как то все это пока отложено.... и эти результаты не вошли в мою
диссертацию....
в ближайшее время вышлю Вам мою статью, посвященную этой задаче...
> >>> Также прямые методы гораздо более устойчивы.
> > интересное замечание, что прямые методы более устойчивы....хотелось бы
> > примеры посмотреть...
>
> Я имел в виду слово robustness. При использовании прямых методов я
> гарантированно получаю решение (или сообщение о том, что матрица
> вырождена) за предсказуемое число операций и это не зависит от матрицы.
> Вы согласны с этим утверждением?
Точные методы более предсказуемы, я бы так сказал....
> >>> Если хотите, их поэтому
> >>> иногда называют black box preconditioner.
> > что-то я не совсем понял, упомянутую вами связь между устойчивым
> > решением СЛАУ прямыми методами и предобусловливанием.... если можно по
> > подробней...
>
> Соответственно, универсальный preconditioner для всех матриц обычно
> называют black box preconditioner и это мечта всех разработчиков
> итерационных методов. Когда они обсуждают этот вопрос, то они говорят,
> что один black box preconditioner - это прямой метод, поскольку здесь
> для всех матриц гарантированно получается решение. Это конечно
> утрировано, поскольку в этом случае итерационной метод и не требуется,
> но если посмотреть формально, то лучший preconditioner - это обратная
> матрица.
да обратная матрица, полученная как можно с меньшими затратами....
> Это я не сам выдумал, это из докладов разработчиков итерационных
> методов, которые я слышал.
вообщем, с такой формулировкой я согласен))
хотелось бы узнать, где это Вы все прослушивали?
> P.S. Мне кажется, что здесь затронуто несколько разных тем, и если вы
> захотите продолжить обсуждение, я бы предложил сделать новые темы,
> скажем, блочные итерационные методы, black box preconditioner, побить
> LAPACK.
думаю пока не стоит....если к нам еще бы присоеденилось несколько
человек, тогда модет и стоило бы...
я тут дописываю свою прогу.... некоторые мелочи....
и хочу сравнить свой не оптимизированный код с MatLab-ом в котором
подключен LAPACK.... только не знаю успею до защиты или нет...
Evgenii Rudnyi
>> Это я не сам выдумал, это из докладов разработчиков итерационных
>> методов, которые я слышал.
>
>
> вообщем, с такой формулировкой я согласен))
> хотелось бы узнать, где это Вы все прослушивали?
Мой интерес к итерационным методам был достаточно опосредованным. Я
работал в университете Фрайбурга по теме связанной с model reduction
и здесь тоже используются подпространства Крылова, но теория совсем
другая. Из конференций, где я был, очень интересная была
Latsis Symposium 2002 on Iterative Solvers for Large Linear Systems
("CG50--GG70")
Это событие было посвящено 50 метода conjugated gradients и там были
представлены очень многие из работающих над итерационными методами. К
сожалению в Интернете почти ничего не осталось. Я нашел только объявление
в NA Digest
http://www.netlib.org/na-digest-html/01/v01n42.html#10
и несколько ссылок на страничке профессора Gutknecht
http://www.sam.math.ethz.ch/~mhg/
где-то внизу странички (проще запустить поиск на Latsis Symposium 2002).
Далее интересный доклад про black box preconditioner был на
MACSI-net-Workshop 2004
http://www-num.math.uni-wuppertal.de/conferences/macsinet.html
Frommer, A.: Towards black-box preconditioners for linear systems
http://www-num.math.uni-wuppertal.de/conferences/frommer_talk.pdf
Из близких тем был семинар Workshop on Exponential Integrators 2004
http://techmath.uibk.ac.at/numbau/alex/events/conference2004.html
Из больших конференций мне понравилась PARA 2004
http://www2.imm.dtu.dk/~jw/para04/
Эта конференция проводится каждые два года, и если у вас будет
возможность, я очень ее рекомендую.
С уважением,
Евгений
KSergP
Да, Евгений, весомый список
ничем подобным похвастаться не могу)
чесно говоря, я уже устал от конференций на которых участвую,
поскольку создается впечатление, что никому это не интересно и не
понятно зачем вообще что-то новое придумывать....
тут набросал одну простенькую статью (думаю потом значительно усилить
и опубликовать в "значительном" издании)....
в ней приведены плотные матрицы после предфильтрации....
получились достаточно интересные портреты...
если кому интересно, могу выслать....
да еще инетесная вещь (по кранйней мере для меня).... нашел статью в
которой для решения СЛАУ использовали генетические алгоритмы.....
KSergP
я доделал провел небольшое сравнение
своих результатов с результатами полученными LAPACK-ом, встроенным в
MatLab....
я буду использовать названия, которые я использовал в автореферате и
диссертации.
использованный прямой метод: метод исключения Гаусса без
упорядочивания (GE)
использованный итерационный метод: BiCGStab
способ формирования матрицы предобусловливания: полное LU
способ предфильтрации: (2.1), основанный на нормах строк матрицы
1 пример.
широкодиапазонная антенна, N=2909
MatLab-50c.
TALGAT(GE)-464c.
TALGAT(BiCGStab)-65c.
2.Трапецивидно зубчатая антенна, N=2975
MatLab-54c.
TALGAT(GE)-483c.
TALGAT(BiCGStab)-50c.
3. диполь с углом между лучами 180 градусов, , N=3001
MatLab-56c.
TALGAT(GE)-501c.
TALGAT(BiCGStab)-27c.
Я для реализации методов не использовал никакие "ускорялки".....т.е. я
"не привязывался" к "железу"...
мои вывод: если осуществить привязку к железу итерационные методы
можно значительно ускорить.....
так разница между "моим" не ускоренным Гауссом и "суперускоренным"
LAPACK-овским Гауссом достигает почти 10 раз....
Evgenii Rudnyi
Несколько соображений.
1) Это показывает, что имеет смысл начинать с использования библиотек.
Если что-то уже сделано и в это вложено человеко-года, то, наверное, это
будет работать быстрее. Это, конечно, не означает, что нельзя сделать
лучше, но это требует времени, которого обычно и нет.
2) Выражение ускоренный Гаусс не совсем правильно, поскольку там на самом
деле используются другие алгоритмы, основанные на блочном подходе. Можно
сказать блочный Гаусс, но это уже все равно другой алгоритм. Это уже потом
блочные алгоритмы привязываются к железу. В любом случае основное
ускорение появляется в замене исходного алгоритма.
3) Основное ускорение появляется при использовании уровня BLAS 3, а
итерационные методы могут использовать похоже только BLAS 2. Отсюда и
сомнения, что удастся достигнуть той же эффективности. Здесь вполне может
быть, что вам потребуется меньше операций, но вы не сможете их выполнить
на максимальной мощности процессора, как это будет в LAPACK.
4) Здесь также можно подумать и про другие возможности. Если у вас плотные
матрицы, то значит вы наверное использует метод граничных элементов
(boundary elements). В этом случае можно использовать multipole expansion,
чтобы существенно ускорить вычисление произведения матрица вектор. У Jacob
White
http://www.rle.mit.edu/rleonline/People/JacobK.White.html
есть серий статей и программ в этом направлении
http://www.rle.mit.edu/cpg/research_codes.htm
С другой стороны, а так ли хорош BEM? Может быть решать в конечных
элементах и лучше.
5) Это приводит нас к учению Владимира Ильича о слабом звене. Было бы
неплохо понять, а что вы хотите достигнуть, и что на этом пути является
слабым звеном. Скажем, если есть желание сделать программу, а потом ее
продавать, то обычно слабым звеном является маркетинг. Без него никуда.
Очень частое явление в западных компаниях, это когда вначале находится, а
что можно продать, а только потом это уже делается.
Успехов,
Евгений
KSergP
> 4) Здесь также можно подумать и про другие возможности. Если у вас плотные
> матрицы, то значит вы наверное использует метод граничных элементов
> (boundary elements). В этом случае можно использовать multipole expansion,
> чтобы существенно ускорить вычисление произведения матрица вектор. У Jacob
> White
>
> http://www.rle.mit.edu/rleonline/People/JacobK.White.html
>
> есть серий статей и программ в этом направлении
>
> http://www.rle.mit.edu/cpg/research_codes.htm
>
> С другой стороны, а так ли хорош BEM? Может быть решать в конечных
> элементах и лучше.
нет...у меня метод моментов.... хотя и к нему применима идея
увеличения матрично-векторного перемножения.....
правда это возможно если матрица имеет определенную структуру.....
например тоэплицевскую...
> 5) Это приводит нас к учению Владимира Ильича о слабом звене. Было бы
> неплохо понять, а что вы хотите достигнуть, и что на этом пути является
> слабым звеном. Скажем, если есть желание сделать программу, а потом ее
> продавать, то обычно слабым звеном является маркетинг. Без него никуда.
> Очень частое явление в западных компаниях, это когда вначале находится, а
> что можно продать, а только потом это уже делается.
лично для меня... интересная идея повышения скорости
быстродействия.....и лично я не "болею" манией продажы своего
софта..... хотя реальность иногда заставляет об этом задуматься.....