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

Одним из наиболее важных и распространенных приложений числовой линейной алгебры является решение линейных систем, которое может быть выражено в форме 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 реализует прямые методы через матричные операторы деления / и \, а также такие функции, как decomposition, lsqminnorm, и linsolve.
Итерационные методы дают приблизительное решение линейной системы после конечного числа шагов. Эти методы полезны для больших систем уравнений, где разумно сравнить точность для более короткого времени запуска. Итерационные методы используют матрицу коэффициентов только косвенно, через матрично-векторный продукт или абстрактный линейный оператор. Итерационные методы могут использоваться с любой матрицей, но они обычно применяются к большим разреженным матрицам, для которых прямые решения являются медленными. Скорость решения линейной системы с косвенным методом не так сильно зависит от шаблона заполнения матрицы коэффициентов, как прямой метод. Однако использование итерационного метода обычно требует настройки параметров для каждой конкретной задачи.
Для примера этот код решает большую разреженную линейную систему, которая имеет симметричную положительно определенную матрицу коэффициентов.
```
A = delsq(numgrid('L',400));
b = ones(size(A,1),1);
x = pcg(A,b,[],1000);
norm(b-A*x)
```
```
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.
Применить A*x0 и обновляйте величину и направление вектора x0 на основе значения остаточных и других вычисленных количеств. Это шаг, на котором выполняется большинство расчетов.
Повторите шаги с 2 по 4 до тех пор, пока значение x0 достаточно хорошо, чтобы удовлетворить допуск.

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

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

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

Описание	Примечания
`pcg` (предварительно обусловленные сопряженные градиенты)	Матрица коэффициентов должна быть симметричной положительно определенной. Наиболее эффективный решатель для симметричных положительно определенных систем, поскольку требуется хранение только для ограниченного числа векторов.
`lsqr` (наименьшие квадраты)	Единственный решатель, доступный для прямоугольных систем. Аналитически эквивалентен методу сопряженных градиентов (PCG), применяемому к нормальным уравнениям `(A'A)x = A'*b`.
`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 почти для всех квадратных, несимметричных задач. Существуют некоторые случаи, когда алгоритмы двояких градиентов (bicg, bicgstab, cgs, и так далее) более эффективны, чем gmres, но их непредсказуемое поведение сходимости часто делает gmres лучший первоначальный выбор.

Workflow to choose an appropriate iterative solver to use for a given problem.

Предварительные кондиционеры

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

$A x = b$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Итеративные решатели в MATLAB позволяют вам задать одну матричную M прекондиционера или два матричных фактора прекондиционера, такие что M = M 1 M 2. Это облегчает определение предкондиционера в его факторизованной форме, такой как 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. Для этой системы, pcg не может найти решение, не задав матрицу предварительной подготовки.

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;
b = ones(size(A,1),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 использование матриц preconditioner, допуск 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 умножает вектор, большинство элементов получившегося вектора являются нулями. Ненулевые элементы в результате соответствуют ненулевым триидиагональным элементам 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. Для решателей lsqr, qmr, и bicgфункции линейного оператора должны также вернуть значение для A'*x или M'\x по запросу. См. итерационные страницы с описанием решателя для примеров и описаний линейных функций оператора.

Ссылки

[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Линейные Системы: Базовые блоки for Итерационные Методы, SIAM, Philadelphia, 1994.

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

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