import warnings
import numpy as np
from scipy.sparse.linalg._interface import LinearOperator
from .utils import make_system
from scipy.linalg import get_lapack_funcs
__all__ = ['bicg', 'bicgstab', 'cg', 'cgs', 'gmres', 'qmr']
def _get_atol_rtol(name, b_norm, atol=0., rtol=1e-5):
"""
A helper function to handle tolerance normalization
"""
if atol == 'legacy' or atol is None or atol < 0:
msg = (f"'scipy.sparse.linalg.{name}' called with invalid `atol`={atol}; "
"if set, `atol` must be a real, non-negative number.")
raise ValueError(msg)
atol = max(float(atol), float(rtol) * float(b_norm))
return atol, rtol
def bicg(A, b, x0=None, *, rtol=1e-5, atol=0., maxiter=None, M=None, callback=None):
"""Use BIConjugate Gradient iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse array, ndarray, LinearOperator}
The real or complex N-by-N matrix of the linear system.
Alternatively, `A` can be a linear operator which can
produce ``Ax`` and ``A^T x`` using, e.g.,
``scipy.sparse.linalg.LinearOperator``.
b : ndarray
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : ndarray
Starting guess for the solution.
rtol, atol : float, optional
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``atol=0.`` and ``rtol=1e-5``.
maxiter : integer
Maximum number of iterations. Iteration will stop after maxiter
steps even if the specified tolerance has not been achieved.
M : {sparse array, ndarray, LinearOperator}
Preconditioner for `A`. It should approximate the
inverse of `A` (see Notes). Effective preconditioning dramatically improves the
rate of convergence, which implies that fewer iterations are needed
to reach a given error tolerance.
callback : function
User-supplied function to call after each iteration. It is called
as ``callback(xk)``, where ``xk`` is the current solution vector.
Returns
-------
x : ndarray
The converged solution.
info : integer
Provides convergence information:
0 : successful exit
>0 : convergence to tolerance not achieved, number of iterations
<0 : parameter breakdown
Notes
-----
The preconditioner `M` should be a matrix such that ``M @ A`` has a smaller
condition number than `A`, see [1]_ .
References
----------
.. [1] "Preconditioner", Wikipedia,
https://en.wikipedia.org/wiki/Preconditioner
.. [2] "Biconjugate gradient method", Wikipedia,
https://en.wikipedia.org/wiki/Biconjugate_gradient_method
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import bicg
>>> A = csc_array([[3, 2, 0], [1, -1, 0], [0, 5, 1.]])
>>> b = np.array([2., 4., -1.])
>>> x, exitCode = bicg(A, b, atol=1e-5)
>>> print(exitCode) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
A, M, x, b, postprocess = make_system(A, M, x0, b)
bnrm2 = np.linalg.norm(b)
atol, _ = _get_atol_rtol('bicg', bnrm2, atol, rtol)
if bnrm2 == 0:
return postprocess(b), 0
n = len(b)
dotprod = np.vdot if np.iscomplexobj(x) else np.dot
if maxiter is None:
maxiter = n*10
matvec, rmatvec = A.matvec, A.rmatvec
psolve, rpsolve = M.matvec, M.rmatvec
rhotol = np.finfo(x.dtype.char).eps**2
# Dummy values to initialize vars, silence linter warnings
rho_prev, p, ptilde = None, None, None
r = b - matvec(x) if x.any() else b.copy()
rtilde = r.copy()
for iteration in range(maxiter):
if np.linalg.norm(r) < atol: # Are we done?
return postprocess(x), 0
z = psolve(r)
ztilde = rpsolve(rtilde)
# order matters in this dot product
rho_cur = dotprod(rtilde, z)
if np.abs(rho_cur) < rhotol: # Breakdown case
return postprocess, -10
if iteration > 0:
beta = rho_cur / rho_prev
p *= beta
p += z
ptilde *= beta.conj()
ptilde += ztilde
else: # First spin
p = z.copy()
ptilde = ztilde.copy()
q = matvec(p)
qtilde = rmatvec(ptilde)
rv = dotprod(ptilde, q)
if rv == 0:
return postprocess(x), -11
alpha = rho_cur / rv
x += alpha*p
r -= alpha*q
rtilde -= alpha.conj()*qtilde
rho_prev = rho_cur
if callback:
callback(x)
else: # for loop exhausted
# Return incomplete progress
return postprocess(x), maxiter
def bicgstab(A, b, x0=None, *, rtol=1e-5, atol=0., maxiter=None, M=None,
callback=None):
"""Use BIConjugate Gradient STABilized iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse array, ndarray, LinearOperator}
The real or complex N-by-N matrix of the linear system.
Alternatively, `A` can be a linear operator which can
produce ``Ax`` and ``A^T x`` using, e.g.,
``scipy.sparse.linalg.LinearOperator``.
b : ndarray
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : ndarray
Starting guess for the solution.
rtol, atol : float, optional
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``atol=0.`` and ``rtol=1e-5``.
maxiter : integer
Maximum number of iterations. Iteration will stop after maxiter
steps even if the specified tolerance has not been achieved.
M : {sparse array, ndarray, LinearOperator}
Preconditioner for `A`. It should approximate the
inverse of `A` (see Notes). Effective preconditioning dramatically improves the
rate of convergence, which implies that fewer iterations are needed
to reach a given error tolerance.
callback : function
User-supplied function to call after each iteration. It is called
as ``callback(xk)``, where ``xk`` is the current solution vector.
Returns
-------
x : ndarray
The converged solution.
info : integer
Provides convergence information:
0 : successful exit
>0 : convergence to tolerance not achieved, number of iterations
<0 : parameter breakdown
Notes
-----
The preconditioner `M` should be a matrix such that ``M @ A`` has a smaller
condition number than `A`, see [1]_ .
References
----------
.. [1] "Preconditioner", Wikipedia,
https://en.wikipedia.org/wiki/Preconditioner
.. [2] "Biconjugate gradient stabilized method",
Wikipedia, https://en.wikipedia.org/wiki/Biconjugate_gradient_stabilized_method
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import bicgstab
>>> R = np.array([[4, 2, 0, 1],
... [3, 0, 0, 2],
... [0, 1, 1, 1],
... [0, 2, 1, 0]])
>>> A = csc_array(R)
>>> b = np.array([-1, -0.5, -1, 2])
>>> x, exit_code = bicgstab(A, b, atol=1e-5)
>>> print(exit_code) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
A, M, x, b, postprocess = make_system(A, M, x0, b)
bnrm2 = np.linalg.norm(b)
atol, _ = _get_atol_rtol('bicgstab', bnrm2, atol, rtol)
if bnrm2 == 0:
return postprocess(b), 0
n = len(b)
dotprod = np.vdot if np.iscomplexobj(x) else np.dot
if maxiter is None:
maxiter = n*10
matvec = A.matvec
psolve = M.matvec
# These values make no sense but coming from original Fortran code
# sqrt might have been meant instead.
rhotol = np.finfo(x.dtype.char).eps**2
omegatol = rhotol
# Dummy values to initialize vars, silence linter warnings
rho_prev, omega, alpha, p, v = None, None, None, None, None
r = b - matvec(x) if x.any() else b.copy()
rtilde = r.copy()
for iteration in range(maxiter):
if np.linalg.norm(r) < atol: # Are we done?
return postprocess(x), 0
rho = dotprod(rtilde, r)
if np.abs(rho) < rhotol: # rho breakdown
return postprocess(x), -10
if iteration > 0:
if np.abs(omega) < omegatol: # omega breakdown
return postprocess(x), -11
beta = (rho / rho_prev) * (alpha / omega)
p -= omega*v
p *= beta
p += r
else: # First spin
s = np.empty_like(r)
p = r.copy()
phat = psolve(p)
v = matvec(phat)
rv = dotprod(rtilde, v)
if rv == 0:
return postprocess(x), -11
alpha = rho / rv
r -= alpha*v
s[:] = r[:]
if np.linalg.norm(s) < atol:
x += alpha*phat
return postprocess(x), 0
shat = psolve(s)
t = matvec(shat)
omega = dotprod(t, s) / dotprod(t, t)
x += alpha*phat
x += omega*shat
r -= omega*t
rho_prev = rho
if callback:
callback(x)
else: # for loop exhausted
# Return incomplete progress
return postprocess(x), maxiter
def cg(A, b, x0=None, *, rtol=1e-5, atol=0., maxiter=None, M=None, callback=None):
"""Use Conjugate Gradient iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse array, ndarray, LinearOperator}
The real or complex N-by-N matrix of the linear system.
`A` must represent a hermitian, positive definite matrix.
Alternatively, `A` can be a linear operator which can
produce ``Ax`` using, e.g.,
``scipy.sparse.linalg.LinearOperator``.
b : ndarray
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : ndarray
Starting guess for the solution.
rtol, atol : float, optional
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``atol=0.`` and ``rtol=1e-5``.
maxiter : integer
Maximum number of iterations. Iteration will stop after maxiter
steps even if the specified tolerance has not been achieved.
M : {sparse array, ndarray, LinearOperator}
Preconditioner for `A`. `M` must represent a hermitian, positive definite
matrix. It should approximate the inverse of `A` (see Notes).
Effective preconditioning dramatically improves the
rate of convergence, which implies that fewer iterations are needed
to reach a given error tolerance.
callback : function
User-supplied function to call after each iteration. It is called
as ``callback(xk)``, where ``xk`` is the current solution vector.
Returns
-------
x : ndarray
The converged solution.
info : integer
Provides convergence information:
0 : successful exit
>0 : convergence to tolerance not achieved, number of iterations
Notes
-----
The preconditioner `M` should be a matrix such that ``M @ A`` has a smaller
condition number than `A`, see [2]_.
References
----------
.. [1] "Conjugate Gradient Method, Wikipedia,
https://en.wikipedia.org/wiki/Conjugate_gradient_method
.. [2] "Preconditioner",
Wikipedia, https://en.wikipedia.org/wiki/Preconditioner
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import cg
>>> P = np.array([[4, 0, 1, 0],
... [0, 5, 0, 0],
... [1, 0, 3, 2],
... [0, 0, 2, 4]])
>>> A = csc_array(P)
>>> b = np.array([-1, -0.5, -1, 2])
>>> x, exit_code = cg(A, b, atol=1e-5)
>>> print(exit_code) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
A, M, x, b, postprocess = make_system(A, M, x0, b)
bnrm2 = np.linalg.norm(b)
atol, _ = _get_atol_rtol('cg', bnrm2, atol, rtol)
if bnrm2 == 0:
return postprocess(b), 0
n = len(b)
if maxiter is None:
maxiter = n*10
dotprod = np.vdot if np.iscomplexobj(x) else np.dot
matvec = A.matvec
psolve = M.matvec
r = b - matvec(x) if x.any() else b.copy()
# Dummy value to initialize var, silences warnings
rho_prev, p = None, None
for iteration in range(maxiter):
if np.linalg.norm(r) < atol: # Are we done?
return postprocess(x), 0
z = psolve(r)
rho_cur = dotprod(r, z)
if iteration > 0:
beta = rho_cur / rho_prev
p *= beta
p += z
else: # First spin
p = np.empty_like(r)
p[:] = z[:]
q = matvec(p)
alpha = rho_cur / dotprod(p, q)
x += alpha*p
r -= alpha*q
rho_prev = rho_cur
if callback:
callback(x)
else: # for loop exhausted
# Return incomplete progress
return postprocess(x), maxiter
def cgs(A, b, x0=None, *, rtol=1e-5, atol=0., maxiter=None, M=None, callback=None):
"""Use Conjugate Gradient Squared iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse array, ndarray, LinearOperator}
The real-valued N-by-N matrix of the linear system.
Alternatively, `A` can be a linear operator which can
produce ``Ax`` using, e.g.,
``scipy.sparse.linalg.LinearOperator``.
b : ndarray
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : ndarray
Starting guess for the solution.
rtol, atol : float, optional
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``atol=0.`` and ``rtol=1e-5``.
maxiter : integer
Maximum number of iterations. Iteration will stop after maxiter
steps even if the specified tolerance has not been achieved.
M : {sparse array, ndarray, LinearOperator}
Preconditioner for ``A``. It should approximate the
inverse of `A` (see Notes). Effective preconditioning dramatically improves the
rate of convergence, which implies that fewer iterations are needed
to reach a given error tolerance.
callback : function
User-supplied function to call after each iteration. It is called
as ``callback(xk)``, where ``xk`` is the current solution vector.
Returns
-------
x : ndarray
The converged solution.
info : integer
Provides convergence information:
0 : successful exit
>0 : convergence to tolerance not achieved, number of iterations
<0 : parameter breakdown
Notes
-----
The preconditioner `M` should be a matrix such that ``M @ A`` has a smaller
condition number than `A`, see [1]_.
References
----------
.. [1] "Preconditioner", Wikipedia,
https://en.wikipedia.org/wiki/Preconditioner
.. [2] "Conjugate gradient squared", Wikipedia,
https://en.wikipedia.org/wiki/Conjugate_gradient_squared_method
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import cgs
>>> R = np.array([[4, 2, 0, 1],
... [3, 0, 0, 2],
... [0, 1, 1, 1],
... [0, 2, 1, 0]])
>>> A = csc_array(R)
>>> b = np.array([-1, -0.5, -1, 2])
>>> x, exit_code = cgs(A, b)
>>> print(exit_code) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
A, M, x, b, postprocess = make_system(A, M, x0, b)
bnrm2 = np.linalg.norm(b)
atol, _ = _get_atol_rtol('cgs', bnrm2, atol, rtol)
if bnrm2 == 0:
return postprocess(b), 0
n = len(b)
dotprod = np.vdot if np.iscomplexobj(x) else np.dot
if maxiter is None:
maxiter = n*10
matvec = A.matvec
psolve = M.matvec
rhotol = np.finfo(x.dtype.char).eps**2
r = b - matvec(x) if x.any() else b.copy()
rtilde = r.copy()
bnorm = np.linalg.norm(b)
if bnorm == 0:
bnorm = 1
# Dummy values to initialize vars, silence linter warnings
rho_prev, p, u, q = None, None, None, None
for iteration in range(maxiter):
rnorm = np.linalg.norm(r)
if rnorm < atol: # Are we done?
return postprocess(x), 0
rho_cur = dotprod(rtilde, r)
if np.abs(rho_cur) < rhotol: # Breakdown case
return postprocess, -10
if iteration > 0:
beta = rho_cur / rho_prev
# u = r + beta * q
# p = u + beta * (q + beta * p);
u[:] = r[:]
u += beta*q
p *= beta
p += q
p *= beta
p += u
else: # First spin
p = r.copy()
u = r.copy()
q = np.empty_like(r)
phat = psolve(p)
vhat = matvec(phat)
rv = dotprod(rtilde, vhat)
if rv == 0: # Dot product breakdown
return postprocess(x), -11
alpha = rho_cur / rv
q[:] = u[:]
q -= alpha*vhat
uhat = psolve(u + q)
x += alpha*uhat
# Due to numerical error build-up the actual residual is computed
# instead of the following two lines that were in the original
# FORTRAN templates, still using a single matvec.
# qhat = matvec(uhat)
# r -= alpha*qhat
r = b - matvec(x)
rho_prev = rho_cur
if callback:
callback(x)
else: # for loop exhausted
# Return incomplete progress
return postprocess(x), maxiter
def gmres(A, b, x0=None, *, rtol=1e-5, atol=0., restart=None, maxiter=None, M=None,
callback=None, callback_type=None):
"""
Use Generalized Minimal RESidual iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse array, ndarray, LinearOperator}
The real or complex N-by-N matrix of the linear system.
Alternatively, `A` can be a linear operator which can
produce ``Ax`` using, e.g.,
``scipy.sparse.linalg.LinearOperator``.
b : ndarray
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : ndarray
Starting guess for the solution (a vector of zeros by default).
atol, rtol : float
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``atol=0.`` and ``rtol=1e-5``.
restart : int, optional
Number of iterations between restarts. Larger values increase
iteration cost, but may be necessary for convergence.
If omitted, ``min(20, n)`` is used.
maxiter : int, optional
Maximum number of iterations (restart cycles). Iteration will stop
after maxiter steps even if the specified tolerance has not been
achieved. See `callback_type`.
M : {sparse array, ndarray, LinearOperator}
Inverse of the preconditioner of `A`. `M` should approximate the
inverse of `A` and be easy to solve for (see Notes). Effective
preconditioning dramatically improves the rate of convergence,
which implies that fewer iterations are needed to reach a given
error tolerance. By default, no preconditioner is used.
In this implementation, left preconditioning is used,
and the preconditioned residual is minimized. However, the final
convergence is tested with respect to the ``b - A @ x`` residual.
callback : function
User-supplied function to call after each iteration. It is called
as ``callback(args)``, where ``args`` are selected by `callback_type`.
callback_type : {'x', 'pr_norm', 'legacy'}, optional
Callback function argument requested:
- ``x``: current iterate (ndarray), called on every restart
- ``pr_norm``: relative (preconditioned) residual norm (float),
called on every inner iteration
- ``legacy`` (default): same as ``pr_norm``, but also changes the
meaning of `maxiter` to count inner iterations instead of restart
cycles.
This keyword has no effect if `callback` is not set.
Returns
-------
x : ndarray
The converged solution.
info : int
Provides convergence information:
0 : successful exit
>0 : convergence to tolerance not achieved, number of iterations
See Also
--------
LinearOperator
Notes
-----
A preconditioner, P, is chosen such that P is close to A but easy to solve
for. The preconditioner parameter required by this routine is
``M = P^-1``. The inverse should preferably not be calculated
explicitly. Rather, use the following template to produce M::
# Construct a linear operator that computes P^-1 @ x.
import scipy.sparse.linalg as spla
M_x = lambda x: spla.spsolve(P, x)
M = spla.LinearOperator((n, n), M_x)
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import gmres
>>> A = csc_array([[3, 2, 0], [1, -1, 0], [0, 5, 1]], dtype=float)
>>> b = np.array([2, 4, -1], dtype=float)
>>> x, exitCode = gmres(A, b, atol=1e-5)
>>> print(exitCode) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
if callback is not None and callback_type is None:
# Warn about 'callback_type' semantic changes.
# Probably should be removed only in far future, Scipy 2.0 or so.
msg = ("scipy.sparse.linalg.gmres called without specifying "
"`callback_type`. The default value will be changed in"
" a future release. For compatibility, specify a value "
"for `callback_type` explicitly, e.g., "
"``gmres(..., callback_type='pr_norm')``, or to retain the "
"old behavior ``gmres(..., callback_type='legacy')``"
)
warnings.warn(msg, category=DeprecationWarning, stacklevel=3)
if callback_type is None:
callback_type = 'legacy'
if callback_type not in ('x', 'pr_norm', 'legacy'):
raise ValueError(f"Unknown callback_type: {callback_type!r}")
if callback is None:
callback_type = None
A, M, x, b, postprocess = make_system(A, M, x0, b)
matvec = A.matvec
psolve = M.matvec
n = len(b)
bnrm2 = np.linalg.norm(b)
atol, _ = _get_atol_rtol('gmres', bnrm2, atol, rtol)
if bnrm2 == 0:
return postprocess(b), 0
eps = np.finfo(x.dtype.char).eps
dotprod = np.vdot if np.iscomplexobj(x) else np.dot
if maxiter is None:
maxiter = n*10
if restart is None:
restart = 20
restart = min(restart, n)
Mb_nrm2 = np.linalg.norm(psolve(b))
# ====================================================
# =========== Tolerance control from gh-8400 =========
# ====================================================
# Tolerance passed to GMRESREVCOM applies to the inner
# iteration and deals with the left-preconditioned
# residual.
ptol_max_factor = 1.
ptol = Mb_nrm2 * min(ptol_max_factor, atol / bnrm2)
presid = 0.
# ====================================================
lartg = get_lapack_funcs('lartg', dtype=x.dtype)
# allocate internal variables
v = np.empty([restart+1, n], dtype=x.dtype)
h = np.zeros([restart, restart+1], dtype=x.dtype)
givens = np.zeros([restart, 2], dtype=x.dtype)
# legacy iteration count
inner_iter = 0
for iteration in range(maxiter):
if iteration == 0:
r = b - matvec(x) if x.any() else b.copy()
if np.linalg.norm(r) < atol: # Are we done?
return postprocess(x), 0
v[0, :] = psolve(r)
tmp = np.linalg.norm(v[0, :])
v[0, :] *= (1 / tmp)
# RHS of the Hessenberg problem
S = np.zeros(restart+1, dtype=x.dtype)
S[0] = tmp
breakdown = False
for col in range(restart):
av = matvec(v[col, :])
w = psolve(av)
# Modified Gram-Schmidt
h0 = np.linalg.norm(w)
for k in range(col+1):
tmp = dotprod(v[k, :], w)
h[col, k] = tmp
w -= tmp*v[k, :]
h1 = np.linalg.norm(w)
h[col, col + 1] = h1
v[col + 1, :] = w[:]
# Exact solution indicator
if h1 <= eps*h0:
h[col, col + 1] = 0
breakdown = True
else:
v[col + 1, :] *= (1 / h1)
# apply past Givens rotations to current h column
for k in range(col):
c, s = givens[k, 0], givens[k, 1]
n0, n1 = h[col, [k, k+1]]
h[col, [k, k + 1]] = [c*n0 + s*n1, -s.conj()*n0 + c*n1]
# get and apply current rotation to h and S
c, s, mag = lartg(h[col, col], h[col, col+1])
givens[col, :] = [c, s]
h[col, [col, col+1]] = mag, 0
# S[col+1] component is always 0
tmp = -np.conjugate(s)*S[col]
S[[col, col + 1]] = [c*S[col], tmp]
presid = np.abs(tmp)
inner_iter += 1
if callback_type in ('legacy', 'pr_norm'):
callback(presid / bnrm2)
# Legacy behavior
if callback_type == 'legacy' and inner_iter == maxiter:
break
if presid <= ptol or breakdown:
break
# Solve h(col, col) upper triangular system and allow pseudo-solve
# singular cases as in (but without the f2py copies):
# y = trsv(h[:col+1, :col+1].T, S[:col+1])
if h[col, col] == 0:
S[col] = 0
y = np.zeros([col+1], dtype=x.dtype)
y[:] = S[:col+1]
for k in range(col, 0, -1):
if y[k] != 0:
y[k] /= h[k, k]
tmp = y[k]
y[:k] -= tmp*h[k, :k]
if y[0] != 0:
y[0] /= h[0, 0]
x += y @ v[:col+1, :]
r = b - matvec(x)
rnorm = np.linalg.norm(r)
# Legacy exit
if callback_type == 'legacy' and inner_iter == maxiter:
return postprocess(x), 0 if rnorm <= atol else maxiter
if callback_type == 'x':
callback(x)
if rnorm <= atol:
break
elif breakdown:
# Reached breakdown (= exact solution), but the external
# tolerance check failed. Bail out with failure.
break
elif presid <= ptol:
# Inner loop passed but outer didn't
ptol_max_factor = max(eps, 0.25 * ptol_max_factor)
else:
ptol_max_factor = min(1.0, 1.5 * ptol_max_factor)
ptol = presid * min(ptol_max_factor, atol / rnorm)
info = 0 if (rnorm <= atol) else maxiter
return postprocess(x), info
def qmr(A, b, x0=None, *, rtol=1e-5, atol=0., maxiter=None, M1=None, M2=None,
callback=None):
"""Use Quasi-Minimal Residual iteration to solve ``Ax = b``.
Parameters
----------
A : {sparse array, ndarray, LinearOperator}
The real-valued N-by-N matrix of the linear system.
Alternatively, ``A`` can be a linear operator which can
produce ``Ax`` and ``A^T x`` using, e.g.,
``scipy.sparse.linalg.LinearOperator``.
b : ndarray
Right hand side of the linear system. Has shape (N,) or (N,1).
x0 : ndarray
Starting guess for the solution.
atol, rtol : float, optional
Parameters for the convergence test. For convergence,
``norm(b - A @ x) <= max(rtol*norm(b), atol)`` should be satisfied.
The default is ``atol=0.`` and ``rtol=1e-5``.
maxiter : integer
Maximum number of iterations. Iteration will stop after maxiter
steps even if the specified tolerance has not been achieved.
M1 : {sparse array, ndarray, LinearOperator}
Left preconditioner for A.
M2 : {sparse array, ndarray, LinearOperator}
Right preconditioner for A. Used together with the left
preconditioner M1. The matrix M1@A@M2 should have better
conditioned than A alone.
callback : function
User-supplied function to call after each iteration. It is called
as callback(xk), where xk is the current solution vector.
Returns
-------
x : ndarray
The converged solution.
info : integer
Provides convergence information:
0 : successful exit
>0 : convergence to tolerance not achieved, number of iterations
<0 : parameter breakdown
See Also
--------
LinearOperator
Examples
--------
>>> import numpy as np
>>> from scipy.sparse import csc_array
>>> from scipy.sparse.linalg import qmr
>>> A = csc_array([[3., 2., 0.], [1., -1., 0.], [0., 5., 1.]])
>>> b = np.array([2., 4., -1.])
>>> x, exitCode = qmr(A, b, atol=1e-5)
>>> print(exitCode) # 0 indicates successful convergence
0
>>> np.allclose(A.dot(x), b)
True
"""
A_ = A
A, M, x, b, postprocess = make_system(A, None, x0, b)
bnrm2 = np.linalg.norm(b)
atol, _ = _get_atol_rtol('qmr', bnrm2, atol, rtol)
if bnrm2 == 0:
return postprocess(b), 0
if M1 is None and M2 is None:
if hasattr(A_, 'psolve'):
def left_psolve(b):
return A_.psolve(b, 'left')
def right_psolve(b):
return A_.psolve(b, 'right')
def left_rpsolve(b):
return A_.rpsolve(b, 'left')
def right_rpsolve(b):
return A_.rpsolve(b, 'right')
M1 = LinearOperator(A.shape,
matvec=left_psolve,
rmatvec=left_rpsolve)
M2 = LinearOperator(A.shape,
matvec=right_psolve,
rmatvec=right_rpsolve)
else:
def id(b):
return b
M1 = LinearOperator(A.shape, matvec=id, rmatvec=id)
M2 = LinearOperator(A.shape, matvec=id, rmatvec=id)
n = len(b)
if maxiter is None:
maxiter = n*10
dotprod = np.vdot if np.iscomplexobj(x) else np.dot
rhotol = np.finfo(x.dtype.char).eps
betatol = rhotol
gammatol = rhotol
deltatol = rhotol
epsilontol = rhotol
xitol = rhotol
r = b - A.matvec(x) if x.any() else b.copy()
vtilde = r.copy()
y = M1.matvec(vtilde)
rho = np.linalg.norm(y)
wtilde = r.copy()
z = M2.rmatvec(wtilde)
xi = np.linalg.norm(z)
gamma, eta, theta = 1, -1, 0
v = np.empty_like(vtilde)
w = np.empty_like(wtilde)
# Dummy values to initialize vars, silence linter warnings
epsilon, q, d, p, s = None, None, None, None, None
for iteration in range(maxiter):
if np.linalg.norm(r) < atol: # Are we done?
return postprocess(x), 0
if np.abs(rho) < rhotol: # rho breakdown
return postprocess(x), -10
if np.abs(xi) < xitol: # xi breakdown
return postprocess(x), -15
v[:] = vtilde[:]
v *= (1 / rho)
y *= (1 / rho)
w[:] = wtilde[:]
w *= (1 / xi)
z *= (1 / xi)
delta = dotprod(z, y)
if np.abs(delta) < deltatol: # delta breakdown
return postprocess(x), -13
ytilde = M2.matvec(y)
ztilde = M1.rmatvec(z)
if iteration > 0:
ytilde -= (xi * delta / epsilon) * p
p[:] = ytilde[:]
ztilde -= (rho * (delta / epsilon).conj()) * q
q[:] = ztilde[:]
else: # First spin
p = ytilde.copy()
q = ztilde.copy()
ptilde = A.matvec(p)
epsilon = dotprod(q, ptilde)
if np.abs(epsilon) < epsilontol: # epsilon breakdown
return postprocess(x), -14
beta = epsilon / delta
if np.abs(beta) < betatol: # beta breakdown
return postprocess(x), -11
vtilde[:] = ptilde[:]
vtilde -= beta*v
y = M1.matvec(vtilde)
rho_prev = rho
rho = np.linalg.norm(y)
wtilde[:] = w[:]
wtilde *= - beta.conj()
wtilde += A.rmatvec(q)
z = M2.rmatvec(wtilde)
xi = np.linalg.norm(z)
gamma_prev = gamma
theta_prev = theta
theta = rho / (gamma_prev * np.abs(beta))
gamma = 1 / np.sqrt(1 + theta**2)
if np.abs(gamma) < gammatol: # gamma breakdown
return postprocess(x), -12
eta *= -(rho_prev / beta) * (gamma / gamma_prev)**2
if iteration > 0:
d *= (theta_prev * gamma) ** 2
d += eta*p
s *= (theta_prev * gamma) ** 2
s += eta*ptilde
else:
d = p.copy()
d *= eta
s = ptilde.copy()
s *= eta
x += d
r -= s
if callback:
callback(x)
else: # for loop exhausted
# Return incomplete progress
return postprocess(x), maxiter