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

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 путем эффективного решения системы ${AM}^{−} 1y$ = b для y, $где$ y = Mx. Использование матрицы предварительного кондиционирования может улучшить числовые свойства задачи и эффективность вычисления.

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 дает относительный остаток, меньший, чем предписанный допуск 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:

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, и 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 \begin{array}{c} x=[1020⋯⋯01920⋮01⋱20⋮010⋱⋱⋮0⋱1⋱0⋮⋱⋱⋱20⋯⋯0110 \end{array}] \begin{array}{c} [_{} \\ _{} \\ _{} \\ _{} \end{array} \begin{array}{c} _{}_{} \\ _{}_{}_{} \\ _{}_{}_{} \\ _{}_{} \end{array} x1x2x3⋮⋮x21]=[10x1+2x2x1+9x2+2x3⋮⋮x19+9x20+2x21x20+10x21$ ].

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

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

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

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

$^{AT} \begin{array}{c} _{}_{} \\ _{}_{}_{} \\ _{}_{}_{} \\ _{}_{} \end{array} \begin{array}{c} _{} \\ _{} \\ _{} \end{array} \begin{array}{c} _{} \\ _{} \\ _{} \\ _{} \end{array} \begin{array}{c} _{} \\ _{} \\ _{} \\ x=[10x1+x22x1+9x2+x3⋮⋮2x19+9x20+x212x20+10x21]=2⋅[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

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

Теперь решите линейную систему $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` - Матрицы предварительного условия (как отдельные аргументы)
`eye(size(A))` (по умолчанию) | матрицы | дескрипторы функций

Матрицы предварительного условия, определенные как отдельные аргументы матриц или дескрипторов функций. Можно указать матрицу предварительного условия. M или его матричные коэффициенты M = M1*M2 улучшить численные аспекты линейной системы и упростить ее lsqr чтобы быстро сойтись. Для матриц с квадратными коэффициентами можно использовать неполные функции факторизации матриц ilu и ichol для формирования матриц предварительного условия. Вы также можете использовать 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. relres вывод содержит значение сходящегося остатка, либо относительного остатка, либо остатка наименьших квадратов:

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

Подробнее