Решение квадратной линейной системы с помощью gmres с настройками по умолчанию, а затем скорректировать допуск и количество итераций, используемых в процессе решения.

Создание случайной разреженной матрицы A с 50% плотностью и ненулем на главной диагонали. Также создать случайный вектор b для правой стороны $\mathrm{Ax}\mathit{=}$b.

rng default
A = sprandn(400,400,0.5) + 12*speye(400);
b = rand(400,1);

Решить $\mathrm{Ax}\mathit{=}$ b с помощьюgmres. Выходной дисплей включает в себя значение относительной остаточной ошибки $\frac{\Vert \mathit{}\mathrm{}b-Ax}{\mathit{\Vert }}$‖ b ‖.

x = gmres(A,b);

gmres stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.35.

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

Повторное решение системы с использованием допуска 1e-4 и 100 итераций.

tol = 1e-4;
maxit = 100;
x = gmres(A,b,[],tol,maxit);

gmres stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.0045.

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

Повторное решение системы с помощью restart значение 100 и a maxit значение 20. Вместо того, чтобы делать 100 итераций один раз, gmres выполняет 100 итераций между перезапусками и повторяет их 20 раз.

restart = 100;
maxit = 20;
x = gmres(A,b,restart,tol,maxit);

gmres(100) converged at outer iteration 2 (inner iteration 75) to a solution with relative residual 9.3e-05.

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

Анализ влияния использования матрицы предварительного кондиционирования без перезапуска gmres для решения линейной системы.

Загрузка west0479, вещественная несимметричная разреженная матрица 479 на 479.

load west0479
A = west0479;

Определить b так что истинное решение $\mathrm{Ax}\mathit{=}$ b является вектором всех.

b = sum(A,2);

Задайте допуск и максимальное количество итераций.

tol = 1e-12;
maxit = 20;

Использовать gmres для поиска решения по требуемому допуску и количеству итераций. Укажите пять выходов для возврата информации о процессе решения:

[x,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0

fl0 = 1

rr0

rr0 = 0.7603

it0

it0 = 1×2

     1    20

fl0 равно 1, потому что gmres не сходится с требуемым допуском 1e-12 в пределах запрошенных 20 итераций. Лучшее приблизительное решение, которое gmres возвращает последнее (как указано it0(2) = 20). MATLAB сохраняет остаточную историю в rv0.

Чтобы помочь с медленной сходимостью, можно указать матрицу предварительного условия. С тех пор A несимметричен, используется ilu для создания устройства предварительного кондиционирования $\mathit{M}\mathit{=}\mathit{L}$U. Задание допуска перетаскивания для игнорирования недиагональных записей со значениями меньше 1e-6. Определите предварительно кондиционированную систему ${\mathit{}}^{\mathrm{}}\mathit{M-1A}\mathit{x}{\mathit{=}}^{\mathrm{}}\mathit{}$M-1b, указавL и U в качестве входных данных для gmres. Обратите внимание, что при указании предварительного условия gmres вычисляет остаточную норму системы предварительного кондиционирования для выходных данных rr1 и rv1.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1

fl1 = 0

rr1

rr1 = 7.2303e-14

it1

it1 = 1×2

     1     6

Использование ilu средство предварительного кондиционирования производит относительный остаток, меньший, чем предписанный допуск 1e-12 на шестой итерации. Первый остаток rv1(1) является norm(U\(L\b)), где M = L*U. Последний остаток rv1(end) является norm(U\(L\(b-A*x1))).

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

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Использование предварительного условия с перезапущенным gmres.

Загрузка west0479, вещественная несимметричная разреженная матрица 479 на 479.

load west0479
A = west0479;

Определить b так что истинное решение $\mathrm{Ax}\mathit{=}$ b является вектором всех.

b = sum(A,2);

Постройте неполное предварительное кондиционирующее устройство логической единицы с допуском сброса 1e-6.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

Преимущество перезапуска gmres означает ограничение объема памяти, необходимого для выполнения метода. Без перезапуска, gmres требует maxit векторы хранения для сохранения основы крыловского подпространства. Также, gmres должен ортогонализироваться по всем предыдущим векторам на каждом шаге. Перезапуск ограничивает объем используемой рабочей области и объем работы на одну внешнюю итерацию.

Выполнить gmres(3), gmres(4), и gmres(5) использование неполных коэффициентов логической единицы в качестве предварительных условий. Использовать допуск 1e-12 и максимум 20 внешних итераций.

tol = 1e-12;
maxit = 20;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,3,tol,maxit,L,U);
[x4,fl4,rr4,it4,rv4] = gmres(A,b,4,tol,maxit,L,U);
[x5,fl5,rr5,it5,rv5] = gmres(A,b,5,tol,maxit,L,U);
fl3

fl3 = 0

fl4

fl4 = 0

fl5

fl5 = 0

fl3, fl4, и fl5 равны 0, потому что в каждом случае перезапускаются gmres приводит относительный остаток к меньшему, чем предписанный допуск 1e-12.

На следующем графике показана история сходимости каждого перезапускаемого gmres способ. gmres(3) сходится на внешней итерации 5, внутренней итерации 3 (it3 = [5, 3]), которая будет такой же, как внешняя итерация 6, внутренняя итерация 0, следовательно, маркировка 6 на конечной отметке.

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

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

Проверить влияние подачи gmres с первоначальным предположением решения.

Создайте тридиагональную разреженную матрицу. Используйте сумму каждой строки в качестве вектора для правой стороны $\mathrm{Ax}\mathit{=}$ b, чтобы ожидаемое решение $\mathit{для}$ x было вектором единиц.

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

Использовать gmres для решения $\mathrm{Ax}\mathit{=}$ b дважды: один раз с начальным предположением по умолчанию и один раз с хорошим начальным предположением решения. Используйте 200 итераций и допуск по умолчанию для обоих решений. Укажите начальное предположение во втором решении как вектор со всеми элементами, равными 0.99.

maxit = 200;
x1 = gmres(A,b,[],[],maxit);

gmres converged at iteration 27 to a solution with relative residual 9.5e-07.

x0 = 0.99*e;
x2 = gmres(A,b,[],[],maxit,[],[],x0);

gmres converged at iteration 7 to a solution with relative residual 6.7e-07.

В этом случае предоставление начального предположения позволяет gmres для более быстрого сближения.

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = gmres(A,b,[],tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

Решение линейной системы путем предоставления gmres с дескриптором функции, который вычисляет A*x вместо матрицы коэффициентов A.

Одна из тестовых матриц Уилкинсона, сгенерированных gallery представляет собой тридиагональную матрицу 21 на 21. Предварительный просмотр матрицы.

A = gallery('wilk',21)

A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

Матрица Уилкинсона имеет специальную структуру, поэтому можно представить операцию A*x с дескриптором функции. Когда A умножает вектор, большинство элементов в полученном векторе равны нулю. Ненулевые элементы в результате соответствуют ненулевым тридиагональным элементам A. При этом только основная диагональ имеет ненулевые значения, не равные 1.

Выражение $\mathrm{Ax}$ становится следующим:

$\mathrm{}\begin{array}{cccccccccc}\mathrm{}& \mathrm{}& \mathrm{}& & & & & & \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & & & \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & & \\ & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & & \\ & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & \\ & & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & \\ & & & & \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ & & & & & \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}\\ \mathrm{}& & & & & & \mathrm{}& & & \mathrm{}\\ \mathrm{}& \mathrm{}& & & & & & \mathrm{}& \mathrm{}& \mathrm{Ax=[1010\cdots \cdots \cdots 001910001810⋮⋮0171001610⋮⋮0151001410⋮⋮013\ddots 000\ddots \ddots 100\cdots \cdots \cdots 0110}\end{array}]\begin{array}{c}{\mathit{[}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ \\ \\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\end{array}x1x2x3x4x5⋮⋮x21]=[10x1+x2x1+9x2+x3x2+8x3+x4⋮x19+9x20+x21x20+10x21$].

Результирующий вектор может быть записан как сумма трех векторов:

$\mathrm{}\begin{array}{c}\mathrm{}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathrm{}\end{array}$Ax=[0+10x1+x2x1+9x2+x3x2+8x3+x4⋮x19+9x20+x21x20+10x21+0]=[0x1⋮x20]+[10x19x2⋮10x21]+[x2⋮x210$\begin{array}{c}\mathrm{}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}\mathrm{}{\mathit{}}_{\mathrm{}}\\ \mathrm{}{\mathit{}}_{\mathrm{}}\\ \\ \mathrm{}{\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\\ \mathrm{}\end{array}$].

В MATLAB ® запишите функцию, которая создает эти векторы и добавляет их вместе, давая таким образом значение A*x:

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

(Эта функция сохраняется как локальная функция в конце примера.)

Теперь решите линейную систему $\mathrm{Ax}\mathit{=}$ b, предоставивgmres с дескриптором функции, который вычисляет A*x. Использовать допуск 1e-12, 15 внешних итераций и 10 внутренних итераций перед перезапуском.

b = ones(21,1);
tol = 1e-12;  
maxit = 15;
restart = 10;
x1 = gmres(@afun,b,restart,tol,maxit)

gmres(10) converged at outer iteration 5 (inner iteration 10) to a solution with relative residual 5.3e-13.

x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

Проверьте, что afun(x1) производит вектор из единиц.

afun(x1)

ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end