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