Как оператор обратной косой черты в MATLAB решает

36

Я сравнивал несколько своих кодов с «стандартными» кодами MATLAB. Я удивлен результатами.

Я запустил пример кода (разреженная матрица)

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

Полученные результаты :

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

Для плотной матрицы:

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

Полученные результаты:

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

Как, черт возьми, это так круто?

расследование
источник
1
Встроенная обратная косая черта в MATLAB, другими словами, прямой решатель для системы линейных уравнений, использует многофронтальный метод для разреженной матрицы, поэтому A \ B является настолько классным.
Шухао Цао
1
Он использует код Тима Дэвиса, доступный по адресу cise.ufl.edu/research/sparse . Кроме того, удивительность уходит, когда у вас нетривиальная проблема
Стали
1
Что такое "LULU"? Как вы думаете, почему это хорошая реализация факторизации LU и последующего прямого решения?
Джед Браун
3
@Nunoxic Какая реализация? Ты сам написал? Высокопроизводительная плотная линейная алгебра, хотя обычно и хорошо понятна алгоритмически, нелегко реализовать эффективно на современном оборудовании. Лучшие реализации BLAS / Lapack должны быть близки к пиковым значениям для матрицы такого размера. Кроме того, из ваших комментариев у меня складывается впечатление, что вы думаете, что LU и метод исключения Гаусса - это разные алгоритмы.
Джед Браун
1
Он вызывает код Fortran, написанный с использованием Intel MKL.
расследование

Ответы:

37

В Matlab команда '\' вызывает алгоритм, который зависит от структуры матрицы A и включает проверки (небольшие накладные расходы) на свойства A.

  1. Если A разреженный и полосатый, используйте решатель полос.
  2. Если A - верхняя или нижняя треугольная матрица, используйте алгоритм обратной подстановки.
  3. Если A симметрична и имеет реальные положительные диагональные элементы, попробуйте факторизацию Холецкого. Если A мало, сначала используйте переупорядочение, чтобы минимизировать заполнение.
  4. Если ни один из вышеуказанных критериев не выполняется, выполните общую треугольную факторизацию с использованием исключения Гаусса с частичным поворотом.
  5. Если A немного, используйте библиотеку UMFPACK.
  6. Если A не является квадратом, используйте алгоритмы, основанные на факторизации QR для неопределенных систем.

Чтобы уменьшить накладные расходы, можно использовать команду linsolve в Matlab и выбрать подходящий решатель среди этих вариантов самостоятельно.

Аллан П. Энгсиг-Каруп
источник
Предполагая, что я имею дело с неструктурированной плотной матрицей размером 10000x10000 со всеми элементами, отличными от нуля (высокий уровень плотности), какой будет мой лучший выбор? Я хочу выделить тот 1 алгоритм, который работает для плотных матриц. Это LU, QR или Гауссово исключение?
Дознание
1
Похоже на Шаг 4, где вызывается удаление по Гауссу, что соответствует наиболее общему случаю, когда структура А не может быть использована для повышения производительности. Таким образом, в основном это факторизация LU и последующий шаг вперед, за которым следует шаг обратной замены.
Аллан П. Энгсиг-Каруп
Благодарность! Я думаю, что это дает мне направление думать. В настоящее время устранение Гаусса - лучшее, что у нас есть для решения таких неструктурированных проблем, верно?
Дознание
37

Если вы хотите увидеть, что a\bподходит для вашей конкретной матрицы, вы можете установить spparms('spumoni',1)и точно определить , какой алгоритм вас впечатлил. Например:

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

будет выводить

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

так что я вижу, что "\" в конечном итоге использовало "CHOLMOD" в этом случае.

dranxo
источник
3
+1 за новые настройки MATLAB, о которых я никогда не слышал. Хорошо сыграно, сэр.
Джефф Оксберри
2
Эй спасибо! Это в help mldivide.
dranxo
16

Для разреженных матриц Matlab использует UMFPACK для \операции " ", которая, в вашем примере, в основном пробегает значения a, инвертирует их и умножает их на значения b. Для этого примера, однако, вы должны использовать b./diag(a), что намного быстрее.

Для плотных систем оператор обратной косой черты немного сложнее. Краткое описание того, что делается, когда дается здесь . Согласно этому описанию, в вашем примере Matlab решит, a\bиспользуя обратную подстановку. Для общих квадратных матриц используется LU-разложение.

Pedro
источник
Заказной. Разреженность, inv матрицы diag будет просто обратной величиной диагональных элементов, так что b./diag(a) будет работать, но a \ b отлично работает и для общих разреженных матриц. Почему Linsolve или LULU (моя оптимизированная версия LU) не быстрее, чем a \ b для плотных матриц в этом случае.
Дознание
@Nunoxic Делает ли ваш LULU что-нибудь, чтобы обнаружить диагональность или треугольность плотной матрицы? Если нет, то это займет столько же времени для каждой матрицы, независимо от ее содержимого или структуры.
Педро
В некотором роде. Но с помощью флагов linsolve OPT я определил все, что нужно определить о структуре. Тем не менее, a \ b быстрее.
Дознание
@Nunoxic, каждый раз, когда вы вызываете пользовательскую функцию, вы вызываете накладные расходы. Matlab делает все для внутренней обратной косой черты, например, для разложения и последующего применения правой стороны, с очень небольшими накладными расходами, и, следовательно, будет быстрее. Кроме того, в ваших тестах вы должны использовать более одного вызова для получения надежных таймингов, например tic; for k=1:100, a\b; end; toc.
Педро
5

Как правило, если у вас есть разреженная матрица разумной сложности (т. Е. Это не обязательно должен быть шаблон из 5 точек, но на самом деле это может быть дискретизация уравнений Стокса, для которых число ненулевых элементов в строке равно намного больше, чем 5), тогда разреженный прямой решатель, такой как UMFPACK, обычно побеждает итеративный решатель Крылова, если проблема не превышает приблизительно 100 000 неизвестных.

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

Вольфганг Бангерт
источник
3
Мне не ясно, как это отвечает на вопрос, но я также не согласен с предпосылкой. Это правда, что прямые решатели, как правило, хорошо работают для 2D задач небольшого размера, но прямые решатели, как правило, страдают в 3D задолго до неизвестных 100k (разделители вершин намного больше, чем в 2D). Кроме того, я утверждаю, что в большинстве случаев (например, эллиптические операторы) итерационные решатели могут быть созданы быстрее, чем прямые решатели, даже для 2D-задач среднего размера, но для этого может потребоваться значительное усилие (например, использование правильных ингредиентов для обеспечения работы многосетки) ,
Джед Браун
1
Если ваши проблемы достаточно малы, и вы не хотите думать о разработке неявных решателей, или если ваши проблемы очень сложны (например, высокочастотный Максвелл), и вы не хотите посвятить свою карьеру разработке хороших решателей, тогда я согласен, что разреженные прямые решатели - отличный выбор.
Джед Браун
На самом деле моя проблема заключается в разработке хорошего плотного решателя. У меня нет конкретного приложения как такового (следовательно, нет структуры). Я хотел настроить несколько своих кодов и сделать их конкурентоспособными с тем, что используется в настоящее время. Мне было просто любопытно, что такое a \ b и как он всегда выбивает из моего кода дерьмо.
расследование
@JedBrown: Да, может быть, если приложить немало усилий, можно найти итеративные решатели, которые побеждают прямой решатель даже для небольших двумерных задач. Но зачем это делать? Для задач с <100k неизвестными достаточно разреженных прямых решателей в 2d. Я хотел сказать, что для решения таких небольших проблем не тратьте время на то, чтобы поработать и выяснить наилучшую комбинацию параметров - просто возьмите черный ящик.
Вольфганг Бангерт
Уже существует порядок различий во времени выполнения между разреженной геометрической многосеткой Холецкого и (на основе матрицы) для «простых» задач 2D с 100k степеней свободы с использованием 5-точечного трафарета (~ 0,5 секунды по сравнению с ~ 0,05 секунды). Если ваш трафарет использует второго соседа (например, дискретизация четвертого порядка, Ньютон для некоторых вариантов нелинейной реологии, стабилизации и т. Д.), Размер разделителей вершин примерно удваивается, поэтому стоимость прямого решения возрастает примерно в 8 раз (стоимость больше проблемный для MG). Со многими временными шагами или циклами оптимизации / UQ эти различия значительны.
Джед Браун