Быстрое решение плотной линейной системы фиксированной размерности (N = 9), симметричной, положительно-полуопределенной

Какой алгоритм вы бы порекомендовали для быстрого решения плотной линейной системы фиксированной размерности (N = 9) (матрица симметрична, положительно-полуопределена)?

Гауссово исключениеLU разложениеРазложение холесскоготак далее?

Типы 32 и 64 бит с плавающей запятой.

Такие системы будут решаться миллионы раз, поэтому алгоритм должен быть достаточно быстрым по отношению к размерности (n = 9).

Постскриптум примерыкрепкий Реализация C ++ для предложенного алгоритма приветствуется.

1) Что вы подразумеваете под "решено миллион раз? Одна и та же матрица коэффициентов с миллионом различных правых слагаемых или миллион различных матриц?

Миллион различных матриц.

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

Ошибка повышения в порядке.

 qble13 нояб. 2012 г., 23:07
Спасибо за примечание
 biziclop13 нояб. 2012 г., 23:00
Ваш вопрос может быть лучшеscicomp.stackexchange.com
 Stefano M14 нояб. 2012 г., 00:10
Я бы определенно предложил опубликовать вопрос на scicomp: однако вы должны указать несколько дополнительных сведений. 1) Что вы подразумеваете под "решено миллион раз? Одна и та же матрица коэффициентов с миллионом различных правых слагаемых или миллион различных матриц? 2) Положительный _semi_definite означает, что матрица может быть единственной (для точности машины). Как бы вы хотели иметь дело с этим делом? Просто поднять ошибку или попытаться дать какой-нибудь толковый ответ?

Ответы на вопрос(6)

Как и другие выше, я рекомендую cholesky. Я'Мы обнаружили, что увеличение количества сложений, вычитаний и обращений к памяти означает, что LDLt медленнее, чем cholesky.

На самом деле существует ряд вариаций cholesky, и какая из них будет самой быстрой, зависит от представления, которое вы выбираете для своих матриц. Я обычно использую представление в стиле фортрана, то есть матрица M является двойной * M, где M (i, j) равно m [i + dim * j]; для этого я считаю, что верхний треугольный холеский (немного) самый быстрый, то есть верхний треугольный U с U '* U = M.

Для фиксированного, небольшого размера, безусловно, стоит подумать о написании версии, в которой нет циклов. Относительно простой способ сделать это - написать программу для этого. Насколько я помню, используя подпрограмму, которая работает с общим случаем в качестве шаблона, только утро заняло написание программы, которая написала бы конкретную версию с фиксированным измерением. Экономия может быть значительной. Например, моя общая версия занимает 0,47 секунды, чтобы выполнить миллион факторизаций 9x9, в то время как версия без петель занимает 0,17 секунды - эти временные интервалы работают на однопоточном компьютере с частотой 2,6 ГГц.

Чтобы показать, что это не главная задача, яНиже приведен источник такой программы. Он включает в себя общую версию факторизации в качестве комментария. Я'мы использовали этот код в обстоятельствах, когда матрицы не близки к единственному, и я считаю, что он хорошо работает там; однако это может быть слишком грубым для более деликатной работы.

/*  ----------------------------------------------------------------
**  to write fixed dimension ut cholesky routines
**  ----------------------------------------------------------------
*/
#include    
#include    
#include    
#include    
#include    
/*  ----------------------------------------------------------------
*/
#if 0
static  inline  double  vec_dot_1_1( int dim, const double* x, const double* y)
{   
double  d = 0.0;
    while( --dim >= 0)
    {   d += *x++ * *y++;
    }
    return d;
}

   /*   ----------------------------------------------------------------
    **  ut cholesky: solve U'*U = P for ut U in P (only ut of P accessed)
    **  ----------------------------------------------------------------
    */   

int mat_ut_cholesky( int dim, double* P)
{
int i, j;
double  d;
double* Ucoli;

    for( Ucoli=P, i=0; i
 qble14 нояб. 2012 г., 13:20
Спасибо! Я проверю.
 Stefano M14 нояб. 2012 г., 15:35
+1 за предоставленный код. Некоторый уровень развертывания цикла может быть выполнен компилятором за счет стоимости, иногда плохого распределения регистров. Внутренний цикл Chol содержит 720 умножений и добавляет: такая полная развертка возможна, как показывает ваш код, но я не уверен, что она будет лучше, чем тщательно закодированная версия цикла ...
 dmuir14 нояб. 2012 г., 18:05
я не был впечатлен GCC 'развернуть петлю. Если я возьму код измерения переменной, удалим аргумент dim и заменим его на 9, результирующая функция будет выполнена примерно за 0,43 микросекунды, тогда как для переменной с переменным один - 0,47, тогда как "handcoded» версия была 0.17
 David Doria04 апр. 2013 г., 20:44
Возможно, вы можете предоставить пример линейной системы для решения и как ввести ее в качестве матрицы для этого кода.

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

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

(Без выборки ваших матриц трудно дать здравый совет, но я подозреваю, что при таком небольшом порядке N = 9 поворот маленьких диагональных членов к нижней части D действительно не нужен. Поэтому я бы начал с классической Cholesky и просто прервать разложение на множители, если термины diag станут слишком малыми по отношению к некоторой разумно выбранной толерантности.)

Cholesky довольно прост в коде, и если вы стремитесь к действительно быстрому коду, лучше реализовать его самостоятельно.

 qble14 нояб. 2012 г., 09:49
Я волнуюсь, что "классический холецкий " LLT включает в себя N квадратных корней. Может быть, не проблема для большого N, но для 9 - это очень подозрительно.
 Stefano M14 нояб. 2012 г., 15:13
Вы можете использовать LDLT без поворота: вы сохраняете N квадратных корней за счет умножения (N-1) (N-2) / 2 во время факторизации: таким образом, это около 9 квадратных корней против 28 умножений. Это следует сравнить с основной стоимостью, которая составляет около (N ^ 3-N) / 6 Мэддс (madd = умножить и добавить), так что 720 Мэддс. Влияние на обратное решение составляет около N умножается. Если вы беспокоитесь об этих цифрах, нет выбора для реализации очень агрессивной собственной подпрограммы. Для таких задач вы должны оптимизировать на очень низком уровне, принимая во внимание пропускную способность памяти и задержку для различных уровней кэша, сродство процессора, ...

лучше всего использовать существующую библиотеку, а не подход «по-своему», поскольку есть много утомительных деталей, на которые следует обратить внимание в стремлении к быстрой и стабильной числовой реализации.

Вот'Вот несколько, чтобы начать:

Собственная библиотека (мои личные предпочтения):

http://eigen.tuxfamily.org/dox/QuickRefPage.html#QuickRef_Headers

Armadillo:http://arma.sourceforge.net/

Ищите вокруг и вынайду много других.

 linello14 нояб. 2012 г., 09:31
+1 за предложение Эйгена
 qble14 нояб. 2012 г., 08:20
Я уже рассматривал Эйген. Спасибо за пункт Armadillo, я проверю это.
 Stefano M14 нояб. 2012 г., 00:14
Хотя я согласен с тем, что библиотека в большинстве случаев является лучшим выбором, для такой небольшой матрицы издержки, связанные с полноценной библиотекой, могут быть чрезмерными.
 mtall14 нояб. 2012 г., 10:08
броненосец Это хороший выбор, поскольку он использует хорошо протестированный LAPACK для декомпозиции LU, декомпозиции Cholesky и т. д. Вы также можете легко изменить базовую библиотеку LAPACK с помощью высокооптимизированных заменителей, таких как Intel.Библиотеки МКЛ.
Решение Вопроса

Разложение холесского строго превосходит разложение LU. (примерно в два раза быстрее, чем LU, независимо от размера матрицы. Источник: "Численная линейная алгебра " Трефетен и Бау)

Это также де-факто стандарт для небольших плотных матриц (источник: я - доктор философии по вычислительной математике). Итерационные методы менее эффективны, чем прямые методы, если система не становится достаточно большой (быстрое правило, которое ничего не значит, но это всегда приятно иметь: на любом современном компьютере любая матрица меньше 100 * 100 - это определенно небольшая матрица, которая нуждается в прямых методах, а не в итерационных)

Теперь я не рекомендую делать это самостоятельно. Существует множество хороших библиотек, которые были тщательно протестированы. Но если я порекомендую вам один, это было быЭйген:

Установка не требуется (библиотека только для заголовков, поэтому нет библиотеки для ссылок, только #include <>)Надежный и эффективный (у них много тестов на главной странице, и результаты хорошие)Простой в использовании и хорошо документированный

Кстати,здесь в документацииу вас есть различные плюсы и минусы их 7 прямых линейных решателей в красивой, сжатой таблице. Похоже, что в вашем случае победит LDLT (разновидность Cholesky)

 qble14 нояб. 2012 г., 08:42
Может ты'придется самому это реализовать ... » - Я ожидал этого :) Поэтому мой главный вопрос был об алгоритме, а не о библиотеке.
 rkellerm09 сент. 2014 г., 14:21
@Fezvez - верно ли приведенное выше также для больших и разреженных матриц? (большее означает n: ~ 700- ~ 300 с ~ 1/16 матрицы, имеющей ненулевые значения)? Спасибо.
 Fezvez11 сент. 2014 г., 08:52
Извините, я немного заржавел в этом аспекте линейной алгебры. Я надеюсь, что вы найдете свой ответ (хотя, лучший способ - сравнить его самостоятельно)
 qble14 нояб. 2012 г., 08:28
Но если я порекомендую вам один, это будет Эйген " - Да, я рассмотрел это - это довольно круто, и имеет встроенные векторы фиксированного размера и матрицы. Но, к сожалению, документация говорит, что это не• работать на старых компиляторах (таких как GCC 3.x, VS2005 и т. д.). Хотя такое ограничение довольно разумно для очень высокопроизводительных библиотек, но проект, над которым я работаю, должен поддерживать такие старые компиляторы. "ЛПНП (вариация Холецкого) побеждает » - да, я проверил, это не такт требует sqrt по сравнению с LLT. Моя главная проблема была - может быть, LU лучше для маленьких матриц, таких как 9x9?
 Fezvez14 нояб. 2012 г., 08:35
Оооо, черт возьми, я никогда не думал о поддержке старых библиотек ... Ну, я неЯ не знаю библиотеку, которая бы соответствовала вашим потребностям, но я знаю, что когда матрица s.p.d, то Cholesky примерно в два раза быстрее, чем LU (источник: "Численная линейная алгебра " Трефетен и Бау)s независимо от размера. Может ты'придется самому это реализовать ... »

особенно еслирешено миллионы раз действительно означает "решена однажды и применена к миллионам векторов, Вы'создам декомпозицию LU, сохраняем ее и применяем подстановку вперед-назад против стольких r.h.s. векторы, как вы хотите.

Это'Более стабильный перед округлением, если вы используете поворот.

 qble14 нояб. 2012 г., 08:02
Нет, извините, я имел в виду миллионы разных систем Ax = b (где A тоже отличается). В любом случае +1

Из-за простоты использования вы можете взять решатели Eigen только для сравнения. Для конкретного случая использования конкретный решатель может быть быстрее, хотя другой должен быть лучше. Для этого вы можете измерить время выполнения для каждого алгоритма только для выбора. После этого вы можете реализовать желаемый вариант (или найти существующий, который наилучшим образом соответствует вашим потребностям).

Ваш ответ на вопрос