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

x = lsqr(A,b) пытается решить систему линейных уравнений A*x = b для x использование метода наименьших квадратов. lsqr находит решение методом наименьших квадратов для x который минимизирует norm(b-A*x). Когда A является допустимым, решение методом наименьших квадратов также является решением линейной системы. Когда попытка успешна, lsqr отображение сообщения для подтверждения сходимости. Если lsqr не сходится после максимального количества итераций или остановок по какой-либо причине, оно отображает диагностическое сообщение, которое включает относительную невязку norm(b-A*x)/norm(b) и номер итерации, на котором остановился метод.

x = lsqr(A,b,tol) задает допуск для метода. Допуск по умолчанию 1e-6.

x = lsqr(A,b,tol,maxit) задает максимальное количество итераций для использования. lsqr отображает диагностическое сообщение, если оно не сходится в maxit итераций.

x = lsqr(A,b,tol,maxit,M) задает матрицу предварительной подготовки M и вычисляет x путем эффективного решения системы $A M^{- 1} y = b$ для y, где $y = M x$ . Использование матрицы preconditioner может улучшить числовые свойства задачи и эффективность вычисления.

x = lsqr(A,b,tol,maxit,M1,M2) задает факторы матрицы предкондиционера M таким образом M = M1*M2.

x = lsqr(A,b,tol,maxit,M1,M2,x0) задает начальное предположение для вектора решения x. По умолчанию это нулевой вектор.

[x,flag] = lsqr(___) возвращает флаг, который определяет, успешно ли сходился алгоритм. Когда flag = 0конвергенция прошла успешно. Можно использовать этот синтаксис выхода с любой из предыдущих комбинаций входных аргументов. Когда вы задаете flag выход, lsqr не отображает никаких диагностических сообщений.

[x,flag,relres] = lsqr(___) также возвращает остаточную ошибку вычисляемого решения x. Если flag является 0, затем x является решением методом наименьших квадратов, которое минимизирует norm(b-A*x). Если relres является маленьким, тогда x также является последовательным решением, поскольку relres представляет norm(b-A*x)/norm(b).

[x,flag,relres,iter] = lsqr(___) также возвращает число итерации iter при котором x был вычислен.

[x,flag,relres,iter,resvec] = lsqr(___) также возвращает вектор норм невязки при каждой итерации, включая первый остаточный norm(b-A*x0).

[x,flag,relres,iter,resvec,lsvec] = lsqr(___) также возвращается lsvec, которая является оценкой масштабированной ошибки нормального уравнения при каждой итерации.

Примеры

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

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

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

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

Решить $Ax = b$ использование lsqr. Выход отображения включает в себя значение относительной остаточной ошибки $\frac{‖ b - Ax ‖}{‖ b ‖}$ .

x = lsqr(A,b);

lsqr 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.26.

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

Снова решить систему с помощью допуска 1e-4 и 70 итераций. Задайте шесть выходов, чтобы вернуть относительную невязку relres вычисленного решения, а также истории остатков resvec и историю невязок методом наименьших квадратов lsvec.

[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,1e-4,70);
flag

flag = 0

Начиная с flag равным 0, алгоритм смог удовлетворить желаемому допуску ошибок при заданном количестве итераций. Обычно можно настроить допуск и количество итераций вместе, чтобы таким образом добиться компромиссов между скоростью и точностью.

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

relres

relres = 0.2625

lsres = lsvec(end)

lsres = 2.7640e-04

Эти нормы невязки указывают, что x является решением методом наименьших квадратов, потому что relres не меньше заданного допуска 1e-4. Поскольку нет последовательного решения для линейной системы, лучший, что может сделать решатель, это сделать наименьшие квадраты невязки удовлетворить допуск.

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

N = length(resvec);
semilogy(0:N-1,lsvec,'--o',0:N-1,resvec,'-o')
legend("Least-squares residual","Relative residual")

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Least-squares residual, Relative residual.

Использование `lsqr` с предварительным кондиционером

Исследуйте эффект использования матрицы предкондиционера с lsqr для решения линейной системы.

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

load west0479
A = west0479;

Определите b чтобы истинное решение $Ax = b$ является вектором всех таковых.

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

x является вычисляемым решением для A*x = b.
fl является флагом, указывающим, сходился ли алгоритм.
rr - относительная невязка вычисляемого ответа x.
it - число итерации, когда x был вычислен.
rv является вектором остаточной истории для $‖ b - Ax ‖$ .
lsrv является вектором остаточной истории методом наименьших квадратов.

[x,fl,rr,it,rv,lsrv] = lsqr(A,b,tol,maxit);
fl

fl = 1

rr

rr = 0.0017

it

it = 20

Начиная с fl = 1алгоритм не сходился к заданному допуску в пределах максимального количества итераций.

Для помощи с медленной сходимостью можно задать матрицу предкондиционера. Начиная с A несимметричен, используйте ilu чтобы сгенерировать предварительный кондиционер $M = L U$ в факторизованной форме. Задайте допуск удаления, чтобы игнорировать недиагональные значения со значениями меньше 1e-6. Решение предварительно обусловленной системы ${AM}^{- 1} (M x) = b$ для $y = Mx$ путем определения L и U как M1 и M2 входы в lsqr.

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

fl1 = 0

rr1

rr1 = 7.0954e-14

it1

it1 = 13

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

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

semilogy(0:length(rv)-1,rv/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')

Figure contains an axes. The axes contains 3 objects of type line, constantline. These objects represent No preconditioner, ILU preconditioner, Tolerance.

Исходное предположение о поставке

Исследуйте эффект подачи lsqr с начальным угадыванием решения.

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

A = sprand(700,900,0.1);
b = sum(A,2);

Использование lsqr решить $Ax = b$ дважды: один раз с начальным предположением по умолчанию и один раз с хорошим начальным предположением решения. Используйте 75 итераций и допуск по умолчанию для обоих решений. Задайте начальное предположение во втором решении как вектор со всеми элементами, равными 0,99.

maxit = 75;
x1 = lsqr(A,b,[],maxit);

lsqr converged at iteration 64 to a solution with relative residual 8.7e-07.

x0 = 0.99*ones(size(A,2),1);
x2 = lsqr(A,b,[],maxit,[],[],x0);

lsqr converged at iteration 26 to a solution with relative residual 9.6e-07.

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

Возврат промежуточных результатов

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

Для примера этот код выполняет 100 итераций четыре раза и сохраняет вектор решения после каждого прохода в цикле for-loop:

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

X(:,k) вектор решения вычисляется при итерации k системы for-loop и R(k) - относительная невязка этого решения.

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

Решить линейную систему путем предоставления lsqr с указателем на функцию, который вычисляет 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.

Выражение $A x$ становится:

$A x = [\begin{array}{c} 10 & 2 & 0 & \dots & \dots & 0 \\ 1 & 9 & 2 & 0 & ⋮ \\ 0 & 1 & ⋱ & 2 & 0 \\ ⋮ & 0 & 1 & 0 & ⋱ & ⋱ & ⋮ \\ 0 & ⋱ & 1 & ⋱ & 0 \\ ⋮ & ⋱ & ⋱ & ⋱ & 2 \\ 0 & \dots & \dots & 0 & 1 & 10 \end{array}] [\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ ⋮ \\ x_{21} \end{array}] = [\begin{array}{c} {10 x}_{1} + 2 x_{2} \\ x_{1} + 9 x_{2} + 2 x_{3} \\ ⋮ \\ ⋮ \\ x_{19} + 9 x_{20} + 2 x_{21} \\ x_{20} + 10 x_{21} \end{array}]$ .

Получившийся вектор может быть записан как сумма трех векторов:

$A x = [\begin{array}{c} {10 x}_{1} + 2 x_{2} \\ x_{1} + 9 x_{2} + 2 x_{3} \\ ⋮ \\ ⋮ \\ x_{19} + 9 x_{20} + 2 x_{21} \\ x_{20} + 10 x_{21} \end{array}]$ = $[\begin{array}{c} 0 \\ x_{1} \\ x_{2} \\ ⋮ \\ x_{20} \end{array}] + [\begin{array}{c} {10 x}_{1} \\ {9 x}_{2} \\ ⋮ \\ 9 x_{20} \\ 10 x_{21} \end{array}] + 2 \cdot [\begin{array}{c} x_{2} \\ x_{3} \\ ⋮ \\ x_{21} \\ 0 \end{array}]$ .

Аналогично, выражение для $A^{T} x$ становится:

$A^{T} x = [\begin{array}{c} 10 & 1 & 0 & \dots & \dots & 0 \\ 2 & 9 & 1 & 0 & ⋮ \\ 0 & 2 & ⋱ & 1 & 0 \\ ⋮ & 0 & 2 & 0 & ⋱ & ⋱ & ⋮ \\ 0 & ⋱ & 1 & ⋱ & 0 \\ ⋮ & ⋱ & ⋱ & ⋱ & 1 \\ 0 & \dots & \dots & 0 & 2 & 10 \end{array}] [\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ ⋮ \\ ⋮ \\ x_{21} \end{array}] = [\begin{array}{c} {10 x}_{1} + x_{2} \\ {2 x}_{1} + 9 x_{2} + x_{3} \\ ⋮ \\ ⋮ \\ {2 x}_{19} + 9 x_{20} + x_{21} \\ {2 x}_{20} + 10 x_{21} \end{array}]$ .

$A^{T} x = [\begin{array}{c} {10 x}_{1} + x_{2} \\ {2 x}_{1} + 9 x_{2} + x_{3} \\ ⋮ \\ ⋮ \\ {2 x}_{19} + 9 x_{20} + x_{21} \\ {2 x}_{20} + 10 x_{21} \end{array}] = 2 \cdot [\begin{array}{c} 0 \\ x_{1} \\ x_{2} \\ ⋮ \\ x_{20} \end{array}] + [\begin{array}{c} {10 x}_{1} \\ {9 x}_{2} \\ ⋮ \\ 9 x_{20} \\ 10 x_{21} \end{array}] + [\begin{array}{c} x_{2} \\ x_{3} \\ ⋮ \\ x_{21} \\ 0 \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

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

Теперь решите линейную систему $Ax = b$ путем предоставления lsqr с указателем на функцию, который вычисляет A*x и A'*x. Используйте допуск 1e-6 и 25 итераций. Определить $b$ как суммы строк $A$ так что истинное решение для $x$ является вектором таковых.

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

lsqr converged at iteration 21 to a solution with relative residual 5.4e-13.

Локальные функции

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

Входные параметры

`A` - Матрица коэффициентов
матрица | указатель на функцию

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

Определение `A` как указатель на функцию

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

Чтобы использовать указатель на функцию, используйте сигнатуру функции function y = afun(x,opt). Параметризация функций объясняет, как предоставить дополнительные параметры функции afun, при необходимости. Функция afun должны удовлетворять следующим условиям:

afun(x,'notransp') возвращает A*x продукта.
afun(x,'transp') возвращает A'*x продукта.

Примером приемлемой функции является:

function y = afun(x,opt,B,C,n)
if strcmp(opt,'notransp')
    y = [B*x(n+1:end); C*x(1:n)];
else
    y = [C'*x(n+1:end); B'*x(1:n)];
end

Функция afun использует значения в B и C для вычисления любой из A*x или A'*x (в зависимости от заданного флага) без фактического формирования всей матрицы.

Типы данных: double | function_handle
Поддержка комплексного числа: Да

`b` - Правая сторона линейного уравнения
Вектор-столбец

Правая сторона линейного уравнения, заданная как вектор-столбец. Длина b должно быть равно size(A,1).

Типы данных: double
Поддержка комплексного числа: Да

`tol` - Допуск метода
`[]` или `1e-6` (по умолчанию) | положительная скалярная величина

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

Типы данных: double

`maxit` - Максимальное количество итераций
`[]` или `min(size(A,1),20)` (по умолчанию) | положительное скалярное целое число

Максимальное количество итераций, заданное как положительное скалярное целое число. Увеличьте значение maxit чтобы разрешить больше итераций для lsqr для достижения допуска tol. Обычно меньшее значение tol означает, что для успешного завершения вычисления требуется больше итераций.

`M`, `M1`, `M2` - Матрицы Preconditioner (как отдельные аргументы)
`eye(size(A))` (по умолчанию) | матрицы | указатели на функцию

Матрицы Preconditioner, заданные как отдельные аргументы матриц или указателей на функцию. Можно задать матрицу предварительной подготовки M или его матричные множители M = M1*M2 улучшить числовые аспекты линейной системы и облегчить ее lsqr быстро сходиться. Для матриц квадратных коэффициентов можно использовать неполные функции матричного факторизации ilu и ichol чтобы сгенерировать матрицы preconditioner. Вы также можете использовать equilibrate перед факторизацией, чтобы улучшить число обусловленности матрицы коэффициентов. Для получения дополнительной информации о предварительных кондиционерах см. «Итерационные методы для линейных систем».

lsqr рассматривает неопределенные предварительные кондиционеры как единичные матрицы.

Определение `M` как указатель на функцию

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

Чтобы использовать указатель на функцию, сначала создайте функцию с подписью function y = mfun(x,opt). Параметризация функций объясняет, как предоставить дополнительные параметры функции mfun, при необходимости. Функция mfun должны удовлетворять следующим условиям:

mfun(x,'notransp') возвращает значение M\x или M2\(M1\x).
mfun(x,'transp') возвращает значение M'\x или M1'\(M2'\x).

Примером приемлемой функции является:

function y = mfun(x,opt,a,b)  
if strcmp(opt,'notransp')
    y = x.*a;
else
    y = x.*b;
end
end

В этом примере функция mfun использует a и b для вычисления любой из M\x = x*a или M'\x = x*b (в зависимости от заданного флага) без фактического формирования всей матрицы M.

Типы данных: double | function_handle
Поддержка комплексного числа: Да

`x0` - Начальное предположение
`[]` или вектор-столбец из нулей (по умолчанию) | вектор-столбец

Начальное предположение, заданное как вектор-столбец с длиной, равной size(A,2). Если вы можете предоставить lsqr с более разумным исходным предположением x0 чем нулевой вектор по умолчанию, то он может сэкономить время расчета и помочь алгоритму сходиться быстрее.

Типы данных: double
Поддержка комплексного числа: Да

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

`x` - Решение линейной системы
Вектор-столбец

Решение линейной системы, возвращаемое как вектор-столбец. Этот выход дает приблизительное решение линейной системы A*x = b.

Если flag является 0 и relres <= tol, затем x является последовательным решением для A*x = b.
Если flag является 0 но relres > tol, затем x - решение методом наименьших квадратов, которое минимизирует norm(b-A*x). В этом случае lsvec выход содержит масштабированную ошибку нормального уравнения x.

Всякий раз, когда расчет не успешен (flag ~= 0), решение x возвращено lsqr - это тот, с минимальной нормой невязки, вычисленной по всем итерациям.

`flag` - Флаг сходимости
скаляр

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

Значение флага	Сходимость
`0`	Успех - `lsqr` сходился к желаемому допуску `tol` в пределах `maxit` итераций.
`1`	Отказ - `lsqr` итерация `maxit` итерации, но не сходились.
`2`	Отказ - матрица предварительной подготовки `M` или `M = M1*M2` является плохо обусловленной.
`3`	Отказ - `lsqr` застопорившиеся после двух последовательных итераций были одинаковыми.
`4`	Отказ - Одна из скалярных величин, вычисленных `lsqr` алгоритм стал слишком маленьким или слишком большим, чтобы продолжить вычисления.

`relres` - Относительная остаточная ошибка
скаляр

Относительная остаточная ошибка, возвращенная как скаляр. Относительная остаточная ошибка является показателем того, насколько точен возвращенный ответ x является. lsqr отслеживает относительную невязку и невязку методом наименьших квадратов при каждой итерации в процессе решения, и алгоритм сходится, когда любая невязка соответствует заданному допуску tol. The relres выход содержит значение сходящейся невязки, относительной невязки или невязки методом наименьших квадратов:

Относительная остаточная ошибка равна norm(b-A*x)/norm(b) и является ли обычно невязка, которая соответствует допуску tol когда lsqr сходится. The resvec выход отслеживает историю этого невязки по всем итерациям.
Ошибка невязки методом наименьших квадратов равна norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro'). Эта остаточная причина lsqr сходиться реже, чем относительная невязка. The lsvec выход отслеживает историю этого невязки по всем итерациям.

`iter` - Число итерации
скаляр

Число итерации, возвращенное как скаляр. Этот выход указывает номер итерации, при котором вычисленный ответ для x была рассчитана.

Типы данных: double

`resvec` - Остаточная ошибка
вектор

Остаточная ошибка, возвращенная как вектор. Остаточная ошибка norm(b-A*x) показывает, насколько близок алгоритм к сходимости для заданного значения x. Количество элементов в resvec равно количеству итераций. Можно изучить содержимое resvec чтобы помочь решить, изменять ли значения tol или maxit.

Типы данных: double

`lsvec` - Масштабированная ошибка нормального уравнения
вектор

Масштабированная нормальная ошибка уравнения, возвращенная как вектор. Для каждой итерации lsvec содержит оценку масштабированной невязки нормального уравнения norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro'). Количество элементов в lsvec равно количеству итераций.

Подробнее о