Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp_ldl.py: 12%
84 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1from warnings import warn
3import numpy as np
4from numpy import (atleast_2d, ComplexWarning, arange, zeros_like, imag, diag,
5 iscomplexobj, tril, triu, argsort, empty_like)
6from ._decomp import _asarray_validated
7from .lapack import get_lapack_funcs, _compute_lwork
9__all__ = ['ldl']
12def ldl(A, lower=True, hermitian=True, overwrite_a=False, check_finite=True):
13 """ Computes the LDLt or Bunch-Kaufman factorization of a symmetric/
14 hermitian matrix.
16 This function returns a block diagonal matrix D consisting blocks of size
17 at most 2x2 and also a possibly permuted unit lower triangular matrix
18 ``L`` such that the factorization ``A = L D L^H`` or ``A = L D L^T``
19 holds. If `lower` is False then (again possibly permuted) upper
20 triangular matrices are returned as outer factors.
22 The permutation array can be used to triangularize the outer factors
23 simply by a row shuffle, i.e., ``lu[perm, :]`` is an upper/lower
24 triangular matrix. This is also equivalent to multiplication with a
25 permutation matrix ``P.dot(lu)``, where ``P`` is a column-permuted
26 identity matrix ``I[:, perm]``.
28 Depending on the value of the boolean `lower`, only upper or lower
29 triangular part of the input array is referenced. Hence, a triangular
30 matrix on entry would give the same result as if the full matrix is
31 supplied.
33 Parameters
34 ----------
35 A : array_like
36 Square input array
37 lower : bool, optional
38 This switches between the lower and upper triangular outer factors of
39 the factorization. Lower triangular (``lower=True``) is the default.
40 hermitian : bool, optional
41 For complex-valued arrays, this defines whether ``A = A.conj().T`` or
42 ``A = A.T`` is assumed. For real-valued arrays, this switch has no
43 effect.
44 overwrite_a : bool, optional
45 Allow overwriting data in `A` (may enhance performance). The default
46 is False.
47 check_finite : bool, optional
48 Whether to check that the input matrices contain only finite numbers.
49 Disabling may give a performance gain, but may result in problems
50 (crashes, non-termination) if the inputs do contain infinities or NaNs.
52 Returns
53 -------
54 lu : ndarray
55 The (possibly) permuted upper/lower triangular outer factor of the
56 factorization.
57 d : ndarray
58 The block diagonal multiplier of the factorization.
59 perm : ndarray
60 The row-permutation index array that brings lu into triangular form.
62 Raises
63 ------
64 ValueError
65 If input array is not square.
66 ComplexWarning
67 If a complex-valued array with nonzero imaginary parts on the
68 diagonal is given and hermitian is set to True.
70 See Also
71 --------
72 cholesky, lu
74 Notes
75 -----
76 This function uses ``?SYTRF`` routines for symmetric matrices and
77 ``?HETRF`` routines for Hermitian matrices from LAPACK. See [1]_ for
78 the algorithm details.
80 Depending on the `lower` keyword value, only lower or upper triangular
81 part of the input array is referenced. Moreover, this keyword also defines
82 the structure of the outer factors of the factorization.
84 .. versionadded:: 1.1.0
86 References
87 ----------
88 .. [1] J.R. Bunch, L. Kaufman, Some stable methods for calculating
89 inertia and solving symmetric linear systems, Math. Comput. Vol.31,
90 1977. :doi:`10.2307/2005787`
92 Examples
93 --------
94 Given an upper triangular array ``a`` that represents the full symmetric
95 array with its entries, obtain ``l``, 'd' and the permutation vector `perm`:
97 >>> import numpy as np
98 >>> from scipy.linalg import ldl
99 >>> a = np.array([[2, -1, 3], [0, 2, 0], [0, 0, 1]])
100 >>> lu, d, perm = ldl(a, lower=0) # Use the upper part
101 >>> lu
102 array([[ 0. , 0. , 1. ],
103 [ 0. , 1. , -0.5],
104 [ 1. , 1. , 1.5]])
105 >>> d
106 array([[-5. , 0. , 0. ],
107 [ 0. , 1.5, 0. ],
108 [ 0. , 0. , 2. ]])
109 >>> perm
110 array([2, 1, 0])
111 >>> lu[perm, :]
112 array([[ 1. , 1. , 1.5],
113 [ 0. , 1. , -0.5],
114 [ 0. , 0. , 1. ]])
115 >>> lu.dot(d).dot(lu.T)
116 array([[ 2., -1., 3.],
117 [-1., 2., 0.],
118 [ 3., 0., 1.]])
120 """
121 a = atleast_2d(_asarray_validated(A, check_finite=check_finite))
122 if a.shape[0] != a.shape[1]:
123 raise ValueError('The input array "a" should be square.')
124 # Return empty arrays for empty square input
125 if a.size == 0:
126 return empty_like(a), empty_like(a), np.array([], dtype=int)
128 n = a.shape[0]
129 r_or_c = complex if iscomplexobj(a) else float
131 # Get the LAPACK routine
132 if r_or_c is complex and hermitian:
133 s, sl = 'hetrf', 'hetrf_lwork'
134 if np.any(imag(diag(a))):
135 warn('scipy.linalg.ldl():\nThe imaginary parts of the diagonal'
136 'are ignored. Use "hermitian=False" for factorization of'
137 'complex symmetric arrays.', ComplexWarning, stacklevel=2)
138 else:
139 s, sl = 'sytrf', 'sytrf_lwork'
141 solver, solver_lwork = get_lapack_funcs((s, sl), (a,))
142 lwork = _compute_lwork(solver_lwork, n, lower=lower)
143 ldu, piv, info = solver(a, lwork=lwork, lower=lower,
144 overwrite_a=overwrite_a)
145 if info < 0:
146 raise ValueError('{} exited with the internal error "illegal value '
147 'in argument number {}". See LAPACK documentation '
148 'for the error codes.'.format(s.upper(), -info))
150 swap_arr, pivot_arr = _ldl_sanitize_ipiv(piv, lower=lower)
151 d, lu = _ldl_get_d_and_l(ldu, pivot_arr, lower=lower, hermitian=hermitian)
152 lu, perm = _ldl_construct_tri_factor(lu, swap_arr, pivot_arr, lower=lower)
154 return lu, d, perm
157def _ldl_sanitize_ipiv(a, lower=True):
158 """
159 This helper function takes the rather strangely encoded permutation array
160 returned by the LAPACK routines ?(HE/SY)TRF and converts it into
161 regularized permutation and diagonal pivot size format.
163 Since FORTRAN uses 1-indexing and LAPACK uses different start points for
164 upper and lower formats there are certain offsets in the indices used
165 below.
167 Let's assume a result where the matrix is 6x6 and there are two 2x2
168 and two 1x1 blocks reported by the routine. To ease the coding efforts,
169 we still populate a 6-sized array and fill zeros as the following ::
171 pivots = [2, 0, 2, 0, 1, 1]
173 This denotes a diagonal matrix of the form ::
175 [x x ]
176 [x x ]
177 [ x x ]
178 [ x x ]
179 [ x ]
180 [ x]
182 In other words, we write 2 when the 2x2 block is first encountered and
183 automatically write 0 to the next entry and skip the next spin of the
184 loop. Thus, a separate counter or array appends to keep track of block
185 sizes are avoided. If needed, zeros can be filtered out later without
186 losing the block structure.
188 Parameters
189 ----------
190 a : ndarray
191 The permutation array ipiv returned by LAPACK
192 lower : bool, optional
193 The switch to select whether upper or lower triangle is chosen in
194 the LAPACK call.
196 Returns
197 -------
198 swap_ : ndarray
199 The array that defines the row/column swap operations. For example,
200 if row two is swapped with row four, the result is [0, 3, 2, 3].
201 pivots : ndarray
202 The array that defines the block diagonal structure as given above.
204 """
205 n = a.size
206 swap_ = arange(n)
207 pivots = zeros_like(swap_, dtype=int)
208 skip_2x2 = False
210 # Some upper/lower dependent offset values
211 # range (s)tart, r(e)nd, r(i)ncrement
212 x, y, rs, re, ri = (1, 0, 0, n, 1) if lower else (-1, -1, n-1, -1, -1)
214 for ind in range(rs, re, ri):
215 # If previous spin belonged already to a 2x2 block
216 if skip_2x2:
217 skip_2x2 = False
218 continue
220 cur_val = a[ind]
221 # do we have a 1x1 block or not?
222 if cur_val > 0:
223 if cur_val != ind+1:
224 # Index value != array value --> permutation required
225 swap_[ind] = swap_[cur_val-1]
226 pivots[ind] = 1
227 # Not.
228 elif cur_val < 0 and cur_val == a[ind+x]:
229 # first neg entry of 2x2 block identifier
230 if -cur_val != ind+2:
231 # Index value != array value --> permutation required
232 swap_[ind+x] = swap_[-cur_val-1]
233 pivots[ind+y] = 2
234 skip_2x2 = True
235 else: # Doesn't make sense, give up
236 raise ValueError('While parsing the permutation array '
237 'in "scipy.linalg.ldl", invalid entries '
238 'found. The array syntax is invalid.')
239 return swap_, pivots
242def _ldl_get_d_and_l(ldu, pivs, lower=True, hermitian=True):
243 """
244 Helper function to extract the diagonal and triangular matrices for
245 LDL.T factorization.
247 Parameters
248 ----------
249 ldu : ndarray
250 The compact output returned by the LAPACK routing
251 pivs : ndarray
252 The sanitized array of {0, 1, 2} denoting the sizes of the pivots. For
253 every 2 there is a succeeding 0.
254 lower : bool, optional
255 If set to False, upper triangular part is considered.
256 hermitian : bool, optional
257 If set to False a symmetric complex array is assumed.
259 Returns
260 -------
261 d : ndarray
262 The block diagonal matrix.
263 lu : ndarray
264 The upper/lower triangular matrix
265 """
266 is_c = iscomplexobj(ldu)
267 d = diag(diag(ldu))
268 n = d.shape[0]
269 blk_i = 0 # block index
271 # row/column offsets for selecting sub-, super-diagonal
272 x, y = (1, 0) if lower else (0, 1)
274 lu = tril(ldu, -1) if lower else triu(ldu, 1)
275 diag_inds = arange(n)
276 lu[diag_inds, diag_inds] = 1
278 for blk in pivs[pivs != 0]:
279 # increment the block index and check for 2s
280 # if 2 then copy the off diagonals depending on uplo
281 inc = blk_i + blk
283 if blk == 2:
284 d[blk_i+x, blk_i+y] = ldu[blk_i+x, blk_i+y]
285 # If Hermitian matrix is factorized, the cross-offdiagonal element
286 # should be conjugated.
287 if is_c and hermitian:
288 d[blk_i+y, blk_i+x] = ldu[blk_i+x, blk_i+y].conj()
289 else:
290 d[blk_i+y, blk_i+x] = ldu[blk_i+x, blk_i+y]
292 lu[blk_i+x, blk_i+y] = 0.
293 blk_i = inc
295 return d, lu
298def _ldl_construct_tri_factor(lu, swap_vec, pivs, lower=True):
299 """
300 Helper function to construct explicit outer factors of LDL factorization.
302 If lower is True the permuted factors are multiplied as L(1)*L(2)*...*L(k).
303 Otherwise, the permuted factors are multiplied as L(k)*...*L(2)*L(1). See
304 LAPACK documentation for more details.
306 Parameters
307 ----------
308 lu : ndarray
309 The triangular array that is extracted from LAPACK routine call with
310 ones on the diagonals.
311 swap_vec : ndarray
312 The array that defines the row swapping indices. If the kth entry is m
313 then rows k,m are swapped. Notice that the mth entry is not necessarily
314 k to avoid undoing the swapping.
315 pivs : ndarray
316 The array that defines the block diagonal structure returned by
317 _ldl_sanitize_ipiv().
318 lower : bool, optional
319 The boolean to switch between lower and upper triangular structure.
321 Returns
322 -------
323 lu : ndarray
324 The square outer factor which satisfies the L * D * L.T = A
325 perm : ndarray
326 The permutation vector that brings the lu to the triangular form
328 Notes
329 -----
330 Note that the original argument "lu" is overwritten.
332 """
333 n = lu.shape[0]
334 perm = arange(n)
335 # Setup the reading order of the permutation matrix for upper/lower
336 rs, re, ri = (n-1, -1, -1) if lower else (0, n, 1)
338 for ind in range(rs, re, ri):
339 s_ind = swap_vec[ind]
340 if s_ind != ind:
341 # Column start and end positions
342 col_s = ind if lower else 0
343 col_e = n if lower else ind+1
345 # If we stumble upon a 2x2 block include both cols in the perm.
346 if pivs[ind] == (0 if lower else 2):
347 col_s += -1 if lower else 0
348 col_e += 0 if lower else 1
349 lu[[s_ind, ind], col_s:col_e] = lu[[ind, s_ind], col_s:col_e]
350 perm[[s_ind, ind]] = perm[[ind, s_ind]]
352 return lu, argsort(perm)