Как ефективно да се реши линейна система с Лапласиан + диагонална матрица?

В моята реализация на алгоритъм за обработка на изображения, трябва да реша голяма линейна система от формата A*x=b, където:

  • Матрица A=L+D е сумата от лапласова матрица L и диагонална матрица D
  • Лапласианската матрица L е разредена, с около 25 ненули на ред
  • Системата е голяма, с толкова неизвестни, колкото има пиксели във входното изображение (обикновено > 1 милион).

Матрицата на Лаплас L не се променя между последователните изпълнения на алгоритъма; Мога да конструирам тази матрица при предварителна обработка и евентуално да изчисля нейното факторизиране. Диагоналната матрица D и десният вектор b се променят при всяко изпълнение на алгоритъма.

Опитвам се да разбера кой би бил най-бързият метод за решаване на системата по време на изпълнение; Нямам нищо против да отделям време за предварителна обработка (за изчисляване на факторизация на L, например).

Първоначалната ми идея беше предварително да изчисля факторизация на Cholesky на L, след това да актуализирам факторизацията по време на изпълнение със стойности от D (ранг-1 актуализация с cholupdate) и да разреша бързо проблема с обратно заместване. За съжаление факторизацията на Cholesky не е толкова рядка като оригиналната L матрица и самото й зареждане от диска вече отнема 5,48 s; за сравнение, отнема 8,30 секунди за директно решаване на системата с обратна наклонена черта.

Като се има предвид формата на моите матрици, има ли друг метод, който бихте препоръчали за ускоряване на решаването по време на изпълнение, без значение колко време отнема при предварителната обработка?


matlab sparse-matrix system factorization
person Asmo    schedule 29.01.2013    source източник
comment
Опитайте SVD, но не знам колко е бърз за толкова голямо количество данни...   -  person Ander Biguri    schedule 29.01.2013
comment
Не виждам защо прилагането на SVD тук би било добра идея, предвид формата на моите матрици.   -  person Asmo    schedule 30.01.2013
comment
да, вероятно си прав, не съм експерт. Само в случай, че не си го взел предвид :)   -  person Ander Biguri    schedule 31.01.2013
comment
Разгледайте това: scicomp.stackexchange.com/questions/1001/   -  person johnish    schedule 31.01.2013

Отговори (1)


arrow_upward
2
arrow_downward

Ако приемем, че работите върху мрежа (тъй като споменавате изображения - въпреки че това не е гарантирано), че се интересувате повече от скоростта, отколкото от прецизността (тъй като 5s вече изглежда твърде бавно за 1 милион неизвестни), виждам няколко опции.

Първо, забравете за точните методи като Cholesky (+пренареждане). Дори ако те позволяват да съхранявате факторизацията и да я използвате повторно за множество rhs, вероятно ще трябва да съхранявате гигантски матрици, които изглеждат неразрешими във вашия случай (надявам се, че пренареждате редове/колони с обратен Cuthill McKee или нещо подобно в противен случай обаче - това разрежда факторизацията много).

В зависимост от вашите гранични условия, първо бих опитал Matlab poisolv, който решава проблем на Поасон с помощта на FFT, и възможни препроекции, ако искате гранични условия на Дирихле вместо периодични. Много е бърз, но може да не е подходящ за вашия проблем (споменавате, че имате 25 nnz за лапласова матрица+идентичност: защо? това е матрица на Лаплас от висок ред, в който случай може да се интересувате повече от точността, отколкото това, което предполагам ? или всъщност това е различен проблем от този, който описвате?).

След това можете да опитате многомрежови решаващи програми, които са много бързи за изображения и гладки проблеми. Можете да използвате прост метод за релаксация за всяка итерация и всяко ниво на мултирешетката или да използвате по-фантастични методи (например предварително кондиционирано ниво на спрегнат градиент). Като алтернатива можете да направите по-прост предварително кондициониран конюгиран градиент (или дори SSOR) без многорешетка и ако се интересувате само от приблизително решение, можете да спрете итерациите преди пълната конвергенция.

Моите аргументи за итеративни решаващи програми са:

  • можете да спрете преди конвергенцията, ако искате приблизителен проблем
  • все още можете да използвате повторно други резултати, за да инициализирате вашето решение (например, ако вашите различни изпълнения съответстват на различни кадри от видео, тогава използването на решението на предишния кадър като инициализация на следващия би имало известен смисъл).

Разбира се, директен решаващ инструмент, за който можете предварително да изчислите, съхраните и запазите факторизацията, също има смисъл (въпреки че не разбирам аргумента ви за актуализация на ранг-1, ако вашата матрица е постоянна), тъй като остава да се направи само обратното заместване на време на изпълнение. Но като се има предвид, че това игнорира структурата на проблема (правилна решетка, възможен интерес към резултати с ограничена точност и т.н.), бих избрал методи, които са проектирани за тези случаи, като методи, подобни на Фурие, или многомрежови. И двата метода могат да бъдат приложени на GPU за по-бързи резултати (припомнете си, че GPU са по-скоро пригодени за работа с изображения/текстури!).

И накрая, можете да получите интересни отговори от scicomp.stackexchange, който е по-насочен към числения анализ.

person nbonneel    schedule 01.02.2013
comment
+1 за споменаването на пренареждането! Просто бих добавил непълни факторизации на Cholesky към вашия списък, перфектно иначе! - person Dr_Sam; 04.02.2013
comment
благодаря :) Бях с впечатлението, че непълните факторизации на Cholesky за прекондициониране са много по-бавни от обикновен прекондиционер на Якоби само за незначително подобрение в скоростта на конвергенция; не е ли - person nbonneel; 06.02.2013