Повторное решение

12

Я использую MATLAB для решения проблемы, которая включает в себя решение на каждом временном шаге, где b изменяется со временем. Прямо сейчас я делаю это, используя MATLAB :Ax=bbmldivide

x = A\b

У меня есть возможность делать столько предварительных вычислений, сколько нужно, поэтому мне интересно, есть ли более быстрый и / или более точный метод, чем mldivide. Что обычно делается здесь? Спасибо всем!

сомнение
источник
1
У вас есть конкретные знания о структуре ? Например, это симметрично? Положительно определен? Трехдиагональная? Ортогональные? A
Доминик
Матрица является плотной квадратной матрицей. A
Сомнение
3
Если у вас нет других знаний об , то лучшим вариантом будет факторизация L U, как описано в ответе ниже. ALU
Доминик

Ответы:

14

Самое очевидное, что вы можете сделать, это предварительно вычислить

[L,U] = lu(A) ~ O (n ^ 3)

Тогда вы просто вычислите

x = U \ (L \ b) ~ O (2 n ^ 2)

Это значительно сократит стоимость и сделает это быстрее. Точность была бы такой же.

Milind R
источник
1
Обратите внимание, из документации , L не обязательно является нижней треугольной. Этот ответ, скорее всего, будет быстрее, чем прямое решение, однако я бы позаботился о том, чтобы команда L \ b была достаточно умной, чтобы знать, чтобы решить L в правильном порядке (скорее всего, но это не говорит наверняка в документации).
Годрик Провидец
Да, вы правы, L является произведением нижней треугольной матрицы и матрицы перестановок. Но я буду проклят, если он не признает, что все, что ему нужно, это обратная замена L\b. Потому что я видел, как именно эту строку используют в высокопроизводительном коде те, кого я считаю экспертами.
Milind R
8
O(n2)
1
A
3
@ BrianBorcher Насколько я знаю, лучший способ отслеживать перестановку - [L,U,p] = lu(A,'vector'); x = U\(L\b(p));см. Пример 3 в lu документации .
Стефано М
5

В наших научных компьютерных курсах по этой теме мы провели несколько обширных компьютерных классов. Для «маленьких» вычислений, которые мы там делали, оператор обратной косой черты в Matlab всегда был быстрее, чем что-либо еще, даже после того, как мы максимально оптимизировали наш код и предварительно переупорядочили все матрицы (например, с помощью обратного порядка Cuthill McKee для разреженных матриц) ,

Вы можете проверить одну из наших лабораторных инструкций . Ответ на ваш вопрос описан (в ближайшее время) на стр. 4.

Хорошая книга на эту тему написана, например, Чейни .

СЕБ
источник
4

An×n Axi=bii=1mm

V = inv(A);
...
x = V*b;

O(n3)inv(A)O(n2)V*bm

>> n = 5000;
>> A = randn(n,n);
>> x = randn(n,1);
>> b = A*x;
>> rcond(A)
ans =
   1.3837e-06
>> tic, xm = A\b; toc
Elapsed time is 1.907102 seconds.
>> tic, [L,U] = lu(A); toc
Elapsed time is 1.818247 seconds.
>> tic, xl = U\(L\b); toc
Elapsed time is 0.399051 seconds.
>> tic, [L,U,p] = lu(A,'vector'); toc
Elapsed time is 1.581756 seconds.
>> tic, xp = U\(L\b(p)); toc
Elapsed time is 0.060203 seconds.
>> tic, V=inv(A); toc
Elapsed time is 7.614582 seconds.
>> tic, xv = V*b; toc     
Elapsed time is 0.011499 seconds.
>> [norm(xm-x), norm(xp-x), norm(xl-x), norm(xv-x)] ./ norm(x)
ans =
   1.0e-11 *
    0.1912    0.1912    0.1912    0.6183

A1LUm>125

Некоторые заметки

Для анализа стабильности и ошибок см. Комментарии к этому другому ответу , особенно к VictorLiu.

mn

Синхронизация была выполнена с Matlab R2011b на 12-ядерном компьютере с довольно постоянной средней нагрузкой UNIX 5; лучшее tic, tocвремя из трех зондов.

Стефано М
источник
Действительно, в матричном умножении намного больше параллелизма, чем в треугольном решателе, так что это должно быть еще более очевидным, если вычисления выполняются параллельно (многоядерный / графический процессор / и т. Д.) Любым способом.
Арон Ахмадиа
@AronAhmadia Я согласен: оценки точки безубыточности, основанные только на количестве операций, имеют смысл только для последовательной реализации.
Стефано М
1
Обратите внимание, что если матрица A будет разреженной, ситуация будет сильно отличаться - обратное будет, как правило, довольно плотным, в то время как коэффициенты LU, как правило, достаточно малы, что приводит к ускорению движения назад в направлении LU.
Брайан Борхерс
1
A
1
inv(A)Ax=bbBA\B
2

Взгляните на этот вопрос , ответы показывают, что mldivideэто довольно умно, а также дает советы о том, как увидеть, что Matlab использует для решения A\b. Это может дать вам подсказку относительно параметров оптимизации.

Psirus
источник
0

Использование обратной косой черты более или менее эквивалентно тому inv(A)*B, что если вы кодируете ее свободно, последняя может быть более интуитивной. Они примерно одинаковы (только отличаются в том, как выполняются вычисления), хотя вы должны проверить документацию Matlab для пояснения.

Чтобы ответить на ваш вопрос, обратная косая черта, как правило, хорошо, но она зависит от свойств матрицы масс.

AlanH
источник
1
Математически inv (A) * b - это то же самое, что и \ \ \ n \ n \ n \ n \ n в численном отношении, фактически формирование обратного является менее эффективным и менее точным. Если вы работаете над изучением линейной алгебры, это может быть приемлемым, но я бы сказал, что вам нужна очень веская причина для формирования обратного.
Годрик Провидец
Но зачем вам когда-нибудь вычислять, inv(A)так как это дороже A\b?
Доминик
7
@Godric: недавно появилась статья, в которой обсуждается «миф» о том, что inv (A) * b менее точен: об ArXiv . Не говоря о том, что обычно есть причина вычислять реальное обратное, а просто сказать.
Виктор Лю
3
@Dominique: Треугольные решения гораздо менее распараллеливаемы, чем умножение матрицы на вектор, и сложные предварительно подготовленные итерационные методы часто используют прямые методы на поддоменах. Часто полезно явно сформировать инверсии нескольких плотных треугольных матриц скромного размера, чтобы улучшить параллелизм.
Джек Поулсон
@VictorLiu: Спасибо за статью. Я исправлен в своем заявлении о точности (по крайней мере, для умных реализаций inv (A)).
Годрик Провидец