Coverage for /usr/lib/python3/dist-packages/scipy/linalg/_decomp_lu.py: 13%
91 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""LU decomposition functions."""
3from warnings import warn
5from numpy import asarray, asarray_chkfinite
6import numpy as np
7from itertools import product
9# Local imports
10from ._misc import _datacopied, LinAlgWarning
11from .lapack import get_lapack_funcs
12from ._decomp_lu_cython import lu_dispatcher
14lapack_cast_dict = {x: ''.join([y for y in 'fdFD' if np.can_cast(x, y)])
15 for x in np.typecodes['All']}
17__all__ = ['lu', 'lu_solve', 'lu_factor']
20def lu_factor(a, overwrite_a=False, check_finite=True):
21 """
22 Compute pivoted LU decomposition of a matrix.
24 The decomposition is::
26 A = P L U
28 where P is a permutation matrix, L lower triangular with unit
29 diagonal elements, and U upper triangular.
31 Parameters
32 ----------
33 a : (M, N) array_like
34 Matrix to decompose
35 overwrite_a : bool, optional
36 Whether to overwrite data in A (may increase performance)
37 check_finite : bool, optional
38 Whether to check that the input matrix contains only finite numbers.
39 Disabling may give a performance gain, but may result in problems
40 (crashes, non-termination) if the inputs do contain infinities or NaNs.
42 Returns
43 -------
44 lu : (M, N) ndarray
45 Matrix containing U in its upper triangle, and L in its lower triangle.
46 The unit diagonal elements of L are not stored.
47 piv : (N,) ndarray
48 Pivot indices representing the permutation matrix P:
49 row i of matrix was interchanged with row piv[i].
51 See Also
52 --------
53 lu : gives lu factorization in more user-friendly format
54 lu_solve : solve an equation system using the LU factorization of a matrix
56 Notes
57 -----
58 This is a wrapper to the ``*GETRF`` routines from LAPACK. Unlike
59 :func:`lu`, it outputs the L and U factors into a single array
60 and returns pivot indices instead of a permutation matrix.
62 Examples
63 --------
64 >>> import numpy as np
65 >>> from scipy.linalg import lu_factor
66 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
67 >>> lu, piv = lu_factor(A)
68 >>> piv
69 array([2, 2, 3, 3], dtype=int32)
71 Convert LAPACK's ``piv`` array to NumPy index and test the permutation
73 >>> piv_py = [2, 0, 3, 1]
74 >>> L, U = np.tril(lu, k=-1) + np.eye(4), np.triu(lu)
75 >>> np.allclose(A[piv_py] - L @ U, np.zeros((4, 4)))
76 True
77 """
78 if check_finite:
79 a1 = asarray_chkfinite(a)
80 else:
81 a1 = asarray(a)
82 overwrite_a = overwrite_a or (_datacopied(a1, a))
83 getrf, = get_lapack_funcs(('getrf',), (a1,))
84 lu, piv, info = getrf(a1, overwrite_a=overwrite_a)
85 if info < 0:
86 raise ValueError('illegal value in %dth argument of '
87 'internal getrf (lu_factor)' % -info)
88 if info > 0:
89 warn("Diagonal number %d is exactly zero. Singular matrix." % info,
90 LinAlgWarning, stacklevel=2)
91 return lu, piv
94def lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True):
95 """Solve an equation system, a x = b, given the LU factorization of a
97 Parameters
98 ----------
99 (lu, piv)
100 Factorization of the coefficient matrix a, as given by lu_factor
101 b : array
102 Right-hand side
103 trans : {0, 1, 2}, optional
104 Type of system to solve:
106 ===== =========
107 trans system
108 ===== =========
109 0 a x = b
110 1 a^T x = b
111 2 a^H x = b
112 ===== =========
113 overwrite_b : bool, optional
114 Whether to overwrite data in b (may increase performance)
115 check_finite : bool, optional
116 Whether to check that the input matrices contain only finite numbers.
117 Disabling may give a performance gain, but may result in problems
118 (crashes, non-termination) if the inputs do contain infinities or NaNs.
120 Returns
121 -------
122 x : array
123 Solution to the system
125 See Also
126 --------
127 lu_factor : LU factorize a matrix
129 Examples
130 --------
131 >>> import numpy as np
132 >>> from scipy.linalg import lu_factor, lu_solve
133 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
134 >>> b = np.array([1, 1, 1, 1])
135 >>> lu, piv = lu_factor(A)
136 >>> x = lu_solve((lu, piv), b)
137 >>> np.allclose(A @ x - b, np.zeros((4,)))
138 True
140 """
141 (lu, piv) = lu_and_piv
142 if check_finite:
143 b1 = asarray_chkfinite(b)
144 else:
145 b1 = asarray(b)
146 overwrite_b = overwrite_b or _datacopied(b1, b)
147 if lu.shape[0] != b1.shape[0]:
148 raise ValueError("Shapes of lu {} and b {} are incompatible"
149 .format(lu.shape, b1.shape))
151 getrs, = get_lapack_funcs(('getrs',), (lu, b1))
152 x, info = getrs(lu, piv, b1, trans=trans, overwrite_b=overwrite_b)
153 if info == 0:
154 return x
155 raise ValueError('illegal value in %dth argument of internal gesv|posv'
156 % -info)
159def lu(a, permute_l=False, overwrite_a=False, check_finite=True,
160 p_indices=False):
161 """
162 Compute LU decomposition of a matrix with partial pivoting.
164 The decomposition satisfies::
166 A = P @ L @ U
168 where ``P`` is a permutation matrix, ``L`` lower triangular with unit
169 diagonal elements, and ``U`` upper triangular. If `permute_l` is set to
170 ``True`` then ``L`` is returned already permuted and hence satisfying
171 ``A = L @ U``.
173 Parameters
174 ----------
175 a : (M, N) array_like
176 Array to decompose
177 permute_l : bool, optional
178 Perform the multiplication P*L (Default: do not permute)
179 overwrite_a : bool, optional
180 Whether to overwrite data in a (may improve performance)
181 check_finite : bool, optional
182 Whether to check that the input matrix contains only finite numbers.
183 Disabling may give a performance gain, but may result in problems
184 (crashes, non-termination) if the inputs do contain infinities or NaNs.
185 p_indices : bool, optional
186 If ``True`` the permutation information is returned as row indices.
187 The default is ``False`` for backwards-compatibility reasons.
189 Returns
190 -------
191 **(If `permute_l` is ``False``)**
193 p : (..., M, M) ndarray
194 Permutation arrays or vectors depending on `p_indices`
195 l : (..., M, K) ndarray
196 Lower triangular or trapezoidal array with unit diagonal.
197 ``K = min(M, N)``
198 u : (..., K, N) ndarray
199 Upper triangular or trapezoidal array
201 **(If `permute_l` is ``True``)**
203 pl : (..., M, K) ndarray
204 Permuted L matrix.
205 ``K = min(M, N)``
206 u : (..., K, N) ndarray
207 Upper triangular or trapezoidal array
209 Notes
210 -----
211 Permutation matrices are costly since they are nothing but row reorder of
212 ``L`` and hence indices are strongly recommended to be used instead if the
213 permutation is required. The relation in the 2D case then becomes simply
214 ``A = L[P, :] @ U``. In higher dimensions, it is better to use `permute_l`
215 to avoid complicated indexing tricks.
217 In 2D case, if one has the indices however, for some reason, the
218 permutation matrix is still needed then it can be constructed by
219 ``np.eye(M)[P, :]``.
221 Examples
222 --------
224 >>> import numpy as np
225 >>> from scipy.linalg import lu
226 >>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
227 >>> p, l, u = lu(A)
228 >>> np.allclose(A, p @ l @ u)
229 True
230 >>> p # Permutation matrix
231 array([[0., 1., 0., 0.], # Row index 1
232 [0., 0., 0., 1.], # Row index 3
233 [1., 0., 0., 0.], # Row index 0
234 [0., 0., 1., 0.]]) # Row index 2
235 >>> p, _, _ = lu(A, p_indices=True)
236 >>> p
237 array([1, 3, 0, 2]) # as given by row indices above
238 >>> np.allclose(A, l[p, :] @ u)
239 True
241 We can also use nd-arrays, for example, a demonstration with 4D array:
243 >>> rng = np.random.default_rng()
244 >>> A = rng.uniform(low=-4, high=4, size=[3, 2, 4, 8])
245 >>> p, l, u = lu(A)
246 >>> p.shape, l.shape, u.shape
247 ((3, 2, 4, 4), (3, 2, 4, 4), (3, 2, 4, 8))
248 >>> np.allclose(A, p @ l @ u)
249 True
250 >>> PL, U = lu(A, permute_l=True)
251 >>> np.allclose(A, PL @ U)
252 True
254 """
255 a1 = np.asarray_chkfinite(a) if check_finite else np.asarray(a)
256 if a1.ndim < 2:
257 raise ValueError('The input array must be at least two-dimensional.')
259 # Also check if dtype is LAPACK compatible
260 if a1.dtype.char not in 'fdFD':
261 dtype_char = lapack_cast_dict[a1.dtype.char]
262 if not dtype_char: # No casting possible
263 raise TypeError(f'The dtype {a1.dtype} cannot be cast '
264 'to float(32, 64) or complex(64, 128).')
266 a1 = a1.astype(dtype_char[0]) # makes a copy, free to scratch
267 overwrite_a = True
269 *nd, m, n = a1.shape
270 k = min(m, n)
271 real_dchar = 'f' if a1.dtype.char in 'fF' else 'd'
273 # Empty input
274 if min(*a1.shape) == 0:
275 if permute_l:
276 PL = np.empty(shape=[*nd, m, k], dtype=a1.dtype)
277 U = np.empty(shape=[*nd, k, n], dtype=a1.dtype)
278 return PL, U
279 else:
280 P = (np.empty([*nd, 0], dtype=np.int32) if p_indices else
281 np.empty([*nd, 0, 0], dtype=real_dchar))
282 L = np.empty(shape=[*nd, m, k], dtype=a1.dtype)
283 U = np.empty(shape=[*nd, k, n], dtype=a1.dtype)
284 return P, L, U
286 # Scalar case
287 if a1.shape[-2:] == (1, 1):
288 if permute_l:
289 return np.ones_like(a1), (a1 if overwrite_a else a1.copy())
290 else:
291 P = (np.zeros(shape=[*nd, m], dtype=int) if p_indices
292 else np.ones_like(a1))
293 return P, np.ones_like(a1), (a1 if overwrite_a else a1.copy())
295 # Then check overwrite permission
296 if not _datacopied(a1, a): # "a" still alive through "a1"
297 if not overwrite_a:
298 # Data belongs to "a" so make a copy
299 a1 = a1.copy(order='C')
300 # else: Do nothing we'll use "a" if possible
301 # else: a1 has its own data thus free to scratch
303 # Then layout checks, might happen that overwrite is allowed but original
304 # array was read-only or non-contiguous.
306 if not (a1.flags['C_CONTIGUOUS'] and a1.flags['WRITEABLE']):
307 a1 = a1.copy(order='C')
309 if not nd: # 2D array
311 p = np.empty(m, dtype=np.int32)
312 u = np.zeros([k, k], dtype=a1.dtype)
313 lu_dispatcher(a1, u, p, permute_l)
314 P, L, U = (p, a1, u) if m > n else (p, u, a1)
316 else: # Stacked array
318 # Prepare the contiguous data holders
319 P = np.empty([*nd, m], dtype=np.int32) # perm vecs
321 if m > n: # Tall arrays, U will be created
322 U = np.zeros([*nd, k, k], dtype=a1.dtype)
323 for ind in product(*[range(x) for x in a1.shape[:-2]]):
324 lu_dispatcher(a1[ind], U[ind], P[ind], permute_l)
325 L = a1
327 else: # Fat arrays, L will be created
328 L = np.zeros([*nd, k, k], dtype=a1.dtype)
329 for ind in product(*[range(x) for x in a1.shape[:-2]]):
330 lu_dispatcher(a1[ind], L[ind], P[ind], permute_l)
331 U = a1
333 # Convert permutation vecs to permutation arrays
334 # permute_l=False needed to enter here to avoid wasted efforts
335 if (not p_indices) and (not permute_l):
336 if nd:
337 Pa = np.zeros([*nd, m, m], dtype=real_dchar)
338 # An unreadable index hack - One-hot encoding for perm matrices
339 nd_ix = np.ix_(*([np.arange(x) for x in nd]+[np.arange(m)]))
340 Pa[(*nd_ix, P)] = 1
341 P = Pa
342 else: # 2D case
343 Pa = np.zeros([m, m], dtype=real_dchar)
344 Pa[np.arange(m), P] = 1
345 P = Pa
347 return (L, U) if permute_l else (P, L, U)