Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp_cholesky.py: 14%
71 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""Cholesky decomposition functions."""
3from numpy import asarray_chkfinite, asarray, atleast_2d
5# Local imports
6from ._misc import LinAlgError, _datacopied
7from .lapack import get_lapack_funcs
9__all__ = ['cholesky', 'cho_factor', 'cho_solve', 'cholesky_banded',
10 'cho_solve_banded']
13def _cholesky(a, lower=False, overwrite_a=False, clean=True,
14 check_finite=True):
15 """Common code for cholesky() and cho_factor()."""
17 a1 = asarray_chkfinite(a) if check_finite else asarray(a)
18 a1 = atleast_2d(a1)
20 # Dimension check
21 if a1.ndim != 2:
22 raise ValueError('Input array needs to be 2D but received '
23 'a {}d-array.'.format(a1.ndim))
24 # Squareness check
25 if a1.shape[0] != a1.shape[1]:
26 raise ValueError('Input array is expected to be square but has '
27 'the shape: {}.'.format(a1.shape))
29 # Quick return for square empty array
30 if a1.size == 0:
31 return a1.copy(), lower
33 overwrite_a = overwrite_a or _datacopied(a1, a)
34 potrf, = get_lapack_funcs(('potrf',), (a1,))
35 c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean)
36 if info > 0:
37 raise LinAlgError("%d-th leading minor of the array is not positive "
38 "definite" % info)
39 if info < 0:
40 raise ValueError('LAPACK reported an illegal value in {}-th argument'
41 'on entry to "POTRF".'.format(-info))
42 return c, lower
45def cholesky(a, lower=False, overwrite_a=False, check_finite=True):
46 """
47 Compute the Cholesky decomposition of a matrix.
49 Returns the Cholesky decomposition, :math:`A = L L^*` or
50 :math:`A = U^* U` of a Hermitian positive-definite matrix A.
52 Parameters
53 ----------
54 a : (M, M) array_like
55 Matrix to be decomposed
56 lower : bool, optional
57 Whether to compute the upper- or lower-triangular Cholesky
58 factorization. Default is upper-triangular.
59 overwrite_a : bool, optional
60 Whether to overwrite data in `a` (may improve performance).
61 check_finite : bool, optional
62 Whether to check that the input matrix contains only finite numbers.
63 Disabling may give a performance gain, but may result in problems
64 (crashes, non-termination) if the inputs do contain infinities or NaNs.
66 Returns
67 -------
68 c : (M, M) ndarray
69 Upper- or lower-triangular Cholesky factor of `a`.
71 Raises
72 ------
73 LinAlgError : if decomposition fails.
75 Examples
76 --------
77 >>> import numpy as np
78 >>> from scipy.linalg import cholesky
79 >>> a = np.array([[1,-2j],[2j,5]])
80 >>> L = cholesky(a, lower=True)
81 >>> L
82 array([[ 1.+0.j, 0.+0.j],
83 [ 0.+2.j, 1.+0.j]])
84 >>> L @ L.T.conj()
85 array([[ 1.+0.j, 0.-2.j],
86 [ 0.+2.j, 5.+0.j]])
88 """
89 c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True,
90 check_finite=check_finite)
91 return c
94def cho_factor(a, lower=False, overwrite_a=False, check_finite=True):
95 """
96 Compute the Cholesky decomposition of a matrix, to use in cho_solve
98 Returns a matrix containing the Cholesky decomposition,
99 ``A = L L*`` or ``A = U* U`` of a Hermitian positive-definite matrix `a`.
100 The return value can be directly used as the first parameter to cho_solve.
102 .. warning::
103 The returned matrix also contains random data in the entries not
104 used by the Cholesky decomposition. If you need to zero these
105 entries, use the function `cholesky` instead.
107 Parameters
108 ----------
109 a : (M, M) array_like
110 Matrix to be decomposed
111 lower : bool, optional
112 Whether to compute the upper or lower triangular Cholesky factorization
113 (Default: upper-triangular)
114 overwrite_a : bool, optional
115 Whether to overwrite data in a (may improve performance)
116 check_finite : bool, optional
117 Whether to check that the input matrix contains only finite numbers.
118 Disabling may give a performance gain, but may result in problems
119 (crashes, non-termination) if the inputs do contain infinities or NaNs.
121 Returns
122 -------
123 c : (M, M) ndarray
124 Matrix whose upper or lower triangle contains the Cholesky factor
125 of `a`. Other parts of the matrix contain random data.
126 lower : bool
127 Flag indicating whether the factor is in the lower or upper triangle
129 Raises
130 ------
131 LinAlgError
132 Raised if decomposition fails.
134 See Also
135 --------
136 cho_solve : Solve a linear set equations using the Cholesky factorization
137 of a matrix.
139 Examples
140 --------
141 >>> import numpy as np
142 >>> from scipy.linalg import cho_factor
143 >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
144 >>> c, low = cho_factor(A)
145 >>> c
146 array([[3. , 1. , 0.33333333, 1.66666667],
147 [3. , 2.44948974, 1.90515869, -0.27216553],
148 [1. , 5. , 2.29330749, 0.8559528 ],
149 [5. , 1. , 2. , 1.55418563]])
150 >>> np.allclose(np.triu(c).T @ np. triu(c) - A, np.zeros((4, 4)))
151 True
153 """
154 c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False,
155 check_finite=check_finite)
156 return c, lower
159def cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True):
160 """Solve the linear equations A x = b, given the Cholesky factorization of A.
162 Parameters
163 ----------
164 (c, lower) : tuple, (array, bool)
165 Cholesky factorization of a, as given by cho_factor
166 b : array
167 Right-hand side
168 overwrite_b : bool, optional
169 Whether to overwrite data in b (may improve performance)
170 check_finite : bool, optional
171 Whether to check that the input matrices contain only finite numbers.
172 Disabling may give a performance gain, but may result in problems
173 (crashes, non-termination) if the inputs do contain infinities or NaNs.
175 Returns
176 -------
177 x : array
178 The solution to the system A x = b
180 See Also
181 --------
182 cho_factor : Cholesky factorization of a matrix
184 Examples
185 --------
186 >>> import numpy as np
187 >>> from scipy.linalg import cho_factor, cho_solve
188 >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
189 >>> c, low = cho_factor(A)
190 >>> x = cho_solve((c, low), [1, 1, 1, 1])
191 >>> np.allclose(A @ x - [1, 1, 1, 1], np.zeros(4))
192 True
194 """
195 (c, lower) = c_and_lower
196 if check_finite:
197 b1 = asarray_chkfinite(b)
198 c = asarray_chkfinite(c)
199 else:
200 b1 = asarray(b)
201 c = asarray(c)
202 if c.ndim != 2 or c.shape[0] != c.shape[1]:
203 raise ValueError("The factored matrix c is not square.")
204 if c.shape[1] != b1.shape[0]:
205 raise ValueError("incompatible dimensions ({} and {})"
206 .format(c.shape, b1.shape))
208 overwrite_b = overwrite_b or _datacopied(b1, b)
210 potrs, = get_lapack_funcs(('potrs',), (c, b1))
211 x, info = potrs(c, b1, lower=lower, overwrite_b=overwrite_b)
212 if info != 0:
213 raise ValueError('illegal value in %dth argument of internal potrs'
214 % -info)
215 return x
218def cholesky_banded(ab, overwrite_ab=False, lower=False, check_finite=True):
219 """
220 Cholesky decompose a banded Hermitian positive-definite matrix
222 The matrix a is stored in ab either in lower-diagonal or upper-
223 diagonal ordered form::
225 ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
226 ab[ i - j, j] == a[i,j] (if lower form; i >= j)
228 Example of ab (shape of a is (6,6), u=2)::
230 upper form:
231 * * a02 a13 a24 a35
232 * a01 a12 a23 a34 a45
233 a00 a11 a22 a33 a44 a55
235 lower form:
236 a00 a11 a22 a33 a44 a55
237 a10 a21 a32 a43 a54 *
238 a20 a31 a42 a53 * *
240 Parameters
241 ----------
242 ab : (u + 1, M) array_like
243 Banded matrix
244 overwrite_ab : bool, optional
245 Discard data in ab (may enhance performance)
246 lower : bool, optional
247 Is the matrix in the lower form. (Default is upper form)
248 check_finite : bool, optional
249 Whether to check that the input matrix contains only finite numbers.
250 Disabling may give a performance gain, but may result in problems
251 (crashes, non-termination) if the inputs do contain infinities or NaNs.
253 Returns
254 -------
255 c : (u + 1, M) ndarray
256 Cholesky factorization of a, in the same banded format as ab
258 See Also
259 --------
260 cho_solve_banded :
261 Solve a linear set equations, given the Cholesky factorization
262 of a banded Hermitian.
264 Examples
265 --------
266 >>> import numpy as np
267 >>> from scipy.linalg import cholesky_banded
268 >>> from numpy import allclose, zeros, diag
269 >>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]])
270 >>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1)
271 >>> A = A + A.conj().T + np.diag(Ab[2, :])
272 >>> c = cholesky_banded(Ab)
273 >>> C = np.diag(c[0, 2:], k=2) + np.diag(c[1, 1:], k=1) + np.diag(c[2, :])
274 >>> np.allclose(C.conj().T @ C - A, np.zeros((5, 5)))
275 True
277 """
278 if check_finite:
279 ab = asarray_chkfinite(ab)
280 else:
281 ab = asarray(ab)
283 pbtrf, = get_lapack_funcs(('pbtrf',), (ab,))
284 c, info = pbtrf(ab, lower=lower, overwrite_ab=overwrite_ab)
285 if info > 0:
286 raise LinAlgError("%d-th leading minor not positive definite" % info)
287 if info < 0:
288 raise ValueError('illegal value in %d-th argument of internal pbtrf'
289 % -info)
290 return c
293def cho_solve_banded(cb_and_lower, b, overwrite_b=False, check_finite=True):
294 """
295 Solve the linear equations ``A x = b``, given the Cholesky factorization of
296 the banded Hermitian ``A``.
298 Parameters
299 ----------
300 (cb, lower) : tuple, (ndarray, bool)
301 `cb` is the Cholesky factorization of A, as given by cholesky_banded.
302 `lower` must be the same value that was given to cholesky_banded.
303 b : array_like
304 Right-hand side
305 overwrite_b : bool, optional
306 If True, the function will overwrite the values in `b`.
307 check_finite : bool, optional
308 Whether to check that the input matrices contain only finite numbers.
309 Disabling may give a performance gain, but may result in problems
310 (crashes, non-termination) if the inputs do contain infinities or NaNs.
312 Returns
313 -------
314 x : array
315 The solution to the system A x = b
317 See Also
318 --------
319 cholesky_banded : Cholesky factorization of a banded matrix
321 Notes
322 -----
324 .. versionadded:: 0.8.0
326 Examples
327 --------
328 >>> import numpy as np
329 >>> from scipy.linalg import cholesky_banded, cho_solve_banded
330 >>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]])
331 >>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1)
332 >>> A = A + A.conj().T + np.diag(Ab[2, :])
333 >>> c = cholesky_banded(Ab)
334 >>> x = cho_solve_banded((c, False), np.ones(5))
335 >>> np.allclose(A @ x - np.ones(5), np.zeros(5))
336 True
338 """
339 (cb, lower) = cb_and_lower
340 if check_finite:
341 cb = asarray_chkfinite(cb)
342 b = asarray_chkfinite(b)
343 else:
344 cb = asarray(cb)
345 b = asarray(b)
347 # Validate shapes.
348 if cb.shape[-1] != b.shape[0]:
349 raise ValueError("shapes of cb and b are not compatible.")
351 pbtrs, = get_lapack_funcs(('pbtrs',), (cb, b))
352 x, info = pbtrs(cb, b, lower=lower, overwrite_b=overwrite_b)
353 if info > 0:
354 raise LinAlgError("%dth leading minor not positive definite" % info)
355 if info < 0:
356 raise ValueError('illegal value in %dth argument of internal pbtrs'
357 % -info)
358 return x