Когда ортогональные преобразования превосходят исключение Гаусса?

Как мы знаем, методы ортогональных преобразований (повороты Гивенса и отражения Хаусхолдера) для систем линейных уравнений более дороги, чем устранение по Гауссу, но теоретически обладают более хорошими свойствами устойчивости в том смысле, что они не изменяют число условий системы. Хотя я знаю только один академический пример матрицы, которая испорчена исключением Гаусса с частичным поворотом. И есть общее мнение, что на практике такого рода поведение вряд ли встретится (см. Примечания к этой лекции [pdf] ).

Итак, где мы будем искать ответ по теме? Параллельные реализации? Обновление? ..

точность

Трефетен и Шрайбер написали отличную статью « Стабильность по Гауссу в среднем случае» , в которой обсуждается вопрос о точности вашего вопроса. Вот несколько его выводов:

1. «Для QR разложения с или без колонки поворота, средний максимальный элемент остаточной матрицы , тогда как для исключения Гаусса это O ( п ) . Это сравнение показывает , что исключение Гаусса умеренно неустойчивым, но нестабильность может быть обнаружена только для очень больших матричных задач, решаемых с низкой точностью. Для большинства практических задач элиминация по Гауссу в среднем очень стабильна. "(Выделение мое)$O(n^{1/2})$$O(n)$

2. «После первых нескольких этапов исключения из Гаусса остальные матричные элементы примерно нормально распределены, независимо от того, были ли они начаты таким образом».

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

Спектакль

Для квадратных вещественных матриц, LU с частичной поворота требуется примерно провалов, тогда как Хаусхолдера на основе QR - требуется примерно 4 / 3 л $2/3 n^3$$4/3 n^3$ провалов. Таким образом, для достаточно больших квадратных матриц QR-факторизация будет в два раза дороже, чем LU-факторизация.

Для матриц, где т ≥ п , LU с частичным поворотом требует м п 2 - п 3 / 3 флопа, по сравнению с QR -- х 2 м н 2 - 2 л 3 / 3 (который по - прежнему вдвое больше , чем LU факторизации). Тем не менее , неожиданно часто приложения создают очень высокие узкие матрицы ( m ≫ $m \times n$$m \ge n$$mn^2 - n^3/3$$2mn^2 - 2n^3/3$$m \gg n$ ), и Demmel et al. приятной работы, избегая общения параллельной и последовательной QR-факторизации, который (в разделе 4) обсуждает умный алгоритм, который требует, чтобы сообщения отправлялись только при использовании p процессоров, в отличие от n log p сообщений традиционных подходов. Суть в том, что O ( n 3 log p ) выполняются дополнительные провалы, но для очень малых n это часто предпочтительнее, чем стоимость задержки отправки большего количества сообщений (по крайней мере, когда необходимо выполнить только одну факторизацию QR).$\log p$$p$$n \log p$$O(n^3 \log p)$$n$

Я удивлен, что никто не упомянул линейные задачи наименьших квадратов , которые часто встречаются в научных вычислениях. Если вы хотите использовать метод исключения Гаусса, вы должны сформировать и решить нормальные уравнения, которые выглядят так:

где - матрица точек данных, соответствующих наблюдениям независимых переменных, x - вектор параметров, которые нужно найти, и $A$$x$$b$ - вектор точек данных, соответствующих наблюдениям зависимой переменной.

Как часто отмечает Джек Полсон, число условий является квадратом числа условий A , поэтому нормальные уравнения могут быть катастрофически плохо обусловлены. В таких случаях, хотя QR-и SVD-подходы медленнее, они дают гораздо более точные результаты.$A^{T}A$$A$

Как вы измеряете производительность? Скорость? Точность? Стабильность? Быстрый тест в Matlab дает следующее:

>> N = 100;
>> A = randn(N); b = randn(N,1);
>> tic, for k=1:10000, [L,U,p] = lu(A,'vector'); x = U\(L\b(p)); end; norm(A*x-b), toc
ans =
   1.4303e-13
Elapsed time is 2.232487 seconds.
>> tic, for k=1:10000, [Q,R] = qr(A); x = R\(Q'*b); end; norm(A*x-b), toc             
ans =
   5.0311e-14
Elapsed time is 7.563242 seconds.

Таким образом, решение одной системы с LU-декомпозицией примерно в три раза быстрее, чем ее решение с помощью QR-декомпозиции, за счет половины точности до десятичного знака (этот пример!).

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

Рассмотрим пространство всех матриц. Эти методы прекрасно работают практически везде. То есть 99,999 ...% всех матриц, которые можно создать, не будут иметь проблем с нестабильными методами. Существует только очень небольшая часть матриц, для которых GE и другие будут испытывать трудности.

Проблемы, которые волнуют исследователей, как правило, в этой небольшой части.

Мы не строим матрицы случайно. Мы строим матрицы с очень особыми свойствами, которые соответствуют совершенно особым неслучайным системам. Эти матрицы часто плохо обусловлены.

Геометрически вы можете рассмотреть линейное пространство всех матриц. Существует подпространство нулевого объема / меры сингулярных матриц, прорезающих это пространство. Многие проблемы, которые мы строим, сгруппированы вокруг этого подпространства. Они не распределены случайно.

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