22.3 Итеративные методы, применяемые к разреженным матрицам

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

: x = pcg (A, b, tol, maxit, m1, m2, x0, …) ¶

: x = pcg (A, b, tol, maxit, M, [], x0, …) ¶

: [x, flag, relres, iter, resvec, eigest] = pcg (A, b, …) ¶

Решает линейную систему уравнений A * x = b с помощью итерационного метода предварительно обусловленного сопряженного градиента.

Входные аргументы:

• A — матрица линейной системы, должна быть квадратной. A может быть передана как матрица, дескриптор функции или встроенная функция Afun таким образом, чтобы Afun(x) = A * x. Дополнительные параметры для Afun могут быть переданы после x0.

  A должна быть эрмитовой и положительно определённой (ЭПД). Если pcg обнаруживает, что A не является положительно определённой, выводится предупреждение, и выходной параметр flag устанавливается.

• b — вектор правой части.
• tol — требуемая относительная погрешность остатка, b - A * x. Итерация прекращается, если norm (b - A * x) ≤ tol * norm (b). Если tol опущено или пусто, используется значение 1e-6.
• maxit — максимальное количество итераций; если maxit опущено или пусто, используется значение 20.
• m — матрица предварительного обуславливания ЭПД. Для любого разложения m = p1 * p2 такое, что inv (p1) * A * inv (p2) является ЭПД, метод сопряжённых градиентов формально применяется к линейной системе inv (p1) * A * inv (p2) * y = inv (p1) * b, с x = inv (p2) * y (раздельное предварительное обуславливание). На практике на каждой итерации метода сопряжённых градиентов решается линейная система с матрицей m с помощью mldivide. Если доступно конкретное разложение m = m1 * m2 (например, неполное разложение Холецкого для a), можно передать две матрицы m1 и m2, и соответствующие линейные системы решаются с помощью оператора mldivide. Обратите внимание, что правильный выбор предварительного обуславливателя может значительно улучшить общую производительность метода. Вместо матриц m1 и m2 пользователь может передать две функции, которые возвращают результаты применения обратной m1 и m2 к вектору. Если m1 опущено или пусто [], предварительное обуславливание не применяется. Если разложение m недоступно, m2 можно опустить или задать [], а переменная m1 может использоваться для передачи предварительного обуславливателя m.
• x0 — начальное приближение. Если x0 опущено или пусто, функция по умолчанию устанавливает x0 в нулевой вектор.

Аргументы, следующие за x0, обрабатываются как параметры и передаются соответствующим образом в любую из функций (A или m1 или m2), которые были переданы pcg. Дополнительные сведения см. в примерах ниже.

Выходные аргументы:

• x — вычисленное приближение решения A * x = b. Если алгоритм не сошёлся, то x — итерация с минимальным остатком.
• flag сообщает о сходимости:
  • 0: Алгоритм сошёлся в пределах заданной погрешности.
  • 1: Алгоритм не сошёлся и достиг максимального количества итераций.
  • 2: Матрица предварительного обуславливания вырожденная.
  • 3: Алгоритм застыл, т. е., абсолютная величина разницы между текущей итерацией x и предыдущей меньше eps * norm (x,2).
  • 4: Алгоритм обнаружил, что входная (предварительно обусловленная) матрица не является ЭПД.
• relres — отношение конечного остатка к начальному значению, измеренное в евклидовой норме.
• iter указывает итерацию x, при которой он был вычислен. Поскольку выходной x соответствует решению с минимальным остатком, общее количество итераций, проведённых методом, определяется по length(resvec) - 1.
• resvec описывает историю сходимости метода. resvec (i, 1) — евклидова норма остатка, а resvec (i, 2) — норма предварительно обусловленного остатка после (i-1)-ой итерации, i = 1, 2, …, iter+1. Норма предварительно обусловленного остатка определяется как r' * (m \ r) , где r = b - A * x, см. также описание m. Если eigest не требуется, возвращается только resvec (:, 1).
• eigest возвращает оценку наименьшего eigest(1) и наибольшего eigest(2) собственных значений предварительно обусловленной матрицы P = m \ A. В частности, если предварительное обуславливание не используется, возвращаются оценки крайних собственных значений A. eigest(1) является завышенной оценкой, а eigest(2) — заниженной оценкой, так что eigest(2) / eigest(1) является нижней границей для cond (P, 2), которая, тем не менее, в пределе теоретически должна быть равна фактическому значению числа обусловленности.

Рассмотрим тривиальную задачу с треугольной матрицей

n = 10;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
b = A * ones (n, 1);
M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
M2 = M1';
M = M1 * M2;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

ПРИМЕР 1: Простейшее использование pcg

x = pcg (A, b)

ПРИМЕР 2: pcg с функцией, вычисляющей A * x

x = pcg (Afun, b)

ПРИМЕР 3: pcg с матрицей предварительного обуславливания M

x = pcg (A, b, 1e-06, 100, M)

ПРИМЕР 4: pcg с функцией в качестве предварительного обуславливателя

x = pcg (Afun, b, 1e-6, 100, Mfun)

ПРИМЕР 5: pcg с матрицами предварительного обуславливания M1 и M2

x = pcg (A, b, 1e-6, 100, M1, M2)

ПРИМЕР 6: pcg с функциями в качестве предварительных обуславливателей

x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun)

ПРИМЕР 7: pcg с входом в функцию, требующую аргумент

function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfun = @(x, p) Ap (A, x, p);
x = pcg (Apfun, b, [], [], [], [], [], 2);

ПРИМЕР 8: явный пример, показывающий, что pcg использует раздельное предварительное обуславливание

M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
M2 = M1';
M = M1 * M2;

## reference solution computed by pcg after two iterations
[x_ref, fl] = pcg (A, b, [], 2, M)

## split preconditioning
[y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
x = M2 \ y # compare x and x_ref

Ссылки:

1. C.T. Kelley, Итерационные методы для линейных и нелинейных уравнений, SIAM, 1995. (основной алгоритм PCG)
2. Y. Saad, Итерационные методы для разреженных линейных систем, PWS 1996. (оценка числа обусловленности из PCG) Пересмотренный вариант этой книги доступен онлайн по адресу https://www-users.cs.umn.edu/~saad/books.html

См. также: sparse, pcr, gmres, bicg, bicgstab, cgs.

: x = pcr (A, b, tol, maxit, m, x0, …) ¶

: [x, flag, relres, iter, resvec] = pcr (…) ¶

Решить систему линейных уравнений A * x = b с помощью итерационного метода Предварительно обусловленных сопряженных градиентов.

Входные аргументы:

• A может быть квадратной (желательно разреженной) матрицей или функцией-обработчиком, встроенной функцией или строкой, содержащей имя функции, которая вычисляет A * x. В принципе A должна быть симметричной и невырожденной; если pcr обнаружит, что A является вырожденной, вы получите сообщение об ошибке, и параметр вывода flag будет установлен.
• b — вектор правой части.
• tol — требуемая относительная погрешность остатка, b - A * x. Итерация останавливается, если norm (b - A * x) <= tol * norm (b - A * x0). Если tol пуста или опущена, функция устанавливает значение tol = 1e-6 по умолчанию.
• maxit — максимальное допустимое количество итераций; если [] задано для maxit, или pcr содержит меньше аргументов, используется значение по умолчанию, равное 20.
• m — матрица (левого) предварительного обусловливания, так что итерация (теоретически) эквивалентна решению с помощью pcr P * x = m \ b, при P = m \ A. Обратите внимание, что правильный выбор предварительного обусловливания может значительно улучшить общую производительность метода. Вместо матрицы m пользователь может передать функцию, возвращающую результаты применения обратной матрицы m к вектору (обычно это предпочтительный способ использования предварительного обусловливания). Если [] задано для m, или m опущено, предварительное обусловливание не применяется.
• x0 — начальное приближение. Если x0 пуста или опущена, функция устанавливает x0 в нулевой вектор по умолчанию.

Аргументы, которые следуют за x0, обрабатываются как параметры и передаются должным образом в любые функции (A или m), которые передаются в pcr. См. примеры ниже для получения дополнительной информации.

Выходные аргументы:

• x — вычисленное приближение решения A * x = b.
• flag сообщает о сходимости. flag = 0 означает, что решение сошлось, и критерий точности, заданный tol, выполнен. flag = 1 означает, что предел maxit для числа итераций был достигнут. flag = 3 сообщает об ошибке в расчете, см. [1] для получения подробной информации.
• relres — отношение конечного остатка к его начальному значению, измеренное в евклидовой норме.
• iter — фактическое количество выполненных итераций.
• resvec описывает историю сходимости метода, так что resvec (i) содержит евклидовы нормы остатка после (i-1)-й итерации, i = 1,2, …, iter+1.

Рассмотрим тривиальную задачу с диагональной матрицей (мы используем разреженность A)

n = 10;
A = sparse (diag (1:n));
b = rand (N, 1);

ПРИМЕР 1: Простейшее использование pcr

x = pcr (A, b)

ПРИМЕР 2: pcr с функцией, которая вычисляет A * x.

function y = apply_a (x)
  y = [1:10]' .* x;
endfunction

x = pcr ("apply_a", b)

ПРИМЕР 3: Предварительно обусловленная итерация с полным диагностированием. Предварительное обусловливание (довольно странное, потому что даже исходная матрица A тривиальна) определяется как функция

function y = apply_m (x)
  k = floor (length (x) - 2);
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], "apply_m")
semilogy ([1:iter+1], resvec);

ПРИМЕР 4: Наконец, предварительное обусловливание, зависящее от параметра k.

function y = apply_m (x, varargin)
  k = varargin{1};
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], "apply_m"', [], 3)

Ссылки:

В. Хакубуш, Итеративное решение больших разреженных систем уравнений, раздел 9.5.4; Springer, 1994

См. также: sparse, pcg.

Скорость сходимости итерационного решателя к решению можно ускорить с помощью матрицы предварительного обусловливания M. В этом случае вместо этого решается линейное уравнение M^-1 * x = M^-1 * A \ b. Типичные матрицы предварительного обусловливания являются частичными факторизациями исходной матрицы.

: L = ichol (A) ¶

: L = ichol (A, opts) ¶

Вычислить неполную факторизацию Холецкого разреженной квадратной матрицы A.

По умолчанию ichol использует только нижнюю треугольную часть A и создает нижнюю треугольную факторизацию L такую, что L*L' приближает A.

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

Факторизацию можно изменить, передав параметры в структуре opts. Имя параметра — поле структуры, а значение — значение поля. Названия и спецификаторы чувствительны к регистру.

type

Тип факторизации.

"nofill" (по умолчанию)

Неполная факторизация Холецкого без добавления элементов (IC(0)).

"ict"

Неполная факторизация Холецкого с отбрасыванием элементов по порогу (ICT).

diagcomp

Неотрицательное скалярное значение alpha для неполной факторизации Холецкого для A + alpha * diag (diag (A)) вместо A. Это может быть полезно, когда A не является положительно определенной. Значение по умолчанию равно 0.

droptol

Неотрицательное скалярное значение, задающее порог отбрасывания для факторизации, если выполняется ICT. Значение по умолчанию равно 0, что дает полную факторизацию Холецкого.

Недиагональные элементы L устанавливаются в 0, если

abs (L(i,j)) >= droptol * norm (A(j:end, j), 1).

michol

Модифицированная неполная факторизация Холецкого:

"off" (по умолчанию)

Суммы по строкам и столбцам не обязательно сохраняются.

"on"

Диагональ L модифицируется так, что суммы по строкам (и столбцам) сохраняются даже при отбрасывании элементов во время факторизации. Сохраненное соотношение: A * e = L * L' * e, где e — вектор единиц.

shape

"lower" (по умолчанию)

Используйте только нижнюю треугольную часть A и возвращайте нижнюю треугольную факторизацию L такую, что L*L' приближает A.

"upper"

Используйте только верхнюю треугольную часть A и возвращайте верхнюю треугольную факторизацию U такую, что U'*U приближает A.

ПРИМЕРЫ

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

A = [ 0.37, -0.05,  -0.05,  -0.07;
     -0.05,  0.116,  0.0,   -0.05;
     -0.05,  0.0,    0.116, -0.05;
     -0.07, -0.05,  -0.05,   0.202];
A = sparse (A);
nnz (tril (A))
ans =  9
L = chol (A, "lower");
nnz (L)
ans =  10
norm (A - L * L', "fro") / norm (A, "fro")
ans =  1.1993e-16
opts.type = "nofill";
L = ichol (A, opts);
nnz (L)
ans =  9
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.019736

Еще один пример факторизации — матрица конечных разностей, используемая для решения задачи на собственные значения на единичном квадрате.

nx = 400; ny = 200;
hx = 1 / (nx + 1); hy = 1 / (ny + 1);
Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
               [-1 0 1 ], nx, nx) / (hx ^ 2);
Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
               [-1 0 1 ], ny, ny) / (hy ^ 2);
A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
nnz (tril (A))
ans =  239400
opts.type = "nofill";
L = ichol (A, opts);
nnz (tril (A))
ans =  239400
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.062327

Ссылки на реализованные алгоритмы:

[1] Ю. Саад. "Методы предварительного обусловливания." Итерационные методы для разреженных систем линейных уравнений, PWS Publishing Company, 1996.

[2] М. Джонс, П. Плассманн: Улучшенная неполная факторизация Холецкого, 1992.

См. также: chol, ilu, pcg.

: ilu (A) ¶

: ilu (A, opts) ¶

: [L, U] = ilu (…) ¶

: [L, U, P] = ilu (…) ¶

Вычислить неполное LU-разложение разреженной квадратной матрицы A.

ilu возвращает единичную нижнюю треугольную матрицу L, верхнюю треугольную матрицу U и необязательно матрицу перестановок P, такую что L*U приближает P*A.

Факторы, полученные с помощью этой функции, могут быть полезны в качестве предварительных условий для системы линейных уравнений, решаемых итерационными методами, такими как BICG (Бисопряжённые градиенты) или GMRES (Обобщённый метод наименьших остатков).

Разложение можно модифицировать, передавая опции в структуре opts. Имя опции — поле структуры, а значение — значение поля. Имена и спецификаторы чувствительны к регистру.

type

Тип разложения.

"nofill" (по умолчанию)

ILU-разложение без заполнения (ILU(0)).

Дополнительные поддерживаемые опции: milu.

"crout"

Версия разложения ILU по методу Кроута (ILUC).

Дополнительные поддерживаемые опции: milu, droptol.

"ilutp"

ILU-разложение с порогом и выбором опорных элементов.

Дополнительные поддерживаемые опции: milu, droptol, udiag, thresh.

droptol

Неотрицательное скалярное значение, задающее порог отбрасывания для разложения. Значение по умолчанию равно 0, что приводит к полному LU-разложению.

Недиагональные элементы U устанавливаются в 0, если

abs (U(i,j)) >= droptol * norm (A(:,j)).

Недиагональные элементы L устанавливаются в 0, если

abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j).

milu

Модифицированное неполное LU-разложение:

"row"

Модифицированное неполное LU-разложение по строкам. Разложение сохраняет суммы по строкам: A * e = L * U * e, где e — вектор из единиц.

"col"

Модифицированное неполное LU-разложение по столбцам. Разложение сохраняет суммы по столбцам: e' * A = e' * L * U.

"off" (по умолчанию)

Суммы по строкам и столбцам не обязательно сохраняются.

udiag

Если значение true, все нули на диагонали верхней треугольной матрицы заменяются локальным порогом отбрасывания droptol * norm (A(:,j))/U(j,j). Значение по умолчанию — false.

thresh

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

Если ilu вызывается с одним выходным аргументом, возвращаемая матрица — L + U - speye (size (A)), где L — единичная нижняя треугольная матрица, а U — верхняя треугольная матрица.

С двумя выходными аргументами ilu возвращает единичную нижнюю треугольную матрицу L и верхнюю треугольную матрицу U. Для opts.type == "ilutp", один из факторов переставляется в соответствии со значением opts.milu. Когда opts.milu == "row", U является переставленной по столбцам верхней треугольной матрицей. В противном случае, L является переставленной по строкам единичной нижней треугольной матрицей.

Если указано три выходных аргумента и opts.milu != "row", возвращается P, такая что L и U — неполные факторы от P*A. Когда opts.milu == "row", возвращается P, такая что L и U — неполные факторы от A*P.

ПРИМЕРЫ

A = gallery ("neumann", 1600) + speye (1600);
opts.type = "nofill";
nnz (A)
ans = 7840

nnz (lu (A))
ans = 126478

nnz (ilu (A, opts))
ans = 7840

Это демонстрирует, что A содержит 7840 ненулевых элементов, полное LU-разложение — 126478, а неполное LU-разложение с 0 уровнем заполнения — 7840, что совпадает с количеством ненулевых элементов в A. Взято с: https://www.mathworks.com/help/matlab/ref/ilu.html

A = gallery ("wathen", 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = "crout";
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)

Этот пример использует ILU в качестве предварительного условия для случайной матрицы конечных элементов, которая имеет большое число обусловленности. Без L и U BICG не сойдётся.

См. также: lu, ichol, bicg, gmres.