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

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

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

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

x = qmr(A,b);

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

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

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

x = qmr(A,b,1e-4,100);

qmr 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.06.

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

Рассчитать неполную факторизацию Холеского Aи используйте L' фактор в качестве предварительного условия ввода в qmr.

L = ichol(A);
x = qmr(A,b,1e-4,100,L');

qmr converged at iteration 76 to a solution with relative residual 6.7e-05.

Использование предварительного условия улучшает численные свойства проблемы настолько, что qmr способен сходиться.

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

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

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

fl0 = 1

rr0

rr0 = 0.7984

it0

it0 = 17

fl0 равно 1, потому что qmr не сходится с требуемым допуском 1e-12 в пределах запрошенных 20 итераций. 17-й итерат является наилучшим приближенным решением и является возвращенным, как указано it0 = 17.

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

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

fl1 = 0

rr1

rr1 = 4.1557e-14

it1

it1 = 6

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

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

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

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

Создайте тридиагональную разреженную матрицу. Используйте сумму каждой строки в качестве вектора для правой стороны $\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);

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

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

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

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

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

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

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

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

Создайте несимметричную тридиагональную матрицу. Предварительный просмотр матрицы.

A = gallery('wilk',21) + diag(ones(20,1),1)

A = 21×21

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

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

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

$\mathit{A}\mathit{}\begin{array}{ccccccc}\mathrm{}& \mathrm{}& \mathrm{}& & & & \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & \\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& & \\ & \mathrm{}& \mathrm{}& \mathrm{}& & & \\ & & \mathrm{}& & \mathrm{}& & \mathrm{}\\ & & & & & & \mathrm{}\\ \mathrm{}& & & & \mathrm{}& \mathrm{}& \mathrm{x=[1020\cdots \cdots 01920⋮01\ddots 20⋮010\ddots \ddots ⋮0\ddots 1\ddots 0⋮\ddots \ddots \ddots 20\cdots \cdots 0110}\end{array}]\begin{array}{c}{\mathit{[}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\\ \\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\end{array}x1x2x3⋮⋮x21]=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21$].

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

$\mathit{A}\mathit{}\begin{array}{c}{\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\\ \\ \\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\end{array}$x=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21]=[0x1x2⋮x20]+[10x19x2⋮9x2010x21]+2⋅[x2x3⋮x210$\begin{array}{c}\mathrm{}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathrm{}\mathit{}}_{\mathrm{}}\\ {\mathrm{}\mathit{}}_{\mathrm{}}\\ \\ \mathrm{}{\mathit{}}_{\mathrm{}}\\ \mathrm{}{\mathit{}}_{\mathrm{}}\end{array}\mathrm{}\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\\ \mathrm{}\end{array}$].

Аналогично, выражение для ${\mathit{}}^{\mathit{AT}}\mathit{x}$ становится следующим:

${\mathit{}}^{\mathit{AT}}\mathit{}\begin{array}{ccccccc}\mathrm{}& \mathrm{}& \mathrm{}& & & & \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& & & \\ \mathrm{}& \mathrm{}& & \mathrm{}& \mathrm{}& & \\ & \mathrm{}& \mathrm{}& \mathrm{}& & & \\ & & \mathrm{}& & \mathrm{}& & \mathrm{}\\ & & & & & & \mathrm{}\\ \mathrm{}& & & & \mathrm{}& \mathrm{}& \mathrm{x=[1010\cdots \cdots 02910⋮02\ddots 10⋮020\ddots \ddots ⋮0\ddots 1\ddots 0⋮\ddots \ddots \ddots 10\cdots \cdots 0210}\end{array}]\begin{array}{c}{\mathit{[}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathrm{}\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ \\ {\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\end{array}x1x2x3⋮⋮x21]=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21$].

${\mathit{}}^{\mathit{AT}}\mathit{}\begin{array}{c}{\mathrm{}\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ \\ \\ {\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}{\mathit{}}_{\mathrm{}}\\ {\mathrm{}\mathit{}}_{\mathrm{}}\mathrm{}{\mathit{}}_{\mathrm{}}\end{array}\mathrm{}\begin{array}{c}\mathrm{}\\ {\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathrm{}\mathit{}}_{\mathrm{}}\\ {\mathrm{}\mathit{}}_{\mathrm{}}\\ \\ \mathrm{}{\mathit{}}_{\mathrm{}}\\ \mathrm{}{\mathit{}}_{\mathrm{}}\end{array}\begin{array}{c}{\mathit{}}_{\mathrm{}}\\ {\mathit{}}_{\mathrm{}}\\ \\ {\mathit{}}_{\mathrm{}}\\ \mathrm{x=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21]=2\cdot [0x1x2⋮x20]+[10x19x2⋮9x2010x21]+[x2x3⋮x210}\end{array}$].

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

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

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

Теперь решите линейную систему $\mathrm{Ax}\mathit{=}$ b, предоставивqmr с дескриптором функции, который вычисляет A*x и A'*x. Использовать допуск 1e-6 и 25 итераций. Укажите $\mathit{b}$ в качестве сумм строк $\mathit{A}$ так, чтобы истинное решение для $\mathit{x}$ было вектором единиц.

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = qmr(@afun,b,tol,maxit)

qmr converged at iteration 19 to a solution with relative residual 4.7e-07.

x1 = 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,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end