scipy.sparse.linalg.

lsqr#

scipy.sparse.linalg.lsqr(A, b, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, iter_lim=None, показать=False, calc_var=False, x0=None)[источник]#

Найти решение методом наименьших квадратов для большой, разреженной, линейной системы уравнений.

Функция решает Ax = b или min ||Ax - b||^2 или min ||Ax - b||^2 + d^2 ||x - x0||^2.

Матрица A может быть квадратной или прямоугольной (переопределённой или недоопределённой) и может иметь любой ранг.

1. Unsymmetric equations --    solve  Ax = b

2. Linear least squares  --    solve  Ax = b
                               in the least-squares sense

3. Damped least squares  --    solve  (   A    )*x = (    b    )
                                      ( damp*I )     ( damp*x0 )
                               in the least-squares sense

Параметры:

A{sparse array, ndarray, LinearOperator}: Представление матрицы размером m на n. Альтернативно, A может быть линейным оператором, который может производить Ax и A^T x используя, например, scipy.sparse.linalg.LinearOperator.
barray_like, shape (m,): Вектор правой части b.
dampfloat: Коэффициент демпфирования. По умолчанию 0.
atol, btolfloat, опционально: Допуски остановки. lsqr продолжает итерации до тех пор, пока некоторый оценка обратной ошибки не станет меньше некоторой величины, зависящей от atol и btol. Пусть r = b - Ax будет вектором невязки для текущего приближенного решения x. Если Ax = b кажется последовательным, lsqr завершается, когда norm(r) <= atol * norm(A) * norm(x) + btol * norm(b). В противном случае, lsqr завершается, когда norm(A^H r) <= atol * norm(A) * norm(r). Если оба допуска равны 1.0e-6 (по умолчанию), итоговый norm(r) должна быть точной примерно до 6 знаков. (Конечный x обычно будет иметь меньше верных цифр, в зависимости от cond(A) и размер LAMBDA.) Если atol или btol если None, будет использовано значение по умолчанию 1.0e-6. В идеале они должны быть оценками относительной ошибки в элементах A и b соответственно. Например, если элементы из A иметь 7 верных цифр, установить atol = 1e-7. Это предотвращает выполнение алгоритмом лишней работы за пределами неопределённости входных данных.
conlimfloat, опционально: Еще один критерий остановки. lsqr завершает работу, если оценка cond(A) превышает conlim. Для совместимых систем Ax = b, conlim может достигать 1.0e+12 (например). Для задач метода наименьших квадратов conlim должен быть меньше 1.0e+8. Максимальная точность может быть достигнута установкой atol = btol = conlim = zero, но количество итераций может быть чрезмерным. По умолчанию 1e8.
iter_limint, необязательный: Явное ограничение на количество итераций (для безопасности).
показатьbool, необязательно: Отображать журнал итераций. По умолчанию False.
calc_varbool, необязательно: Оценивать ли диагонали (A'A + damp^2*I)^{-1}.
x0array_like, shape (n,), optional: Начальное предположение для x, если None, используются нули. По умолчанию None.

Добавлено в версии 1.0.0.

Возвращает:

xndarray из float: Финальное решение.
istopint: Указывает причину завершения. 1 означает, что x является приближенным решением Ax = b. 2 означает, что x приближенно решает задачу наименьших квадратов.
itnint: Номер итерации при завершении.
r1normfloat: norm(r), где r = b - Ax.
r2normfloat: sqrt( norm(r)^2 + damp^2 * norm(x - x0)^2 ). Равно r1norm if damp == 0.
anormfloat: Оценка нормы Фробениуса для Abar = [[A]; [damp*I]].
acondfloat: Оценка cond(Abar).
arnormfloat: Оценка norm(A'@r - damp^2*(x - x0)).
xnormfloat: norm(x)
varndarray из float: Если calc_var is True, оценивает все диагонали (A'A)^{-1} (если damp == 0) или в более общем случае (A'A + damp^2*I)^{-1}. Это хорошо определено, если A имеет полный ранг по столбцам или damp > 0. (Не уверен, что означает var, если rank(A) < n и damp = 0.)

Примечания

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

Например, в задаче 1 решение не изменяется при масштабировании строк. Если строка A очень мала или велика по сравнению с другими строками A, соответствующая строка ( A b ) должна быть масштабирована вверх или вниз.

В задачах 1 и 2 решение x легко восстановить после масштабирования по столбцам. Если нет лучшей информации, ненулевые столбцы A должны быть масштабированы так, чтобы все они имели одинаковую евклидову норму (например, 1.0).

В задаче 3 нет свободы для изменения масштаба, если damp ненулевой. Однако значение damp должно назначаться только после внимательного рассмотрения масштабирования A.

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

Если какая-то начальная оценка x0 известно и если damp == 0, можно действовать следующим образом:

Вычислить вектор остатков r0 = b - A@x0.

Используйте LSQR для решения системы A@dx = r0.

Добавить поправку dx для получения окончательного решения x = x0 + dx.

Это требует, чтобы x0 должны быть доступны до и после вызова LSQR. Чтобы оценить преимущества, предположим, что LSQR требует k1 итераций для решения A@x = b и k2 итераций для решения A@dx = r0. Если x0 «хороший», norm(r0) будет меньше norm(b). Если одинаковые допуски остановки atol и btol используются для каждой системы, k1 и k2 будут похожи, но итоговое решение x0 + dx должно быть более точным. Единственный способ уменьшить общую работу - использовать больший допуск остановки для второй системы. Если некоторое значение btol подходит для A@x = b, большее значение btol*norm(b)/norm(r0) должно подходить для A@dx = r0.

Предобусловливание — ещё один способ уменьшить количество итераций. Если возможно решить связанную систему M@x = b эффективно, где M аппроксимирует A полезным образом (например, M - A имеет низкий ранг или его элементы малы относительно элементов A), LSQR может сходиться быстрее на системе A@M(inverse)@z = b, после чего x можно восстановить, решив M@x = z.

Если A симметрична, LSQR не следует использовать!

Альтернативами являются симметричный метод сопряжённых градиентов (cg) и/или SYMMLQ. SYMMLQ - это реализация симметричного cg, которая применяется к любой симметричной A и будет сходиться быстрее, чем LSQR. Если A положительно определена, существуют другие реализации симметричного cg, которые требуют немного меньше работы на итерацию, чем SYMMLQ (но займут то же количество итераций).

Ссылки

[1]

C. C. Paige и M. A. Saunders (1982a). “LSQR: Алгоритм для разреженных линейных уравнений и разреженных наименьших квадратов”, ACM TOMS 8(1), 43-71.

[2]

C. C. Paige and M. A. Saunders (1982b). “Algorithm 583. LSQR: Sparse linear equations and least squares problems”, ACM TOMS 8(2), 195-209.

[3]

M. A. Saunders (1995). «Решение разреженных прямоугольных систем с использованием LSQR и CRAIG», BIT 35, 588-604.

Примеры

>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import lsqr
>>> A = csc_array([[1., 0.], [1., 1.], [0., 1.]], dtype=float)

Первый пример имеет тривиальное решение [0, 0]

>>> b = np.array([0., 0., 0.], dtype=float)
>>> x, istop, itn, normr = lsqr(A, b)[:4]
>>> istop
0
>>> x
array([ 0.,  0.])

Код остановки istop=0 возвращённое значение указывает, что вектор нулей был найден как решение. Возвращённое решение x действительно содержит [0., 0.]. Следующий пример имеет нетривиальное решение:

>>> b = np.array([1., 0., -1.], dtype=float)
>>> x, istop, itn, r1norm = lsqr(A, b)[:4]
>>> istop
1
>>> x
array([ 1., -1.])
>>> itn
1
>>> r1norm
4.440892098500627e-16

Как указано istop=1, lsqr найдено решение, удовлетворяющее пределам допуска. Данное решение [1., -1.] очевидно решает уравнение. Остальные возвращаемые значения включают информацию о количестве итераций (itn=1) и оставшаяся разница левой и правой частей решенного уравнения. Последний пример демонстрирует поведение в случае, когда нет решения уравнения:

>>> b = np.array([1., 0.01, -1.], dtype=float)
>>> x, istop, itn, r1norm = lsqr(A, b)[:4]
>>> istop
2
>>> x
array([ 1.00333333, -0.99666667])
>>> A.dot(x)-b
array([ 0.00333333, -0.00333333,  0.00333333])
>>> r1norm
0.005773502691896255

istop указывает, что система несовместна и, следовательно, x скорее приближённое решение соответствующей задачи наименьших квадратов. r1norm содержит норму минимального остатка, который был найден.