Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_decomp_cholesky.py: 14%

1"""Cholesky decomposition functions.""" 

2 

3from numpy import asarray_chkfinite, asarray, atleast_2d 

4 

5# Local imports 

6from ._misc import LinAlgError, _datacopied 

7from .lapack import get_lapack_funcs 

8 

9__all__ = ['cholesky', 'cho_factor', 'cho_solve', 'cholesky_banded', 

10 'cho_solve_banded'] 

11 

12 

13def _cholesky(a, lower=False, overwrite_a=False, clean=True, 

14 check_finite=True): 

15 """Common code for cholesky() and cho_factor().""" 

16 

17 a1 = asarray_chkfinite(a) if check_finite else asarray(a) 

18 a1 = atleast_2d(a1) 

19 

20 # Dimension check 

21 if a1.ndim != 2: 

22 raise ValueError(f'Input array needs to be 2D but received a {a1.ndim}d-array.') 

23 # Squareness check 

24 if a1.shape[0] != a1.shape[1]: 

25 raise ValueError('Input array is expected to be square but has ' 

26 f'the shape: {a1.shape}.') 

27 

28 # Quick return for square empty array 

29 if a1.size == 0: 

30 return a1.copy(), lower 

31 

32 overwrite_a = overwrite_a or _datacopied(a1, a) 

33 potrf, = get_lapack_funcs(('potrf',), (a1,)) 

34 c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean) 

35 if info > 0: 

36 raise LinAlgError("%d-th leading minor of the array is not positive " 

37 "definite" % info) 

38 if info < 0: 

39 raise ValueError(f'LAPACK reported an illegal value in {-info}-th argument' 

40 'on entry to "POTRF".') 

41 return c, lower 

42 

43 

44def cholesky(a, lower=False, overwrite_a=False, check_finite=True): 

45 """ 

46 Compute the Cholesky decomposition of a matrix. 

47 

48 Returns the Cholesky decomposition, :math:`A = L L^*` or 

49 :math:`A = U^* U` of a Hermitian positive-definite matrix A. 

50 

51 Parameters 

52 ---------- 

53 a : (M, M) array_like 

54 Matrix to be decomposed 

55 lower : bool, optional 

56 Whether to compute the upper- or lower-triangular Cholesky 

57 factorization. Default is upper-triangular. 

58 overwrite_a : bool, optional 

59 Whether to overwrite data in `a` (may improve performance). 

60 check_finite : bool, optional 

61 Whether to check that the input matrix contains only finite numbers. 

62 Disabling may give a performance gain, but may result in problems 

63 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

64 

65 Returns 

66 ------- 

67 c : (M, M) ndarray 

68 Upper- or lower-triangular Cholesky factor of `a`. 

69 

70 Raises 

71 ------ 

72 LinAlgError : if decomposition fails. 

73 

74 Examples 

75 -------- 

76 >>> import numpy as np 

77 >>> from scipy.linalg import cholesky 

78 >>> a = np.array([[1,-2j],[2j,5]]) 

79 >>> L = cholesky(a, lower=True) 

80 >>> L 

81 array([[ 1.+0.j, 0.+0.j], 

82 [ 0.+2.j, 1.+0.j]]) 

83 >>> L @ L.T.conj() 

84 array([[ 1.+0.j, 0.-2.j], 

85 [ 0.+2.j, 5.+0.j]]) 

86 

87 """ 

88 c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True, 

89 check_finite=check_finite) 

90 return c 

91 

92 

93def cho_factor(a, lower=False, overwrite_a=False, check_finite=True): 

94 """ 

95 Compute the Cholesky decomposition of a matrix, to use in cho_solve 

96 

97 Returns a matrix containing the Cholesky decomposition, 

98 ``A = L L*`` or ``A = U* U`` of a Hermitian positive-definite matrix `a`. 

99 The return value can be directly used as the first parameter to cho_solve. 

100 

101 .. warning:: 

102 The returned matrix also contains random data in the entries not 

103 used by the Cholesky decomposition. If you need to zero these 

104 entries, use the function `cholesky` instead. 

105 

106 Parameters 

107 ---------- 

108 a : (M, M) array_like 

109 Matrix to be decomposed 

110 lower : bool, optional 

111 Whether to compute the upper or lower triangular Cholesky factorization 

112 (Default: upper-triangular) 

113 overwrite_a : bool, optional 

114 Whether to overwrite data in a (may improve performance) 

115 check_finite : bool, optional 

116 Whether to check that the input matrix contains 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. 

119 

120 Returns 

121 ------- 

122 c : (M, M) ndarray 

123 Matrix whose upper or lower triangle contains the Cholesky factor 

124 of `a`. Other parts of the matrix contain random data. 

125 lower : bool 

126 Flag indicating whether the factor is in the lower or upper triangle 

127 

128 Raises 

129 ------ 

130 LinAlgError 

131 Raised if decomposition fails. 

132 

133 See Also 

134 -------- 

135 cho_solve : Solve a linear set equations using the Cholesky factorization 

136 of a matrix. 

137 

138 Examples 

139 -------- 

140 >>> import numpy as np 

141 >>> from scipy.linalg import cho_factor 

142 >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]]) 

143 >>> c, low = cho_factor(A) 

144 >>> c 

145 array([[3. , 1. , 0.33333333, 1.66666667], 

146 [3. , 2.44948974, 1.90515869, -0.27216553], 

147 [1. , 5. , 2.29330749, 0.8559528 ], 

148 [5. , 1. , 2. , 1.55418563]]) 

149 >>> np.allclose(np.triu(c).T @ np. triu(c) - A, np.zeros((4, 4))) 

150 True 

151 

152 """ 

153 c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=False, 

154 check_finite=check_finite) 

155 return c, lower 

156 

157 

158def cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True): 

159 """Solve the linear equations A x = b, given the Cholesky factorization of A. 

160 

161 Parameters 

162 ---------- 

163 (c, lower) : tuple, (array, bool) 

164 Cholesky factorization of a, as given by cho_factor 

165 b : array 

166 Right-hand side 

167 overwrite_b : bool, optional 

168 Whether to overwrite data in b (may improve performance) 

169 check_finite : bool, optional 

170 Whether to check that the input matrices contain only finite numbers. 

171 Disabling may give a performance gain, but may result in problems 

172 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

173 

174 Returns 

175 ------- 

176 x : array 

177 The solution to the system A x = b 

178 

179 See Also 

180 -------- 

181 cho_factor : Cholesky factorization of a matrix 

182 

183 Examples 

184 -------- 

185 >>> import numpy as np 

186 >>> from scipy.linalg import cho_factor, cho_solve 

187 >>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]]) 

188 >>> c, low = cho_factor(A) 

189 >>> x = cho_solve((c, low), [1, 1, 1, 1]) 

190 >>> np.allclose(A @ x - [1, 1, 1, 1], np.zeros(4)) 

191 True 

192 

193 """ 

194 (c, lower) = c_and_lower 

195 if check_finite: 

196 b1 = asarray_chkfinite(b) 

197 c = asarray_chkfinite(c) 

198 else: 

199 b1 = asarray(b) 

200 c = asarray(c) 

201 if c.ndim != 2 or c.shape[0] != c.shape[1]: 

202 raise ValueError("The factored matrix c is not square.") 

203 if c.shape[1] != b1.shape[0]: 

204 raise ValueError(f"incompatible dimensions ({c.shape} and {b1.shape})") 

205 

206 overwrite_b = overwrite_b or _datacopied(b1, b) 

207 

208 potrs, = get_lapack_funcs(('potrs',), (c, b1)) 

209 x, info = potrs(c, b1, lower=lower, overwrite_b=overwrite_b) 

210 if info != 0: 

211 raise ValueError('illegal value in %dth argument of internal potrs' 

212 % -info) 

213 return x 

214 

215 

216def cholesky_banded(ab, overwrite_ab=False, lower=False, check_finite=True): 

217 """ 

218 Cholesky decompose a banded Hermitian positive-definite matrix 

219 

220 The matrix a is stored in ab either in lower-diagonal or upper- 

221 diagonal ordered form:: 

222 

223 ab[u + i - j, j] == a[i,j] (if upper form; i <= j) 

224 ab[ i - j, j] == a[i,j] (if lower form; i >= j) 

225 

226 Example of ab (shape of a is (6,6), u=2):: 

227 

228 upper form: 

229 * * a02 a13 a24 a35 

230 * a01 a12 a23 a34 a45 

231 a00 a11 a22 a33 a44 a55 

232 

233 lower form: 

234 a00 a11 a22 a33 a44 a55 

235 a10 a21 a32 a43 a54 * 

236 a20 a31 a42 a53 * * 

237 

238 Parameters 

239 ---------- 

240 ab : (u + 1, M) array_like 

241 Banded matrix 

242 overwrite_ab : bool, optional 

243 Discard data in ab (may enhance performance) 

244 lower : bool, optional 

245 Is the matrix in the lower form. (Default is upper form) 

246 check_finite : bool, optional 

247 Whether to check that the input matrix contains only finite numbers. 

248 Disabling may give a performance gain, but may result in problems 

249 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

250 

251 Returns 

252 ------- 

253 c : (u + 1, M) ndarray 

254 Cholesky factorization of a, in the same banded format as ab 

255 

256 See Also 

257 -------- 

258 cho_solve_banded : 

259 Solve a linear set equations, given the Cholesky factorization 

260 of a banded Hermitian. 

261 

262 Examples 

263 -------- 

264 >>> import numpy as np 

265 >>> from scipy.linalg import cholesky_banded 

266 >>> from numpy import allclose, zeros, diag 

267 >>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]]) 

268 >>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1) 

269 >>> A = A + A.conj().T + np.diag(Ab[2, :]) 

270 >>> c = cholesky_banded(Ab) 

271 >>> C = np.diag(c[0, 2:], k=2) + np.diag(c[1, 1:], k=1) + np.diag(c[2, :]) 

272 >>> np.allclose(C.conj().T @ C - A, np.zeros((5, 5))) 

273 True 

274 

275 """ 

276 if check_finite: 

277 ab = asarray_chkfinite(ab) 

278 else: 

279 ab = asarray(ab) 

280 

281 pbtrf, = get_lapack_funcs(('pbtrf',), (ab,)) 

282 c, info = pbtrf(ab, lower=lower, overwrite_ab=overwrite_ab) 

283 if info > 0: 

284 raise LinAlgError("%d-th leading minor not positive definite" % info) 

285 if info < 0: 

286 raise ValueError('illegal value in %d-th argument of internal pbtrf' 

287 % -info) 

288 return c 

289 

290 

291def cho_solve_banded(cb_and_lower, b, overwrite_b=False, check_finite=True): 

292 """ 

293 Solve the linear equations ``A x = b``, given the Cholesky factorization of 

294 the banded Hermitian ``A``. 

295 

296 Parameters 

297 ---------- 

298 (cb, lower) : tuple, (ndarray, bool) 

299 `cb` is the Cholesky factorization of A, as given by cholesky_banded. 

300 `lower` must be the same value that was given to cholesky_banded. 

301 b : array_like 

302 Right-hand side 

303 overwrite_b : bool, optional 

304 If True, the function will overwrite the values in `b`. 

305 check_finite : bool, optional 

306 Whether to check that the input matrices contain only finite numbers. 

307 Disabling may give a performance gain, but may result in problems 

308 (crashes, non-termination) if the inputs do contain infinities or NaNs. 

309 

310 Returns 

311 ------- 

312 x : array 

313 The solution to the system A x = b 

314 

315 See Also 

316 -------- 

317 cholesky_banded : Cholesky factorization of a banded matrix 

318 

319 Notes 

320 ----- 

321 

322 .. versionadded:: 0.8.0 

323 

324 Examples 

325 -------- 

326 >>> import numpy as np 

327 >>> from scipy.linalg import cholesky_banded, cho_solve_banded 

328 >>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]]) 

329 >>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1) 

330 >>> A = A + A.conj().T + np.diag(Ab[2, :]) 

331 >>> c = cholesky_banded(Ab) 

332 >>> x = cho_solve_banded((c, False), np.ones(5)) 

333 >>> np.allclose(A @ x - np.ones(5), np.zeros(5)) 

334 True 

335 

336 """ 

337 (cb, lower) = cb_and_lower 

338 if check_finite: 

339 cb = asarray_chkfinite(cb) 

340 b = asarray_chkfinite(b) 

341 else: 

342 cb = asarray(cb) 

343 b = asarray(b) 

344 

345 # Validate shapes. 

346 if cb.shape[-1] != b.shape[0]: 

347 raise ValueError("shapes of cb and b are not compatible.") 

348 

349 pbtrs, = get_lapack_funcs(('pbtrs',), (cb, b)) 

350 x, info = pbtrs(cb, b, lower=lower, overwrite_b=overwrite_b) 

351 if info > 0: 

352 raise LinAlgError("%dth leading minor not positive definite" % info) 

353 if info < 0: 

354 raise ValueError('illegal value in %dth argument of internal pbtrs' 

355 % -info) 

356 return x