Итерационные методы для линейных систем

Одно из большинства важных и распространенных приложений числовой линейной алгебры является решением линейных систем, которые могут быть выражены в форме A*x = b. Когда A большая разреженная матрица, можно решить линейную систему с помощью итерационных методов, которые включают вам к компромиссу между временем выполнения вычисления и точностью решения. Эта тема описывает итерационные методы, доступные в MATLAB^®, чтобы решить уравнение A*x = b.

Прямой по сравнению с итерационными методами

Существует два типа методов для решения линейных уравнений A*x = b:

Прямые методы являются вариантами Исключения Гаусса. Эти методы включают отдельные элементы матрицы непосредственно, посредством операций над матрицей, таких как LU, QR или факторизация Холесского. Можно использовать прямые методы решить линейные уравнения с высоким уровнем точности, но эти методы могут быть медленными при работе с большими разреженными матрицами. Скорость решения линейной системы с прямым методом строго зависит от размера матрицы коэффициентов.
Например, этот код решает маленькую линейную систему.
```
A = magic(5);
b = sum(A,2);
x = A\b;
norm(A*x-b)
```
```
ans =

   1.4211e-14
```
MATLAB реализует прямые методы через матричные операторы деления / и \, а также функционирует, такие как lsqminnormразложение, и linsolve.
Итерационные методы производят приближенное решение линейной системы после конечного числа шагов. Эти методы полезны для больших систем уравнений, где это разумно к точности компромисса в течение более короткого времени выполнения. Эти методы включают матрицу коэффициентов только косвенно посредством матричного векторного произведения или абстрактного линейного оператора. Итерационные методы обычно применяются только к разреженным матрицам, потому что меньшие системы могут быть легко решены с прямыми методами. Скорость решения линейной системы с косвенным методом не зависит так же строго от размера матрицы коэффициентов как прямой метод. Однако использование итерационного метода обычно требует настраивающихся параметров для каждой определенной проблемы.
Например, этот код решает большую разреженную линейную систему, которая имеет симметричную положительную определенную матрицу коэффициентов.
```
A = delsq(numgrid('L',400));
b = ones(size(A,1),1);
x = pcg(A,b,[],1000);
norm(A*x-b)
```
```
pcg converged at iteration 796 to a solution with relative residual 9.9e-07.

ans =

   3.4285e-04
```
MATLAB реализует множество итерационных методов, которые имеют различные достоинства и недостатки в зависимости от свойств матрицы коэффициентов A.

Прямые методы обычно быстрее и в более общем плане применимы, чем косвенные методы, если существует достаточно устройства хранения данных, доступного, чтобы выполнить их. Обычно необходимо попытаться использовать x = A\b сначала. Если прямые решают, является слишком медленным, то можно попытаться использовать итерационные методы.

Типовой итеративный алгоритм

Большинство итеративных алгоритмов, которые решают линейные уравнения, следует за подобным процессом:

Начните с исходного предположения для вектора решения x0. (Это обычно - нулевой вектор, если вы не задаете лучшее предположение.)
Вычислите норму невязки res = norm(b-A*x0).
Сравните невязку с заданным допуском. Если res <= tol, закончите расчет и дайте вычисленный ответ для x0.
Обновите величину и направление векторного x0 на основе значения невязки и других расчетных количеств.
Повторите Шаги 2 - 4 до значения x0 достаточно хорошо, чтобы удовлетворить допуску.

Итерационные методы отличаются по тому, как они обновляют величину и направление x0 на Шаге 4, и у некоторых есть немного отличающиеся критерии сходимости на Шагах 2 и 3, но это получает базовый процесс, за которым следуют все итеративные решатели.

Сводные данные итерационных методов

MATLAB имеет несколько функций, которые реализуют итерационные методы для разреженных систем линейных уравнений. Эти методы спроектированы, чтобы решить A x = b или минимизировать норму || b – A x ||. Несколько из этих методов имеют общие черты и основаны на тех же базовых алгоритмах, но каждый алгоритм обладает преимуществами в определенных ситуациях [1], [2].

Описание	Примечания
`lsqr` Метод наименьших квадратов	Матрица коэффициентов может быть прямоугольной. Единственный решатель, доступный для прямоугольных систем. Аналитически эквивалентный методу сопряженных градиентов (PCG) применился к нормальным уравнениям `(A'A)x = A'*b`, но обладает более благоприятными числовыми свойствами.
`pcg` (предобусловленные методы сопряженных градиентов)	Матрица коэффициентов должна быть симметрична положительный определенный. Самый эффективный решатель для симметричных положительных определенных систем, поскольку устройство хранения данных только для ограниченного количества векторов требуется.
`minres` (минимальная невязка)	Матрица коэффициентов должна быть симметричной, но не должна быть положительна определенный. Каждая итерация минимизирует остаточную ошибку в 2-норме, таким образом, алгоритм, как гарантируют, сделает успехи от шага до шага. Не страдает от отказов (когда алгоритм становится не могущим сделать успехи к решению и остановам).
`symmlq` (симметричный LQ)	Матрица коэффициентов должна быть симметричной, но не должна быть положительна определенный. Решает спроектированную систему и сохраняет невязку ортогональной ко всем предыдущим единицам. Не страдает от отказов (когда алгоритм становится не могущим сделать успехи к решению и остановам).
`bicg` (бисопряженный градиент)	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. `bicg` является в вычислительном отношении дешевым, но сходимость неправильна и ненадежна. `bicg` редко лучший выбор для большинства проблем, но это исторически важно, потому что много других алгоритмов были разработаны как улучшения его.
`bicgstab` (бисопряженный стабилизированный градиент)	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. Использование шаги BiCG, чередующиеся с GMRES (1) шаги для дополнительной устойчивости.
`bicgstabl` (бисопряженный градиент стабилизируется (l)),	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. Использование шаги BiCG, чередующиеся с GMRES (2) шаги для дополнительной устойчивости.
`cgs` (метод сопряженных градиентов придал квадратную форму),	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. Требует того же количества операций на итерацию как `bicg`, но избегает использования транспонирования путем работы с квадратом остатка.
`gmres` (обобщенная минимальная невязка)	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. Один из самых надежных алгоритмов, поскольку норма невязки минимизирована в каждой итерации. Работа и устройство хранения данных требуемое повышение линейно с количеством итерации. Выбор соответствующего `restart` значение важно, чтобы избежать ненужной работы и устройства хранения данных.
`qmr` (квазиминимальная невязка)	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. Наверху на итерацию немного больше, чем `bicg`, но это обеспечивает больше устойчивости.
`tfqmr` (квазиминимальная невязка без транспонирования)	Матрица коэффициентов должна быть квадратной, но не должна быть симметричной. Лучший решатель, чтобы попробовать за симметричные неопределенные системы, когда память ограничивается.

Выбор итеративного решателя

Эта блок-схема итеративных решателей в MATLAB дает общее представление о ситуациях, где каждый решатель полезен. Можно обычно использовать gmres почти для всех квадратных, несимметричных проблем. Существуют некоторые случаи где бисопряженные алгоритмы градиентов (bicgbicgstabcgs, и так далее), более эффективны, чем gmres, но их непредсказуемое поведение сходимости часто делает gmres лучший начальный выбор.

Предварительные формирователи

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

$A x = b$

в эквивалентную систему

$\tilde{A} \tilde{x} = \tilde{b} .$

Идеальный предварительный формирователь преобразовывает матрицу коэффициентов A в единичную матрицу, поскольку любой итерационный метод будет сходиться в одной итерации с таким предварительным формирователем. На практике нахождение хорошего предварительного формирователя требует компромиссов. Преобразование выполняется одним из трех способов: оставленное предварительное создание условий, правильное предварительное создание условий или предварительное создание условий разделения.

Первый случай называется оставленным предварительным созданием условий начиная с матрицы перед формирователем, M появляется слева от A:

$(M^{- 1} A) x = (M^{- 1} b) .$

Эти итеративные решатели использование оставили предварительное создание условий:

В правильном предварительном создании условий M появляется справа от A:

$(A M^{- 1}) (M x) = b .$

Эти итеративные решатели используют правильное предварительное создание условий:

Наконец, для симметричных содействующих матриц, разделяет, предварительное создание условий гарантирует, что новая система все еще имеет симметрическую матрицу. Предварительный формирователь $M = H H^{T}$ появляется с обеих сторон A:

$(H^{- 1} A H^{- T}) H^{T} x = (H^{- 1} b)$

Алгоритм решателя для разделения, которое предобусловленные системы на основе вышеупомянутого уравнения, но на практике нет никакой потребности вычислить H. Алгоритм решателя только когда-либо применяет M или inv(M) непосредственно.

Эти итеративные решатели используют предварительное создание условий разделения:

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

Итеративные решатели в MATLAB позволяют вам задавать одну матрицу перед формирователем M, или две матрицы перед формирователем учитывают таким образом что M = M _1M2. Это дает возможность задавать предварительный формирователь, который вычисляется из факторизации, такой как M = L U. Обратите внимание на то, что в разделении предобусловленный случай, где M = H H ^T также содержит, нет прямой связи между M1 и M2 входные параметры с факторами H.

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

Функция Факторизация Описание

Функция	Факторизация	Описание
`ilu`	$A \approx Лютеций$	Неполная LU-факторизация для квадратных или прямоугольных матриц.
`ichol`	$A \approx L L^{*}$	Неполная факторизация Холесского для симметричных положительных определенных матриц.

ilu

$A \approx Лютеций$

Неполная LU-факторизация для квадратных или прямоугольных матриц.

ichol

$A \approx L L^{*}$

Неполная факторизация Холесского для симметричных положительных определенных матриц.

Смотрите Неполные Факторизации для получения дополнительной информации о ilu и ichol.

Пример перед формирователем

Рассмотрите приближение конечной разности с пятью точками к уравнению Лапласа на квадратной, двумерной области. Следующие команды используют предварительный формирователь метода предобусловленного метода сопряженных градиентов (PCG) M = L*L', где L заполнение нулями неполного Фактора Холесского A.

A = delsq(numgrid('S',250));
b = ones(size(A,1),1);
tol = 1e-3;
maxit = 100;
L = ichol(A);
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 92 to a solution with relative residual 0.00076.

pcg требует, чтобы 92 итерации достигли заданного допуска. Однако использование различного предварительного формирователя может привести к лучшим результатам. Например, использование ichol создать модифицированный неполный Холесский позволяет pcg соответствовать заданному допуску только после 39 итераций.

L = ichol(A,struct('type','nofill','michol','on'));
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 39 to a solution with relative residual 0.00098.

Уравновешивание и переупорядочение

Для в вычислительном отношении жестких проблем вам может быть нужен лучший предварительный формирователь, чем тот, сгенерированный ilu или ichol непосредственно. Например, вы можете хотеть сгенерировать лучший качественный предварительный формирователь или минимизировать объем сделанного расчета. В этих случаях можно использовать уравновешивание, чтобы сделать матрицу коэффициентов более по диагонали доминирующей (который может привести к лучшему качественному предварительному формирователю), и переупорядочивающий, чтобы минимизировать количество ненулей в матричных факторах (который может уменьшать требования к памяти и может повысить эффективность последующих вычислений).

Если вы используете и уравновешивание и переупорядочивающий, чтобы сгенерировать предварительный формирователь, процесс:

Используйте equilibrate на матрице коэффициентов.
Переупорядочьте уравновешенную матрицу с помощью функции переупорядочения разреженной матрицы, такой как dissect или symrcm.
Сгенерируйте итоговый предварительный формирователь с помощью ilu или ichol.

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

Создайте матрицу коэффициентов A и вектор из единиц b для правой стороны линейного уравнения. Вычислите оценку числа обусловленности для A.
```
load west0479
A = west0479;
n = size(A,1);
b = ones(n,1);
condest(A)
```
```
ans =

   1.4244e+12
```
Используйте equilibrate улучшить число обусловленности матрицы коэффициентов.
```
[P,R,C] = equilibrate(A);
Anew = R*P*A*C;
bnew = R*P*b;
condest(Anew)
```
```
ans =

   5.1042e+04
```
Переупорядочьте уравновешенную матрицу с помощью dissect.
```
q = dissect(Anew);
Anew = Anew(q,q);
bnew = bnew(q);
```
Сгенерируйте предварительный формирователь с помощью неполной LU-факторизации.
```
[L,U] = ilu(Anew);
```
Решите линейную систему с gmres с помощью матриц перед формирователем, допуска 1e-10, 50 максимальных внешних итераций и 30 внутренних итераций.
```
tol = 1e-10;
maxit = 50;
restart = 30;
[xnew, flag, relres] = gmres(Anew,bnew,restart,tol,maxit,L,U);
x(q) = xnew;
x = C*x(:);
```
Теперь сравните relres относительная невязка возвращена gmres (который включает предварительные формирователи) к относительной невязке без предварительных формирователей resnew и относительная невязка без уравновешивания res. Результаты показывают, что даже при том, что линейные системы являются всем эквивалентом, различные методы применяют различные веса к каждому элементу, и это может значительно влиять на значение невязки.
```
relres
resnew = norm(Anew*xnew - bnew) / norm(bnew)
res = norm(A*x - b) / norm(b)
```
```
relres =
   8.7537e-11
resnew =
   3.6805e-08
res =
   5.1415e-04
```

Используя линейные операторы вместо матриц

Итеративные решатели в MATLAB не требуют, чтобы вы обеспечили числовую матрицу для A. Поскольку вычисления, выполняемые решателями, в основном, используют результат матричного векторного умножения A*x или A'*x, можно вместо этого обеспечить функцию, которая вычисляет результат тех линейных операций. Функция, которая вычисляет эти количества, часто вызывается linear operator.

В дополнение к использованию линейного оператора вместо матрицы коэффициентов A, можно также использовать линейный оператор вместо матрицы для предварительного формирователя M. В этом случае функция должна вычислить M\x или M'\x.

Используя линейные операторы позволяет вам использовать шаблоны в A или M вычислить значение линейных операций более эффективно, чем если бы решатель выполнил полные матричные векторные умножения. Это также означает, что вам не нужна память, чтобы сохранить коэффициент или матрицы перед формирователем, поскольку линейный оператор обычно вычисляет результат матричного векторного умножения, не формируя матрицу вообще.

Например, рассмотрите матрицу коэффициентов

A = [2 -1  0  0  0  0;
    -1  2 -1  0  0  0;
     0 -1  2 -1  0  0;
     0  0 -1  2 -1  0;
     0  0  0 -1  2 -1;
     0  0  0  0 -1  2];

Когда эта матрица коэффициентов умножает вектор, единственные ненулевые элементы - те, которые умножают трехдиагональные элементы A. Так, для данного векторного x, функция линейного оператора просто должна добавить вместе три вектора, чтобы вычислить значение A*x:

function y = linearOperatorA(x)
y = -1*[0; x(1:end-1)] ...
  + 2*x ...
  + -1*[x(2:end); 0];
end

Большинство итеративных решателей требует функции линейного оператора для A возвращать значение A*x. Аналогично, для матрицы перед формирователем M, функция обычно должна вычислять M\x. Для решателей lsqrqmr, и bicg, функция линейного оператора должна также возвратить значение для A'*x или M'\x когда требуется. Смотрите итеративные страницы с описанием решателя для примеров и описаний функций линейного оператора.

Ссылки

[1] Барретт, R., М. Берри, Т. Ф. Чан, и др., Шаблоны для Решения Линейных систем: Базовые блоки для Итерационных методов, SIAM, Филадельфия, 1994.

[2] Саад, Yousef, итерационные методы для разреженных линейных уравнений. PWS Publishing Company, 1996.

Документация