Coverage for /usr/lib/python3/dist-packages/sympy/matrices/inverse.py: 12%

1from sympy.core.numbers import mod_inverse 

2 

3from .common import MatrixError, NonSquareMatrixError, NonInvertibleMatrixError 

4from .utilities import _iszero 

5 

6 

7def _pinv_full_rank(M): 

8 """Subroutine for full row or column rank matrices. 

9 

10 For full row rank matrices, inverse of ``A * A.H`` Exists. 

11 For full column rank matrices, inverse of ``A.H * A`` Exists. 

12 

13 This routine can apply for both cases by checking the shape 

14 and have small decision. 

15 """ 

16 

17 if M.is_zero_matrix: 

18 return M.H 

19 

20 if M.rows >= M.cols: 

21 return M.H.multiply(M).inv().multiply(M.H) 

22 else: 

23 return M.H.multiply(M.multiply(M.H).inv()) 

24 

25def _pinv_rank_decomposition(M): 

26 """Subroutine for rank decomposition 

27 

28 With rank decompositions, `A` can be decomposed into two full- 

29 rank matrices, and each matrix can take pseudoinverse 

30 individually. 

31 """ 

32 

33 if M.is_zero_matrix: 

34 return M.H 

35 

36 B, C = M.rank_decomposition() 

37 

38 Bp = _pinv_full_rank(B) 

39 Cp = _pinv_full_rank(C) 

40 

41 return Cp.multiply(Bp) 

42 

43def _pinv_diagonalization(M): 

44 """Subroutine using diagonalization 

45 

46 This routine can sometimes fail if SymPy's eigenvalue 

47 computation is not reliable. 

48 """ 

49 

50 if M.is_zero_matrix: 

51 return M.H 

52 

53 A = M 

54 AH = M.H 

55 

56 try: 

57 if M.rows >= M.cols: 

58 P, D = AH.multiply(A).diagonalize(normalize=True) 

59 D_pinv = D.applyfunc(lambda x: 0 if _iszero(x) else 1 / x) 

60 

61 return P.multiply(D_pinv).multiply(P.H).multiply(AH) 

62 

63 else: 

64 P, D = A.multiply(AH).diagonalize( 

65 normalize=True) 

66 D_pinv = D.applyfunc(lambda x: 0 if _iszero(x) else 1 / x) 

67 

68 return AH.multiply(P).multiply(D_pinv).multiply(P.H) 

69 

70 except MatrixError: 

71 raise NotImplementedError( 

72 'pinv for rank-deficient matrices where ' 

73 'diagonalization of A.H*A fails is not supported yet.') 

74 

75def _pinv(M, method='RD'): 

76 """Calculate the Moore-Penrose pseudoinverse of the matrix. 

77 

78 The Moore-Penrose pseudoinverse exists and is unique for any matrix. 

79 If the matrix is invertible, the pseudoinverse is the same as the 

80 inverse. 

81 

82 Parameters 

83 ========== 

84 

85 method : String, optional 

86 Specifies the method for computing the pseudoinverse. 

87 

88 If ``'RD'``, Rank-Decomposition will be used. 

89 

90 If ``'ED'``, Diagonalization will be used. 

91 

92 Examples 

93 ======== 

94 

95 Computing pseudoinverse by rank decomposition : 

96 

97 >>> from sympy import Matrix 

98 >>> A = Matrix([[1, 2, 3], [4, 5, 6]]) 

99 >>> A.pinv() 

100 Matrix([ 

101 [-17/18, 4/9], 

102 [ -1/9, 1/9], 

103 [ 13/18, -2/9]]) 

104 

105 Computing pseudoinverse by diagonalization : 

106 

107 >>> B = A.pinv(method='ED') 

108 >>> B.simplify() 

109 >>> B 

110 Matrix([ 

111 [-17/18, 4/9], 

112 [ -1/9, 1/9], 

113 [ 13/18, -2/9]]) 

114 

115 See Also 

116 ======== 

117 

118 inv 

119 pinv_solve 

120 

121 References 

122 ========== 

123 

124 .. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse 

125 

126 """ 

127 

128 # Trivial case: pseudoinverse of all-zero matrix is its transpose. 

129 if M.is_zero_matrix: 

130 return M.H 

131 

132 if method == 'RD': 

133 return _pinv_rank_decomposition(M) 

134 elif method == 'ED': 

135 return _pinv_diagonalization(M) 

136 else: 

137 raise ValueError('invalid pinv method %s' % repr(method)) 

138 

139 

140def _inv_mod(M, m): 

141 r""" 

142 Returns the inverse of the matrix `K` (mod `m`), if it exists. 

143 

144 Method to find the matrix inverse of `K` (mod `m`) implemented in this function: 

145 

146 * Compute `\mathrm{adj}(K) = \mathrm{cof}(K)^t`, the adjoint matrix of `K`. 

147 

148 * Compute `r = 1/\mathrm{det}(K) \pmod m`. 

149 

150 * `K^{-1} = r\cdot \mathrm{adj}(K) \pmod m`. 

151 

152 Examples 

153 ======== 

154 

155 >>> from sympy import Matrix 

156 >>> A = Matrix(2, 2, [1, 2, 3, 4]) 

157 >>> A.inv_mod(5) 

158 Matrix([ 

159 [3, 1], 

160 [4, 2]]) 

161 >>> A.inv_mod(3) 

162 Matrix([ 

163 [1, 1], 

164 [0, 1]]) 

165 

166 """ 

167 

168 if not M.is_square: 

169 raise NonSquareMatrixError() 

170 

171 N = M.cols 

172 det_K = M.det() 

173 det_inv = None 

174 

175 try: 

176 det_inv = mod_inverse(det_K, m) 

177 except ValueError: 

178 raise NonInvertibleMatrixError('Matrix is not invertible (mod %d)' % m) 

179 

180 K_adj = M.adjugate() 

181 K_inv = M.__class__(N, N, 

182 [det_inv * K_adj[i, j] % m for i in range(N) for j in range(N)]) 

183 

184 return K_inv 

185 

186 

187def _verify_invertible(M, iszerofunc=_iszero): 

188 """Initial check to see if a matrix is invertible. Raises or returns 

189 determinant for use in _inv_ADJ.""" 

190 

191 if not M.is_square: 

192 raise NonSquareMatrixError("A Matrix must be square to invert.") 

193 

194 d = M.det(method='berkowitz') 

195 zero = d.equals(0) 

196 

197 if zero is None: # if equals() can't decide, will rref be able to? 

198 ok = M.rref(simplify=True)[0] 

199 zero = any(iszerofunc(ok[j, j]) for j in range(ok.rows)) 

200 

201 if zero: 

202 raise NonInvertibleMatrixError("Matrix det == 0; not invertible.") 

203 

204 return d 

205 

206def _inv_ADJ(M, iszerofunc=_iszero): 

207 """Calculates the inverse using the adjugate matrix and a determinant. 

208 

209 See Also 

210 ======== 

211 

212 inv 

213 inverse_GE 

214 inverse_LU 

215 inverse_CH 

216 inverse_LDL 

217 """ 

218 

219 d = _verify_invertible(M, iszerofunc=iszerofunc) 

220 

221 return M.adjugate() / d 

222 

223def _inv_GE(M, iszerofunc=_iszero): 

224 """Calculates the inverse using Gaussian elimination. 

225 

226 See Also 

227 ======== 

228 

229 inv 

230 inverse_ADJ 

231 inverse_LU 

232 inverse_CH 

233 inverse_LDL 

234 """ 

235 

236 from .dense import Matrix 

237 

238 if not M.is_square: 

239 raise NonSquareMatrixError("A Matrix must be square to invert.") 

240 

241 big = Matrix.hstack(M.as_mutable(), Matrix.eye(M.rows)) 

242 red = big.rref(iszerofunc=iszerofunc, simplify=True)[0] 

243 

244 if any(iszerofunc(red[j, j]) for j in range(red.rows)): 

245 raise NonInvertibleMatrixError("Matrix det == 0; not invertible.") 

246 

247 return M._new(red[:, big.rows:]) 

248 

249def _inv_LU(M, iszerofunc=_iszero): 

250 """Calculates the inverse using LU decomposition. 

251 

252 See Also 

253 ======== 

254 

255 inv 

256 inverse_ADJ 

257 inverse_GE 

258 inverse_CH 

259 inverse_LDL 

260 """ 

261 

262 if not M.is_square: 

263 raise NonSquareMatrixError("A Matrix must be square to invert.") 

264 if M.free_symbols: 

265 _verify_invertible(M, iszerofunc=iszerofunc) 

266 

267 return M.LUsolve(M.eye(M.rows), iszerofunc=_iszero) 

268 

269def _inv_CH(M, iszerofunc=_iszero): 

270 """Calculates the inverse using cholesky decomposition. 

271 

272 See Also 

273 ======== 

274 

275 inv 

276 inverse_ADJ 

277 inverse_GE 

278 inverse_LU 

279 inverse_LDL 

280 """ 

281 

282 _verify_invertible(M, iszerofunc=iszerofunc) 

283 

284 return M.cholesky_solve(M.eye(M.rows)) 

285 

286def _inv_LDL(M, iszerofunc=_iszero): 

287 """Calculates the inverse using LDL decomposition. 

288 

289 See Also 

290 ======== 

291 

292 inv 

293 inverse_ADJ 

294 inverse_GE 

295 inverse_LU 

296 inverse_CH 

297 """ 

298 

299 _verify_invertible(M, iszerofunc=iszerofunc) 

300 

301 return M.LDLsolve(M.eye(M.rows)) 

302 

303def _inv_QR(M, iszerofunc=_iszero): 

304 """Calculates the inverse using QR decomposition. 

305 

306 See Also 

307 ======== 

308 

309 inv 

310 inverse_ADJ 

311 inverse_GE 

312 inverse_CH 

313 inverse_LDL 

314 """ 

315 

316 _verify_invertible(M, iszerofunc=iszerofunc) 

317 

318 return M.QRsolve(M.eye(M.rows)) 

319 

320def _inv_block(M, iszerofunc=_iszero): 

321 """Calculates the inverse using BLOCKWISE inversion. 

322 

323 See Also 

324 ======== 

325 

326 inv 

327 inverse_ADJ 

328 inverse_GE 

329 inverse_CH 

330 inverse_LDL 

331 """ 

332 from sympy.matrices.expressions.blockmatrix import BlockMatrix 

333 i = M.shape[0] 

334 if i <= 20 : 

335 return M.inv(method="LU", iszerofunc=_iszero) 

336 A = M[:i // 2, :i //2] 

337 B = M[:i // 2, i // 2:] 

338 C = M[i // 2:, :i // 2] 

339 D = M[i // 2:, i // 2:] 

340 try: 

341 D_inv = _inv_block(D) 

342 except NonInvertibleMatrixError: 

343 return M.inv(method="LU", iszerofunc=_iszero) 

344 B_D_i = B*D_inv 

345 BDC = B_D_i*C 

346 A_n = A - BDC 

347 try: 

348 A_n = _inv_block(A_n) 

349 except NonInvertibleMatrixError: 

350 return M.inv(method="LU", iszerofunc=_iszero) 

351 B_n = -A_n*B_D_i 

352 dc = D_inv*C 

353 C_n = -dc*A_n 

354 D_n = D_inv + dc*-B_n 

355 nn = BlockMatrix([[A_n, B_n], [C_n, D_n]]).as_explicit() 

356 return nn 

357 

358def _inv(M, method=None, iszerofunc=_iszero, try_block_diag=False): 

359 """ 

360 Return the inverse of a matrix using the method indicated. Default for 

361 dense matrices is is Gauss elimination, default for sparse matrices is LDL. 

362 

363 Parameters 

364 ========== 

365 

366 method : ('GE', 'LU', 'ADJ', 'CH', 'LDL') 

367 

368 iszerofunc : function, optional 

369 Zero-testing function to use. 

370 

371 try_block_diag : bool, optional 

372 If True then will try to form block diagonal matrices using the 

373 method get_diag_blocks(), invert these individually, and then 

374 reconstruct the full inverse matrix. 

375 

376 Examples 

377 ======== 

378 

379 >>> from sympy import SparseMatrix, Matrix 

380 >>> A = SparseMatrix([ 

381 ... [ 2, -1, 0], 

382 ... [-1, 2, -1], 

383 ... [ 0, 0, 2]]) 

384 >>> A.inv('CH') 

385 Matrix([ 

386 [2/3, 1/3, 1/6], 

387 [1/3, 2/3, 1/3], 

388 [ 0, 0, 1/2]]) 

389 >>> A.inv(method='LDL') # use of 'method=' is optional 

390 Matrix([ 

391 [2/3, 1/3, 1/6], 

392 [1/3, 2/3, 1/3], 

393 [ 0, 0, 1/2]]) 

394 >>> A * _ 

395 Matrix([ 

396 [1, 0, 0], 

397 [0, 1, 0], 

398 [0, 0, 1]]) 

399 >>> A = Matrix(A) 

400 >>> A.inv('CH') 

401 Matrix([ 

402 [2/3, 1/3, 1/6], 

403 [1/3, 2/3, 1/3], 

404 [ 0, 0, 1/2]]) 

405 >>> A.inv('ADJ') == A.inv('GE') == A.inv('LU') == A.inv('CH') == A.inv('LDL') == A.inv('QR') 

406 True 

407 

408 Notes 

409 ===== 

410 

411 According to the ``method`` keyword, it calls the appropriate method: 

412 

413 GE .... inverse_GE(); default for dense matrices 

414 LU .... inverse_LU() 

415 ADJ ... inverse_ADJ() 

416 CH ... inverse_CH() 

417 LDL ... inverse_LDL(); default for sparse matrices 

418 QR ... inverse_QR() 

419 

420 Note, the GE and LU methods may require the matrix to be simplified 

421 before it is inverted in order to properly detect zeros during 

422 pivoting. In difficult cases a custom zero detection function can 

423 be provided by setting the ``iszerofunc`` argument to a function that 

424 should return True if its argument is zero. The ADJ routine computes 

425 the determinant and uses that to detect singular matrices in addition 

426 to testing for zeros on the diagonal. 

427 

428 See Also 

429 ======== 

430 

431 inverse_ADJ 

432 inverse_GE 

433 inverse_LU 

434 inverse_CH 

435 inverse_LDL 

436 

437 Raises 

438 ====== 

439 

440 ValueError 

441 If the determinant of the matrix is zero. 

442 """ 

443 

444 from sympy.matrices import diag, SparseMatrix 

445 

446 if method is None: 

447 method = 'LDL' if isinstance(M, SparseMatrix) else 'GE' 

448 

449 if try_block_diag: 

450 blocks = M.get_diag_blocks() 

451 r = [] 

452 

453 for block in blocks: 

454 r.append(block.inv(method=method, iszerofunc=iszerofunc)) 

455 

456 return diag(*r) 

457 

458 if method == "GE": 

459 rv = M.inverse_GE(iszerofunc=iszerofunc) 

460 elif method == "LU": 

461 rv = M.inverse_LU(iszerofunc=iszerofunc) 

462 elif method == "ADJ": 

463 rv = M.inverse_ADJ(iszerofunc=iszerofunc) 

464 elif method == "CH": 

465 rv = M.inverse_CH(iszerofunc=iszerofunc) 

466 elif method == "LDL": 

467 rv = M.inverse_LDL(iszerofunc=iszerofunc) 

468 elif method == "QR": 

469 rv = M.inverse_QR(iszerofunc=iszerofunc) 

470 elif method == "BLOCK": 

471 rv = M.inverse_BLOCK(iszerofunc=iszerofunc) 

472 else: 

473 raise ValueError("Inversion method unrecognized") 

474 

475 return M._new(rv)