lsqr

Решите систему линейных уравнений — метод наименьших квадратов

Описание

пример

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 путем эффективного решения системы AM1y=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. Выходное отображение включает значение относительной остаточной ошибки b-Axb.

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 object. The axes object contains 2 objects of type line. These objects represent Least-squares residual, Relative residual.

Исследуйте эффект использования матрицы перед формирователем с 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=LU в разложенной на множители форме. Задайте допуск отбрасывания, чтобы проигнорировать недиагональные записи со значениями, меньшими, чем 1e-6. Решите предобусловленную систему AM-1(Mx)=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 object. The axes object 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.

Выражение Ax становится:

Ax=[1020019200120010010200110][x1x2x3x21]=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21].

Итоговый вектор может быть записан как сумма трех векторов:

Ax=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21]=[0x1x2x20]+[10x19x29x2010x21]+2[x2x3x210].

Аналогично, выражение для ATx становится:

ATx=[1010029100210020010100210][x1x2x3x21]=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21].

ATx=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21]=2[0x1x2x20]+[10x19x29x2010x21]+[x2x3x210].

В 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.
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

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

свернуть все

Матрица коэффициентов в виде матрицы или указателя на функцию. Эта матрица является матрицей коэффициентов в линейной системе 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 должно быть равно size(A,1).

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

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

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

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

Матрицы перед формирователем в виде отдельных аргументов матриц или указателей на функцию. Можно задать матрицу перед формирователем 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
Поддержка комплексного числа: Да

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

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

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

свернуть все

Решение для линейной системы, возвращенное как вектор-столбец. Этот выход дает приближенное решение линейной системы 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 тот с минимальной невязкой нормы, вычисленной по всем итерациям.

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

Флаговое значение

Сходимость

0

Успех — lsqr сходившийся к желаемому допуску tol в maxit итерации.

1

Отказ — lsqr выполненный с помощью итераций maxit итерации, но не сходились.

2

Отказ — матрица перед формирователем M или M = M1*M2 isIllConditioned.

3

Отказ — lsqr застоявшийся после того, как две последовательных итерации были тем же самым.

4

Отказ — Один из скаляров, вычисленных lsqr алгоритм стал слишком маленьким или слишком большим, чтобы продолжить вычислять.

Относительная остаточная ошибка, возвращенная как скаляр. Относительная остаточная ошибка является индикацией относительно как точный данный ответ 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 выведите отслеживает историю этой невязки по всем итерациям.

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

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

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

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

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

Больше о

свернуть все

Метод наименьших квадратов

Наименьшие квадраты (LSQR) алгоритм являются адаптацией метода методов сопряженных градиентов (CG) для прямоугольных матриц. Аналитически, LSQR для A*x = b производит те же остаточные значения как CG для нормальных уравнений A'*A*x = A'*b, но LSQR обладает более благоприятными числовыми свойствами и таким образом обычно более надежен [1].

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

Советы

  • Сходимость большинства итерационных методов зависит от числа обусловленности матрицы коэффициентов, cond(A). Можно использовать equilibrate улучшить число обусловленности A, и самостоятельно это облегчает для большинства итеративных решателей сходиться. Однако использование equilibrate также приводит к лучшим качественным матрицам перед формирователем, когда вы впоследствии учитываете уравновешенный матричный B = R*P*A*C.

  • Можно использовать матричные функции переупорядочения такой как dissect и symrcm переставить строки и столбцы матрицы коэффициентов и минимизировать количество ненулей, когда матрица коэффициентов учтена, чтобы сгенерировать предварительный формирователь. Это может уменьшать память и время, требуемое впоследствии решить предобусловленную линейную систему.

Ссылки

[1] Барретт, R., М. Берри, Т. Ф. Чан, и др., Шаблоны для Решения Линейных систем: Базовые блоки для Итерационных методов, SIAM, Филадельфия, 1994.

[2] Пэйдж, C. C. и М. А. Сондерс, "LSQR: Алгоритм для Разреженных Линейных уравнений И Разреженных Наименьших квадратов", Математика Сделки ACM. Мягкий., Vol.8, 1982, стр 43-71.

Расширенные возможности

Представлено до R2006a