ILUT предобусловливание

[KSergP]

Кто-нибудь использовал ILUT?
В Saad=ом описании не понятные строки.

[Nizky]

Вам что-нибудь удалось с ILUT?
Мне кажется очень странным, что при первом dropping элемент сравнивается с 2-нормой строки исходной матрицы. Ведь мы по сути применяем этот критерий к найденному элементу матрицы L, которая не меняется при домножении исходной матрицы на число. В то время как 2-норма строки меняется.

[KSergP]

Nizky писал(а):

Вам что-нибудь удалось с ILUT?
Мне кажется очень странным, что при первом dropping элемент сравнивается с 2-нормой строки исходной матрицы. Ведь мы по сути применяем этот критерий к найденному элементу матрицы L, которая не меняется при домножении исходной матрицы на число. В то время как 2-норма строки меняется.

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

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

[Nizky]

Мне тоже казалось, что логичней отбрасывать их по абсолютной величине. Все-таки, при вычислении элементов L происходит деление на диагональный элемент, в некотором роде нормирование.
Логики в применении к элементам L критериев, учитывающих какие-либо нормы исходной матрицы, не вижу.
Второе отсеивание по модулю тоже вроде не очень нужно: первое отсеивание гарантирует экономию времени, а то, что мы выбираем в каждой строке только наибольшие, экономит нам память.
Исходная система у меня плотная. Но умножение на матрицу производится приближенно и неявно. Явно известна только небольшая часть матрицы. Вот ее и хочется "неполно разложить".
Есть уверенность, что если эту часть матрицы честно обратить, то выйдет хороший предобуславливатель. Так что в этом смысле можно считать, что матрица у меня редкая. Только запас прочности меньше.
Можно продолжить обсуждение по личной почте: mailto:[e-mail].
А можно и здесь - авось кто еще заинтересуется.
С уважением, Низкий Роман.
г. Санкт-Петербург.

[KSergP]

Чесно говоря я тоже не понимаю зачем применять к уже вычисленным значениям пороги основанные на норме исходной матрицы.
В общем я понял что у вас почти такая же проблема как у меня, плотные матрицы и как к ним применнить предобуславливание.
У них же ILUT применяется к разреженным матрицам.
У вас в каом приложении возникают плотные матрицы?
И что значит перемножение матрицы на вектор выполняется не явно?

[Nizky]

Матрица у меня плотная. Чтобы применить неполное разложение и прочие методы редких матриц, я строю предобуславливатель на основе части исходной матрицы, т.е. матрицы, полученной из исходной обнулением большого числа элементов. Полученная матрица, конечно, не настолько приближает исходную, чтобы вместо исходной решать редкую систему. Зато есть надежда, что для построения предобуславливателя она пойдет. Основное отличие ситуации от того, с чем работал Saad, в следующем. У них в идеале, в пределе (при точном обращении, при точном LU-разложении типа ILU(inf) или ILUT(N,inf)) получается предобуславливатель, с которым все сходится за 1 итерацию. Загрубляют разложнеие - число итераций несколько увеличивается. В моем же случае даже идеальное разложение дает существенное число итераций, а при загрублении разложения все может вылиться вообще в отсутствие сходимости.
Основные претензии к ILUT - для матриц L и U применяются принципиально различные критерии отсеивания. А точнее - к матрице L они строже, чем к U. В итоге эти две матрицы заполнены совсем по-разному. У меня эта штука работает, но p при этом настолько велико, что никакого выигрыша по сравнению с точным LU не ощущается.
Сейчас я пробую ILU(p). Во-первых, по одномерному, да еще дискретному семейству параметров бегать гораздо легче, чем по двумерному. Во-вторых, пока я решаю задачи с симметричной матрицей, так что при точном LU получается U=L^T. ILU(p) хотя бы дает матрицы L и U c одинаковым числом ненулевых элементов.
Ключевые слова в моей области - "метод граничных элементов".
Исследуется рассеяние произвольного излучения на идеально проводящих (т.е. идеальных металлических) поверхностях.
Записывается "интегральное уравнение электрического поля", которое связывает индуцированный (получившийся вследствие падающего излучения) электрический ток по поверхности с падающим излучением.
Потом применяется "метод моментов" с базисными функциями предложенными Rao и компанией. В результате получается плотная комплескная система, где переменные соответствуют току, перетекающему через конкретное ребро одного из треугольников, на которые разбита поверхность. В интересующем случае размерность 10000 - 100000 переменных.
Неявное произведение выполняется с помощью так называемого fast multipole method (Vladimr Rokhlin, конец 80х). Используется происхождение матрицы. Грубо говоря, в моем случае сосчитать произведение матрицы на вектор это то же самое, что определить потенциал некоторого распределения зарядов в точках, в которых располагаются заряды. Вектор, на который умножают, соответствует величинам зарядов. Поле каждого заряда можно разложить в ряд, а при складывании полей нескольких зарядов просто складывать коэффициенты ряда. Все достаточно сложно, но позволяет выполнять умножение за O(n log n) операций. Побочный эффект - некоторые элементы матриц нельзя учесть, так как соответствующие пары точек (т.е. ребер треугольников) находятся слишком близко. Эта часть выделяется в отдельное слагаемое - "ближнюю матрицу", и при умножении на вектор это слагаемое умножается честно. Собственно, на основе этой "ближней матрицы" я и строю предобуславливатель.
По этой задаче очень мало написано по-русски, из тех людей, чьи фамилии встречались: Смирнов Ю.Г., Ивахненко В.И., Ильинский А.С..
Сам я этим занимаюсь в ЦНИИ им. А.Н.Крылова (уже года два), но институт занимается судостроением, так что слово "preconditioning" в нем кроме моих начальников (Вишневский А.М., Лаповок А.Я.), почти ничего никому не скажет.
Интересно, где с предобуславливателями столкнулись Вы. И вообще, как найти в России специалистов по данному вопросу.

[KSergP]

Т.е если я правильно понял, вы используете «ближнею» матрицу в качестве предобусловливателя?
Таким образом, вы предобуславливатель получаете, используя структуру решаемой задачи?
Подобный подход я находил в одной статье (на англ. языке) если надо могу выслать. Они там использовали ILUT.
Каким методом вы решаете? (я так понимаю методом сопряженных градиентов?).
И как вы поступили с арифметикой? (просто действительную арифметику заменили на комплексную?)
Скажите можно ли использовать подход Рохлина для несимметричных матриц?
Насколько я понял (из статей) они используют этот подход только при симметричных матрицах.

О своей проблеме. Решение задач излучения или рассеяния электромагнитной волны сложными объектами может быть получено с помощью интегральных уравнений, сводящихся методом моментов к СЛАУ с плотными матрицами J) (ну как?). В качестве примеров использую задачу определения токов в проводной антенне (различной конфигурации). Сами понимаете приходится решать СЛАУ с несимметричной, комплексной матрицей. Поэтому и приходится искать спасение в итерац. методах с предобусловливанием. Похоже проблемы у нас похожи. Проведенный мной обзор показывает, что итерационные методы хорошо изучены для разреженных матриц, но мало работ посвященных решению СЛАУ с плотными матрицами. В отечественной литературе кроме упоминания об BiCG я ничего более не встречал. Ну еще какие-то способы предобусловливания (основанные на обращении матрицы М) можно найти в записках ПОМИ (С.-Петербуржское отделение математического института им. Стеклова). Кстати упоминание о Колотилиной иногда можно найти в «буржуйской» литературе. Но контакт я пока с ними не наладил.

[Nizky]

Все-таки хорошо, что этот форум индексируется в Яндексе - если кто в России заинтересуется, то точно нас найдет.
Да, я использую в качестве предобуславливателя "ближнюю" матрицу, так как другого выхода у меня нет. Как я понимаю, все в моем случае так делают. Вот только согласия нет о том, какой предобуславливатель лучше.
Первый метод, которым я воспользовался, был BiCGSTAB (Van der Vorst). У меня он не заработал: на малых частотах сходился как-то "не до конца", а на больших вообще ерунду показывал. Ни один из семейства CG-методов так нормально и не заработал. Потом были различные вариации QMR (Freund). Начало сходиться. Скопировал свои матрицы в MatLab - с помощью его методов понял, что так и должно быть, но GMRES еще лучше. Реализовал GMRES (Saad), им пока и пользуюсь.
В моей задаче размерность прямо связана с частотой излучения - размер элемента должен быть пропорционален длине волны. Где-то около 5000 переменных сходимость стала очень грустной - 40 dB за 200 итераций, и стало ясно, что потом будет еще хуже. Так я и начал искать предобуславливатели.
В BiCG долго стоял вопрос, какое скалярное умножение использовать: обычное v на v^T или эрмитово v на v^*? Вопрос не был решен, так как все равно не сошлось, но я все же склоняюсь в пользу эрмитова. В описании QMR говорилось о комплексных числах в явном виде. В GMRES надо использовать эрмитово произведение и комплексное вращение Гивенса.
Думаю, что симметрия матрицы для FMM (fast multipole method) не принципиальна. Сам Рохлин это придумывал для потенциалов, а лишь потом выяснилось, что этот трюк можно вносить и под интегралы, и под дифференцирование, и под векторное произведение. Главная идея в том, с каждой неизвестной связана некоторая точка r_i (реально - семейство точек с весами, по которым идет численное интегрирование), а элемент матрицы имеет вид A_{ij} = f(r_j - r_i). Далее нужно существование универсальных разложений f(r_1+r_2)=g(r_1)h(r_2), а А f(x) = f(-x) вроде не требуется. По fast multipole method есть уйма статей в совершенно различных приложениях.
Вопросы по Вашей задаче:
1) У Вас метод конечных элементов или граничных элементов? В смысле интегралы по объему или по поверхности? Конечные элементы что-то мало в России встречаются.
2) Как выглядят базисные функции? В моем случае поверхность делится на треугольники, а область задания каждой функции - два смежных треугольника. Это так называемые RWG-функции. С другими не работал.
3) Где Вы все-таки этим занимаетесь? И под чьим руководством?
Плотные задачи тоже много решают итерационными методами, если переменных много - всяко быстрее, чем прямыми методами. Хотя есть задачи с плохой сходимостью... Не представляю, как отличить плохую сходимость от ошибок в коде.
Я так понимаю, есть два основных направления предобуславливателей. Первое - ILU. Но у него много всяких вариаций, параметров, алгоритмы из статьи в статью переходят с описками. Не всегда помогает, плюс использование специальных структур для хранения редких матриц добавляет проблем с быстродействием.
Второй, который идет от статьи Колотилиной, развился в SAI, sparse approximate inverse. В нем ищут какую-то матрицу с фиксированными местами (?) у элементов, минимизирующую некоторую норму, решая при этом (на этапе построения предобуславливателя) неслабую задачу наименьших квадратов. На этом во французском CERFACS собаку съели (J.Rahola, B.Carpentieri, I.Duff, L.Giraud). Но описания его я пока не добрался.
Есть еще блочный ILU, где матрицу делят на блоки, для блоков делают ILU, а матрицу целиком обращают рекурсивно более точно, используя обращенные блоки (J.Rius, A.Heldring).
С ILUT я связался, прочитав Jeonghwa Lee, Jun Zhang, Cai-Cheng Lu, "Incomplete LU preconditioning for large scale dense complex linear systems from electromagnetic wave scattering problems", Journal of Computational Physics,
Vol. 185, No. 1 (2003), pp. 158 - 175. У них есть и объемные, и поверхностные, и матрица несимметричная.
Но пока ни ILUT, ни ILU(p) ничего не дали.

[Certain]

Привет!
У вас изначально интегральные уравнения, поэтому и матрицы получаются плотными, я правильно понял? Сам я занимаюсь интегральными уравнениями, так что проблемы возникают схожие. Для решения своих СЛАУ использую GMRES, он позволяет достаточно быстро находить решение с точностью до $10^-^6$, что в моём случае и требуется. Матрицы комплексные, без диагонального преобладания, не сопряжённые по Эрмиту. Число итераций варьируется от 20 до 100 (в зависимости от исходных данных), размерность системы - от 500 до 16000. Так что переобусловливатели не нужны. Если есть вопросы, давайте пообщаемся.
Меня интересует, какой метод вы используете для дискретизации своих уравнений. Граничные элементы?

[Vassil]

Можете если хотите посмотреть их по адресу:

http://www.savefile.com/files2.php?fid=5377819

[KSergP]

СЛАУ у меня получается после линейного интегрирования (идея Харингтона).
Используемый вами подход более совершенный, но у лично у меня задача решать готовую СЛАУ.
Я аспирант Томского государственного университета систем управления и радиоэлектроники (ТУСУР). Руководители Газизов Т.Р.
Про QMR: я пробовал его реализовать, но что-то не получился пока. Там упоминается, что M=M1*M2 т.е. можно ли использовать M=LU (M1=L, M2=U)? Откуда вы взяли алгоритм?
Насчет метода: я пробовал BiCG, BiCGStab, GMRES лучше всех у меня BiCGStab. Я просто заменил действительную арифметику комплексной. В GMRES параметр m дает лишние проблемы надо подбирать его оптимальное значение. Реализовал CGS, но пока некогда тестировать.
Знакомы с TFQMR?

[Nizky]

To Vassil:
Уважаемый Васил Георгиев!
Большое спасибо за Ваш обзор! Несмотря на то, что численное решение интегральных уравнений итерационными методами начало развиваться только в середине 80х, в областях математики, существовавших до этого тоже есть интересные и полезные идеи. Кроме задач линейного программирования можно вспомнить еще и аналоговые вычисления (когда система уравнений решается мгновенно включением напряжения в сети, а решением является совокупность токов), методы предобуславливания в которых по-прежнему актуальны.

To Certain:
Метод дискретизации - граничные элементы, метод моментов. Базисные функции соответствуют каждой паре примыкающих треугольников. Описаны в S.M.Rao, D.R.Wilton, A.W.Glisson, "Electromagnetic scatering by surfaces of arbitrary shape", IEEE Transaction on Antennas and Propagation, Vol. 30, No. 3, pp. 409-418 (1982).
Из вопросов: как Вы выбираете параметр n в GMRES(n)? Считаете невязку внутри большого цикла (из n итераций) явно или неявно?

To KSerg:
у меня есть две статьи Roland W. Freund (вторая совместно с Noel M. Nachtigal):
"Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices", SIAM J. Sci. Stat. Comput., Vol. 13, No. 1, pp. 425-448 (1992)
"An implementation of the QMR method based on coupled two-term recurrences", SIAM J. Sci. Comput., Vol. 15, No. 2, pp. 313-337 (1994)
Я так понимаю, что в первой статье матричны должны быть симметричными, а во второй нет, но доступ к траспонированной матрице все равно нужен.
Есть очень хороший способ найти почти любые современные статьи - это сообщество "pdf" в Живом Журнале. Подробнее смотрите по ссылке http://www.livejournal.com/userinfo.bml?user=pdf (или просто www.livejournal.com). Но для этого, видимо, придется там зарегистрироваться.
С GMRES по моим наблюдениям всегда лучше брать большее n. Памяти мне не жалко, зато там возникают при больших n проблема: реальное значение невязки сильно разнится с оцениваемым. Ошибка, что ли, накапливается...

[Certain]

Nizky писал(а):

как Вы выбираете параметр n в GMRES(n)? Считаете невязку внутри большого цикла (из n итераций) явно или неявно?

Я провёл ряд экспериментов и выяснил, что число итераций слабо зависит от размерности СЛАУ. Иногда число итераций совсем не изменяется. Поэтому я заранее выставляю $n=150$ и мне хватает с запасом.
Невязку я считаю неявно, поскольку при реализации GMRES использую вращения Гивенса. Проблемы могут возникнуть, если Вы считаете с одинарной точностью, тогда возможны расхождения между явно и неявно вычисленной невязкой. По крайней мере у меня такое было. После перехода к вычислениям с двойной точностью проблема исчезла.
Интересно узнать, какова исходная дифференциальная задача, которую Вы решаете? Для рассеяния электромагнитных волн исходным уравнением, вероятно, является уравнение Максвелла? Как Вы сводите его к эквивалентной интегральной задаче?
С граничными элементами я только-только начал разбираться, поскольку использую для своего уравнения другой метод дискретизации. Если я правильно понимаю, граничные элементы всегда подразумевают триангуляцию поверхности. Это для меня важный вопрос, потому что используемый мною метод позволяет обходиться без триангуляции. Можно узнать об этом поподробнее? И ещё, как Вы вычисляете двойные интегралы: аналитически или по квадратурным формулам?

[Vassil]

Наверное не пропустили перелистать и следующие две книге, из которых хотя бы первая по моему прямо относится к кругу обсуждаемых проблем.

Дж. Ортега, Введение в паралельные и векторные методы решения линейних систем, Москва, "Мир", 1991

Уилкинсон, Райнш, Справочник алгоритмов на языке АЛГОЛ, Линейная алгебра, Москва, "Машиностроение", 1976

[Nizky]

Всем, кто интересуется итерационными методами решений систем уравнений, методами численного редения дифференциальных уравнений, предобуславливанием, редкими матрицами рекомендую книгу Yousef Saad "Iterative Methods for Sparse Linear Systems", 2е издание, 2000. Легко ищется через www.poiskknig.ru.

Что касается методов перехода от дифференциальных уравнений к интегральным, то рекомендую книгу К. Бреббия, Ж. Теллес, Л. Вроубел "Методы граничных элементов", М., "Мир", 1987.

В вывод уравнения, с которым работаю я (EFIE) особо не вникал. Набор преобразований понятен - интегрирование и интегрирование по частям, а уж способов их комбинировать слишком много.

О интегрировании по треугольникам. В поинтегральных функциях присутсвует 1/r, так что при малых расстояниях между треугольниками интегрируется по квадратурным функциям плохо. Мои предшественники посчитали аналитические интегралы для треугольников, которые касаются друг друга, и такие интегралы теперь считаются аналитически. Для очень далеких друг от друга треугольников интегралы считаются по одной точке, а для близких по специальным треугольным формулам Гаусса, которые состоят из 7 узлов (для каждого треугольника). Проверено, что если все считать по квадратурным формулам, то ошибка будет ощущаться в решении.

Тут я давно заметил странное поведение своего GMRES. Если через некоторое количество итераций модуль логарифма невязки растет все медленнее и медленнее, я считаю это признаком плохой обусловленности системы. Но бывает наоборот. Для некоторых "плохих" предобуславливателей GMRES запрягает очень медленно, после чего иногда ему удается раскочегариться, иногда после 40 и более итераций. После завершения n циклов все возращается на место - невязка как вкопанная и т.д.