|
"""Cholesky decomposition functions."""
# Local imports
'cho_solve_banded']
check_finite=True): """Common code for cholesky() and cho_factor()."""
a1 = asarray_chkfinite(a) if check_finite else asarray(a) a1 = atleast_2d(a1)
# Dimension check if a1.ndim != 2: raise ValueError('Input array needs to be 2 dimensional but received ' 'a {}d-array.'.format(a1.ndim)) # Squareness check if a1.shape[0] != a1.shape[1]: raise ValueError('Input array is expected to be square but has ' 'the shape: {}.'.format(a1.shape))
# Quick return for square empty array if a1.size == 0: return a1.copy(), lower
overwrite_a = overwrite_a or _datacopied(a1, a) potrf, = get_lapack_funcs(('potrf',), (a1,)) c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean) if info > 0: raise LinAlgError("%d-th leading minor of the array is not positive " "definite" % info) if info < 0: raise ValueError('LAPACK reported an illegal value in {}-th argument' 'on entry to "POTRF".'.format(-info)) return c, lower
""" Compute the Cholesky decomposition of a matrix.
Returns the Cholesky decomposition, :math:`A = L L^*` or :math:`A = U^* U` of a Hermitian positive-definite matrix A.
Parameters ---------- a : (M, M) array_like Matrix to be decomposed lower : bool, optional Whether to compute the upper or lower triangular Cholesky factorization. Default is upper-triangular. overwrite_a : bool, optional Whether to overwrite data in `a` (may improve performance). check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- c : (M, M) ndarray Upper- or lower-triangular Cholesky factor of `a`.
Raises ------ LinAlgError : if decomposition fails.
Examples -------- >>> from scipy.linalg import cholesky >>> a = np.array([[1,-2j],[2j,5]]) >>> L = cholesky(a, lower=True) >>> L array([[ 1.+0.j, 0.+0.j], [ 0.+2.j, 1.+0.j]]) >>> L @ L.T.conj() array([[ 1.+0.j, 0.-2.j], [ 0.+2.j, 5.+0.j]])
""" c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True, check_finite=check_finite) return c
""" Compute the Cholesky decomposition of a matrix, to use in cho_solve
Returns a matrix containing the Cholesky decomposition, ``A = L L*`` or ``A = U* U`` of a Hermitian positive-definite matrix `a`. The return value can be directly used as the first parameter to cho_solve.
.. warning:: The returned matrix also contains random data in the entries not used by the Cholesky decomposition. If you need to zero these entries, use the function `cholesky` instead.
Parameters ---------- a : (M, M) array_like Matrix to be decomposed lower : bool, optional Whether to compute the upper or lower triangular Cholesky factorization (Default: upper-triangular) overwrite_a : bool, optional Whether to overwrite data in a (may improve performance) check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- c : (M, M) ndarray Matrix whose upper or lower triangle contains the Cholesky factor of `a`. Other parts of the matrix contain random data. lower : bool Flag indicating whether the factor is in the lower or upper triangle
Raises ------ LinAlgError Raised if decomposition fails.
See also -------- cho_solve : Solve a linear set equations using the Cholesky factorization of a matrix.
Examples -------- >>> from scipy.linalg import cho_factor >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]]) >>> c, low = cho_factor(A) >>> c array([[3. , 1. , 0.33333333, 1.66666667], [3. , 2.44948974, 1.90515869, -0.27216553], [1. , 5. , 2.29330749, 0.8559528 ], [5. , 1. , 2. , 1.55418563]]) >>> np.allclose(np.triu(c).T @ np. triu(c) - A, np.zeros((4, 4))) True
""" c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False, check_finite=check_finite) return c, lower
"""Solve the linear equations A x = b, given the Cholesky factorization of A.
Parameters ---------- (c, lower) : tuple, (array, bool) Cholesky factorization of a, as given by cho_factor b : array Right-hand side overwrite_b : bool, optional Whether to overwrite data in b (may improve performance) check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- x : array The solution to the system A x = b
See also -------- cho_factor : Cholesky factorization of a matrix
Examples -------- >>> from scipy.linalg import cho_factor, cho_solve >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]]) >>> c, low = cho_factor(A) >>> x = cho_solve((c, low), [1, 1, 1, 1]) >>> np.allclose(A @ x - [1, 1, 1, 1], np.zeros(4)) True
""" (c, lower) = c_and_lower if check_finite: b1 = asarray_chkfinite(b) c = asarray_chkfinite(c) else: b1 = asarray(b) c = asarray(c) if c.ndim != 2 or c.shape[0] != c.shape[1]: raise ValueError("The factored matrix c is not square.") if c.shape[1] != b1.shape[0]: raise ValueError("incompatible dimensions.")
overwrite_b = overwrite_b or _datacopied(b1, b)
potrs, = get_lapack_funcs(('potrs',), (c, b1)) x, info = potrs(c, b1, lower=lower, overwrite_b=overwrite_b) if info != 0: raise ValueError('illegal value in %d-th argument of internal potrs' % -info) return x
""" Cholesky decompose a banded Hermitian positive-definite matrix
The matrix a is stored in ab either in lower diagonal or upper diagonal ordered form::
ab[u + i - j, j] == a[i,j] (if upper form; i <= j) ab[ i - j, j] == a[i,j] (if lower form; i >= j)
Example of ab (shape of a is (6,6), u=2)::
upper form: * * a02 a13 a24 a35 * a01 a12 a23 a34 a45 a00 a11 a22 a33 a44 a55
lower form: a00 a11 a22 a33 a44 a55 a10 a21 a32 a43 a54 * a20 a31 a42 a53 * *
Parameters ---------- ab : (u + 1, M) array_like Banded matrix overwrite_ab : bool, optional Discard data in ab (may enhance performance) lower : bool, optional Is the matrix in the lower form. (Default is upper form) check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- c : (u + 1, M) ndarray Cholesky factorization of a, in the same banded format as ab
Examples -------- >>> from scipy.linalg import cholesky_banded >>> from numpy import allclose, zeros, diag >>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]]) >>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1) >>> A = A + A.conj().T + np.diag(Ab[2, :]) >>> c = cholesky_banded(Ab) >>> C = np.diag(c[0, 2:], k=2) + np.diag(c[1, 1:], k=1) + np.diag(c[2, :]) >>> np.allclose(C.conj().T @ C - A, np.zeros((5, 5))) True
""" if check_finite: ab = asarray_chkfinite(ab) else: ab = asarray(ab)
pbtrf, = get_lapack_funcs(('pbtrf',), (ab,)) c, info = pbtrf(ab, lower=lower, overwrite_ab=overwrite_ab) if info > 0: raise LinAlgError("%d-th leading minor not positive definite" % info) if info < 0: raise ValueError('illegal value in %d-th argument of internal pbtrf' % -info) return c
""" Solve the linear equations ``A x = b``, given the Cholesky factorization of the banded hermitian ``A``.
Parameters ---------- (cb, lower) : tuple, (ndarray, bool) `cb` is the Cholesky factorization of A, as given by cholesky_banded. `lower` must be the same value that was given to cholesky_banded. b : array_like Right-hand side overwrite_b : bool, optional If True, the function will overwrite the values in `b`. check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
Returns ------- x : array The solution to the system A x = b
See also -------- cholesky_banded : Cholesky factorization of a banded matrix
Notes -----
.. versionadded:: 0.8.0
Examples -------- >>> from scipy.linalg import cholesky_banded, cho_solve_banded >>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]]) >>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1) >>> A = A + A.conj().T + np.diag(Ab[2, :]) >>> c = cholesky_banded(Ab) >>> x = cho_solve_banded((c, False), np.ones(5)) >>> np.allclose(A @ x - np.ones(5), np.zeros(5)) True
""" (cb, lower) = cb_and_lower if check_finite: cb = asarray_chkfinite(cb) b = asarray_chkfinite(b) else: cb = asarray(cb) b = asarray(b)
# Validate shapes. if cb.shape[-1] != b.shape[0]: raise ValueError("shapes of cb and b are not compatible.")
pbtrs, = get_lapack_funcs(('pbtrs',), (cb, b)) x, info = pbtrs(cb, b, lower=lower, overwrite_b=overwrite_b) if info > 0: raise LinAlgError("%d-th leading minor not positive definite" % info) if info < 0: raise ValueError('illegal value in %d-th argument of internal pbtrs' % -info) return x |