scipy.linalg.solve_circulant#
- scipy.linalg.solve_circulant(c, b, singular='raise', tol=None, caxis=-1, baxis=0, outaxis=0)[source]#
Solve C x = b for x, where C is a circulant matrix.
C is the circulant matrix associated with the vector c.
The system is solved by doing division in Fourier space. The calculation is:
x = ifft(fft(b) / fft(c))
where fft and ifft are the fast Fourier transform and its inverse, respectively. For a large vector c, this is much faster than solving the system with the full circulant matrix.
- Parameters:
- carray_like
The coefficients of the circulant matrix.
- barray_like
Right-hand side matrix in
a x = b
.- singularstr, optional
This argument controls how a near singular circulant matrix is handled. If singular is “raise” and the circulant matrix is near singular, a
LinAlgError
is raised. If singular is “lstsq”, the least squares solution is returned. Default is “raise”.- tolfloat, optional
If any eigenvalue of the circulant matrix has an absolute value that is less than or equal to tol, the matrix is considered to be near singular. If not given, tol is set to:
tol = abs_eigs.max() * abs_eigs.size * np.finfo(np.float64).eps
where abs_eigs is the array of absolute values of the eigenvalues of the circulant matrix.
- caxisint
When c has dimension greater than 1, it is viewed as a collection of circulant vectors. In this case, caxis is the axis of c that holds the vectors of circulant coefficients.
- baxisint
When b has dimension greater than 1, it is viewed as a collection of vectors. In this case, baxis is the axis of b that holds the right-hand side vectors.
- outaxisint
When c or b are multidimensional, the value returned by
solve_circulant
is multidimensional. In this case, outaxis is the axis of the result that holds the solution vectors.
- Returns:
- xndarray
Solution to the system
C x = b
.
- Raises:
- LinAlgError
If the circulant matrix associated with c is near singular.
See also
circulant
circulant matrix
Notes
For a 1-D vector c with length m, and an array b with shape
(m, ...)
,solve_circulant(c, b)
returns the same result as
solve(circulant(c), b)
where
solve
andcirculant
are fromscipy.linalg
.Added in version 0.16.0.
Examples
>>> import numpy as np >>> from scipy.linalg import solve_circulant, solve, circulant, lstsq
>>> c = np.array([2, 2, 4]) >>> b = np.array([1, 2, 3]) >>> solve_circulant(c, b) array([ 0.75, -0.25, 0.25])
Compare that result to solving the system with
scipy.linalg.solve
:>>> solve(circulant(c), b) array([ 0.75, -0.25, 0.25])
A singular example:
>>> c = np.array([1, 1, 0, 0]) >>> b = np.array([1, 2, 3, 4])
Calling
solve_circulant(c, b)
will raise aLinAlgError
. For the least square solution, use the optionsingular='lstsq'
:>>> solve_circulant(c, b, singular='lstsq') array([ 0.25, 1.25, 2.25, 1.25])
Compare to
scipy.linalg.lstsq
:>>> x, resid, rnk, s = lstsq(circulant(c), b) >>> x array([ 0.25, 1.25, 2.25, 1.25])
A broadcasting example:
Suppose we have the vectors of two circulant matrices stored in an array with shape (2, 5), and three b vectors stored in an array with shape (3, 5). For example,
>>> c = np.array([[1.5, 2, 3, 0, 0], [1, 1, 4, 3, 2]]) >>> b = np.arange(15).reshape(-1, 5)
We want to solve all combinations of circulant matrices and b vectors, with the result stored in an array with shape (2, 3, 5). When we disregard the axes of c and b that hold the vectors of coefficients, the shapes of the collections are (2,) and (3,), respectively, which are not compatible for broadcasting. To have a broadcast result with shape (2, 3), we add a trivial dimension to c:
c[:, np.newaxis, :]
has shape (2, 1, 5). The last dimension holds the coefficients of the circulant matrices, so when we callsolve_circulant
, we can use the defaultcaxis=-1
. The coefficients of the b vectors are in the last dimension of the array b, so we usebaxis=-1
. If we use the default outaxis, the result will have shape (5, 2, 3), so we’ll useoutaxis=-1
to put the solution vectors in the last dimension.>>> x = solve_circulant(c[:, np.newaxis, :], b, baxis=-1, outaxis=-1) >>> x.shape (2, 3, 5) >>> np.set_printoptions(precision=3) # For compact output of numbers. >>> x array([[[-0.118, 0.22 , 1.277, -0.142, 0.302], [ 0.651, 0.989, 2.046, 0.627, 1.072], [ 1.42 , 1.758, 2.816, 1.396, 1.841]], [[ 0.401, 0.304, 0.694, -0.867, 0.377], [ 0.856, 0.758, 1.149, -0.412, 0.831], [ 1.31 , 1.213, 1.603, 0.042, 1.286]]])
Check by solving one pair of c and b vectors (cf.
x[1, 1, :]
):>>> solve_circulant(c[1], b[1, :]) array([ 0.856, 0.758, 1.149, -0.412, 0.831])