Coverage for /usr/lib/python3/dist-packages/sympy/matrices/decompositions.py: 7%

1import copy

3from sympy.core import S

4from sympy.core.function import expand_mul

5from sympy.functions.elementary.miscellaneous import Min, sqrt

6from sympy.functions.elementary.complexes import sign

8from .common import NonSquareMatrixError, NonPositiveDefiniteMatrixError

9from .utilities import _get_intermediate_simp, _iszero

10from .determinant import _find_reasonable_pivot_naive

13def _rank_decomposition(M, iszerofunc=_iszero, simplify=False):

14 r"""Returns a pair of matrices (`C`, `F`) with matching rank

15 such that `A = C F`.

17 Parameters

18 ==========

20 iszerofunc : Function, optional

21 A function used for detecting whether an element can

22 act as a pivot. ``lambda x: x.is_zero`` is used by default.

24 simplify : Bool or Function, optional

25 A function used to simplify elements when looking for a

26 pivot. By default SymPy's ``simplify`` is used.

28 Returns

29 =======

31 (C, F) : Matrices

32 `C` and `F` are full-rank matrices with rank as same as `A`,

33 whose product gives `A`.

35 See Notes for additional mathematical details.

37 Examples

38 ========

40 >>> from sympy import Matrix

41 >>> A = Matrix([

42 ... [1, 3, 1, 4],

43 ... [2, 7, 3, 9],

44 ... [1, 5, 3, 1],

45 ... [1, 2, 0, 8]

46 ... ])

47 >>> C, F = A.rank_decomposition()

48 >>> C

49 Matrix([

50 [1, 3, 4],

51 [2, 7, 9],

52 [1, 5, 1],

53 [1, 2, 8]])

54 >>> F

55 Matrix([

56 [1, 0, -2, 0],

57 [0, 1, 1, 0],

58 [0, 0, 0, 1]])

59 >>> C * F == A

60 True

62 Notes

63 =====

65 Obtaining `F`, an RREF of `A`, is equivalent to creating a

66 product

68 .. math::

69 E_n E_{n-1} ... E_1 A = F

71 where `E_n, E_{n-1}, \dots, E_1` are the elimination matrices or

72 permutation matrices equivalent to each row-reduction step.

74 The inverse of the same product of elimination matrices gives

75 `C`:

77 .. math::

78 C = \left(E_n E_{n-1} \dots E_1\right)^{-1}

80 It is not necessary, however, to actually compute the inverse:

81 the columns of `C` are those from the original matrix with the

82 same column indices as the indices of the pivot columns of `F`.

84 References

85 ==========

87 .. [1] https://en.wikipedia.org/wiki/Rank_factorization

89 .. [2] Piziak, R.; Odell, P. L. (1 June 1999).

90 "Full Rank Factorization of Matrices".

91 Mathematics Magazine. 72 (3): 193. doi:10.2307/2690882

93 See Also

94 ========

96 sympy.matrices.matrices.MatrixReductions.rref

97 """

99 F, pivot_cols = M.rref(simplify=simplify, iszerofunc=iszerofunc,

100 pivots=True)

101 rank = len(pivot_cols)

102

103 C = M.extract(range(M.rows), pivot_cols)

104 F = F[:rank, :]

105

106 return C, F

107

108

109def _liupc(M):

110 """Liu's algorithm, for pre-determination of the Elimination Tree of

111 the given matrix, used in row-based symbolic Cholesky factorization.

112

113 Examples

114 ========

115

116 >>> from sympy import SparseMatrix

117 >>> S = SparseMatrix([

118 ... [1, 0, 3, 2],

119 ... [0, 0, 1, 0],

120 ... [4, 0, 0, 5],

121 ... [0, 6, 7, 0]])

122 >>> S.liupc()

123 ([[0], [], [0], [1, 2]], [4, 3, 4, 4])

124

125 References

126 ==========

127

128 .. [1] Symbolic Sparse Cholesky Factorization using Elimination Trees,

129 Jeroen Van Grondelle (1999)

130 https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.39.7582

131 """

132 # Algorithm 2.4, p 17 of reference

133

134 # get the indices of the elements that are non-zero on or below diag

135 R = [[] for r in range(M.rows)]

136

137 for r, c, _ in M.row_list():

138 if c <= r:

139 R[r].append(c)

140

141 inf = len(R) # nothing will be this large

142 parent = [inf]*M.rows

143 virtual = [inf]*M.rows

144

145 for r in range(M.rows):

146 for c in R[r][:-1]:

147 while virtual[c] < r:

148 t = virtual[c]

149 virtual[c] = r

150 c = t

151

152 if virtual[c] == inf:

153 parent[c] = virtual[c] = r

154

155 return R, parent

156

157def _row_structure_symbolic_cholesky(M):

158 """Symbolic cholesky factorization, for pre-determination of the

159 non-zero structure of the Cholesky factororization.

160

161 Examples

162 ========

163

164 >>> from sympy import SparseMatrix

165 >>> S = SparseMatrix([

166 ... [1, 0, 3, 2],

167 ... [0, 0, 1, 0],

168 ... [4, 0, 0, 5],

169 ... [0, 6, 7, 0]])

170 >>> S.row_structure_symbolic_cholesky()

171 [[0], [], [0], [1, 2]]

172

173 References

174 ==========

175

176 .. [1] Symbolic Sparse Cholesky Factorization using Elimination Trees,

177 Jeroen Van Grondelle (1999)

178 https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.39.7582

179 """

180

181 R, parent = M.liupc()

182 inf = len(R) # this acts as infinity

183 Lrow = copy.deepcopy(R)

184

185 for k in range(M.rows):

186 for j in R[k]:

187 while j != inf and j != k:

188 Lrow[k].append(j)

189 j = parent[j]

190

191 Lrow[k] = sorted(set(Lrow[k]))

192

193 return Lrow

194

195

196def _cholesky(M, hermitian=True):

197 """Returns the Cholesky-type decomposition L of a matrix A

198 such that L * L.H == A if hermitian flag is True,

199 or L * L.T == A if hermitian is False.

200

201 A must be a Hermitian positive-definite matrix if hermitian is True,

202 or a symmetric matrix if it is False.

203

204 Examples

205 ========

206

207 >>> from sympy import Matrix

208 >>> A = Matrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))

209 >>> A.cholesky()

210 Matrix([

211 [ 5, 0, 0],

212 [ 3, 3, 0],

213 [-1, 1, 3]])

214 >>> A.cholesky() * A.cholesky().T

215 Matrix([

216 [25, 15, -5],

217 [15, 18, 0],

218 [-5, 0, 11]])

219

220 The matrix can have complex entries:

221

222 >>> from sympy import I

223 >>> A = Matrix(((9, 3*I), (-3*I, 5)))

224 >>> A.cholesky()

225 Matrix([

226 [ 3, 0],

227 [-I, 2]])

228 >>> A.cholesky() * A.cholesky().H

229 Matrix([

230 [ 9, 3*I],

231 [-3*I, 5]])

232

233 Non-hermitian Cholesky-type decomposition may be useful when the

234 matrix is not positive-definite.

235

236 >>> A = Matrix([[1, 2], [2, 1]])

237 >>> L = A.cholesky(hermitian=False)

238 >>> L

239 Matrix([

240 [1, 0],

241 [2, sqrt(3)*I]])

242 >>> L*L.T == A

243 True

244

245 See Also

246 ========

247

248 sympy.matrices.dense.DenseMatrix.LDLdecomposition

249 sympy.matrices.matrices.MatrixBase.LUdecomposition

250 QRdecomposition

251 """

252

253 from .dense import MutableDenseMatrix

254

255 if not M.is_square:

256 raise NonSquareMatrixError("Matrix must be square.")

257 if hermitian and not M.is_hermitian:

258 raise ValueError("Matrix must be Hermitian.")

259 if not hermitian and not M.is_symmetric():

260 raise ValueError("Matrix must be symmetric.")

261

262 L = MutableDenseMatrix.zeros(M.rows, M.rows)

263

264 if hermitian:

265 for i in range(M.rows):

266 for j in range(i):

267 L[i, j] = ((1 / L[j, j])*(M[i, j] -

268 sum(L[i, k]*L[j, k].conjugate() for k in range(j))))

269

270 Lii2 = (M[i, i] -

271 sum(L[i, k]*L[i, k].conjugate() for k in range(i)))

272

273 if Lii2.is_positive is False:

274 raise NonPositiveDefiniteMatrixError(

275 "Matrix must be positive-definite")

276

277 L[i, i] = sqrt(Lii2)

278

279 else:

280 for i in range(M.rows):

281 for j in range(i):

282 L[i, j] = ((1 / L[j, j])*(M[i, j] -

283 sum(L[i, k]*L[j, k] for k in range(j))))

284

285 L[i, i] = sqrt(M[i, i] -

286 sum(L[i, k]**2 for k in range(i)))

287

288 return M._new(L)

289

290def _cholesky_sparse(M, hermitian=True):

291 """

292 Returns the Cholesky decomposition L of a matrix A

293 such that L * L.T = A

294

295 A must be a square, symmetric, positive-definite

296 and non-singular matrix

297

298 Examples

299 ========

300

301 >>> from sympy import SparseMatrix

302 >>> A = SparseMatrix(((25,15,-5),(15,18,0),(-5,0,11)))

303 >>> A.cholesky()

304 Matrix([

305 [ 5, 0, 0],

306 [ 3, 3, 0],

307 [-1, 1, 3]])

308 >>> A.cholesky() * A.cholesky().T == A

309 True

310

311 The matrix can have complex entries:

312

313 >>> from sympy import I

314 >>> A = SparseMatrix(((9, 3*I), (-3*I, 5)))

315 >>> A.cholesky()

316 Matrix([

317 [ 3, 0],

318 [-I, 2]])

319 >>> A.cholesky() * A.cholesky().H

320 Matrix([

321 [ 9, 3*I],

322 [-3*I, 5]])

323

324 Non-hermitian Cholesky-type decomposition may be useful when the

325 matrix is not positive-definite.

326

327 >>> A = SparseMatrix([[1, 2], [2, 1]])

328 >>> L = A.cholesky(hermitian=False)

329 >>> L

330 Matrix([

331 [1, 0],

332 [2, sqrt(3)*I]])

333 >>> L*L.T == A

334 True

335

336 See Also

337 ========

338

339 sympy.matrices.sparse.SparseMatrix.LDLdecomposition

340 sympy.matrices.matrices.MatrixBase.LUdecomposition

341 QRdecomposition

342 """

343

344 from .dense import MutableDenseMatrix

345

346 if not M.is_square:

347 raise NonSquareMatrixError("Matrix must be square.")

348 if hermitian and not M.is_hermitian:

349 raise ValueError("Matrix must be Hermitian.")

350 if not hermitian and not M.is_symmetric():

351 raise ValueError("Matrix must be symmetric.")

352

353 dps = _get_intermediate_simp(expand_mul, expand_mul)

354 Crowstruc = M.row_structure_symbolic_cholesky()

355 C = MutableDenseMatrix.zeros(M.rows)

356

357 for i in range(len(Crowstruc)):

358 for j in Crowstruc[i]:

359 if i != j:

360 C[i, j] = M[i, j]

361 summ = 0

362

363 for p1 in Crowstruc[i]:

364 if p1 < j:

365 for p2 in Crowstruc[j]:

366 if p2 < j:

367 if p1 == p2:

368 if hermitian:

369 summ += C[i, p1]*C[j, p1].conjugate()

370 else:

371 summ += C[i, p1]*C[j, p1]

372 else:

373 break

374 else:

375 break

376

377 C[i, j] = dps((C[i, j] - summ) / C[j, j])

378

379 else: # i == j

380 C[j, j] = M[j, j]

381 summ = 0

382

383 for k in Crowstruc[j]:

384 if k < j:

385 if hermitian:

386 summ += C[j, k]*C[j, k].conjugate()

387 else:

388 summ += C[j, k]**2

389 else:

390 break

391

392 Cjj2 = dps(C[j, j] - summ)

393

394 if hermitian and Cjj2.is_positive is False:

395 raise NonPositiveDefiniteMatrixError(

396 "Matrix must be positive-definite")

397

398 C[j, j] = sqrt(Cjj2)

399

400 return M._new(C)

401

402

403def _LDLdecomposition(M, hermitian=True):

404 """Returns the LDL Decomposition (L, D) of matrix A,

405 such that L * D * L.H == A if hermitian flag is True, or

406 L * D * L.T == A if hermitian is False.

407 This method eliminates the use of square root.

408 Further this ensures that all the diagonal entries of L are 1.

409 A must be a Hermitian positive-definite matrix if hermitian is True,

410 or a symmetric matrix otherwise.

411

412 Examples

413 ========

414

415 >>> from sympy import Matrix, eye

416 >>> A = Matrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))

417 >>> L, D = A.LDLdecomposition()

418 >>> L

419 Matrix([

420 [ 1, 0, 0],

421 [ 3/5, 1, 0],

422 [-1/5, 1/3, 1]])

423 >>> D

424 Matrix([

425 [25, 0, 0],

426 [ 0, 9, 0],

427 [ 0, 0, 9]])

428 >>> L * D * L.T * A.inv() == eye(A.rows)

429 True

430

431 The matrix can have complex entries:

432

433 >>> from sympy import I

434 >>> A = Matrix(((9, 3*I), (-3*I, 5)))

435 >>> L, D = A.LDLdecomposition()

436 >>> L

437 Matrix([

438 [ 1, 0],

439 [-I/3, 1]])

440 >>> D

441 Matrix([

442 [9, 0],

443 [0, 4]])

444 >>> L*D*L.H == A

445 True

446

447 See Also

448 ========

449

450 sympy.matrices.dense.DenseMatrix.cholesky

451 sympy.matrices.matrices.MatrixBase.LUdecomposition

452 QRdecomposition

453 """

454

455 from .dense import MutableDenseMatrix

456

457 if not M.is_square:

458 raise NonSquareMatrixError("Matrix must be square.")

459 if hermitian and not M.is_hermitian:

460 raise ValueError("Matrix must be Hermitian.")

461 if not hermitian and not M.is_symmetric():

462 raise ValueError("Matrix must be symmetric.")

463

464 D = MutableDenseMatrix.zeros(M.rows, M.rows)

465 L = MutableDenseMatrix.eye(M.rows)

466

467 if hermitian:

468 for i in range(M.rows):

469 for j in range(i):

470 L[i, j] = (1 / D[j, j])*(M[i, j] - sum(

471 L[i, k]*L[j, k].conjugate()*D[k, k] for k in range(j)))

472

473 D[i, i] = (M[i, i] -

474 sum(L[i, k]*L[i, k].conjugate()*D[k, k] for k in range(i)))

475

476 if D[i, i].is_positive is False:

477 raise NonPositiveDefiniteMatrixError(

478 "Matrix must be positive-definite")

479

480 else:

481 for i in range(M.rows):

482 for j in range(i):

483 L[i, j] = (1 / D[j, j])*(M[i, j] - sum(

484 L[i, k]*L[j, k]*D[k, k] for k in range(j)))

485

486 D[i, i] = M[i, i] - sum(L[i, k]**2*D[k, k] for k in range(i))

487

488 return M._new(L), M._new(D)

489

490def _LDLdecomposition_sparse(M, hermitian=True):

491 """

492 Returns the LDL Decomposition (matrices ``L`` and ``D``) of matrix

493 ``A``, such that ``L * D * L.T == A``. ``A`` must be a square,

494 symmetric, positive-definite and non-singular.

495

496 This method eliminates the use of square root and ensures that all

497 the diagonal entries of L are 1.

498

499 Examples

500 ========

501

502 >>> from sympy import SparseMatrix

503 >>> A = SparseMatrix(((25, 15, -5), (15, 18, 0), (-5, 0, 11)))

504 >>> L, D = A.LDLdecomposition()

505 >>> L

506 Matrix([

507 [ 1, 0, 0],

508 [ 3/5, 1, 0],

509 [-1/5, 1/3, 1]])

510 >>> D

511 Matrix([

512 [25, 0, 0],

513 [ 0, 9, 0],

514 [ 0, 0, 9]])

515 >>> L * D * L.T == A

516 True

517

518 """

519

520 from .dense import MutableDenseMatrix

521

522 if not M.is_square:

523 raise NonSquareMatrixError("Matrix must be square.")

524 if hermitian and not M.is_hermitian:

525 raise ValueError("Matrix must be Hermitian.")

526 if not hermitian and not M.is_symmetric():

527 raise ValueError("Matrix must be symmetric.")

528

529 dps = _get_intermediate_simp(expand_mul, expand_mul)

530 Lrowstruc = M.row_structure_symbolic_cholesky()

531 L = MutableDenseMatrix.eye(M.rows)

532 D = MutableDenseMatrix.zeros(M.rows, M.cols)

533

534 for i in range(len(Lrowstruc)):

535 for j in Lrowstruc[i]:

536 if i != j:

537 L[i, j] = M[i, j]

538 summ = 0

539

540 for p1 in Lrowstruc[i]:

541 if p1 < j:

542 for p2 in Lrowstruc[j]:

543 if p2 < j:

544 if p1 == p2:

545 if hermitian:

546 summ += L[i, p1]*L[j, p1].conjugate()*D[p1, p1]

547 else:

548 summ += L[i, p1]*L[j, p1]*D[p1, p1]

549 else:

550 break

551 else:

552 break

553

554 L[i, j] = dps((L[i, j] - summ) / D[j, j])

555

556 else: # i == j

557 D[i, i] = M[i, i]

558 summ = 0

559

560 for k in Lrowstruc[i]:

561 if k < i:

562 if hermitian:

563 summ += L[i, k]*L[i, k].conjugate()*D[k, k]

564 else:

565 summ += L[i, k]**2*D[k, k]

566 else:

567 break

568

569 D[i, i] = dps(D[i, i] - summ)

570

571 if hermitian and D[i, i].is_positive is False:

572 raise NonPositiveDefiniteMatrixError(

573 "Matrix must be positive-definite")

574

575 return M._new(L), M._new(D)

576

577

578def _LUdecomposition(M, iszerofunc=_iszero, simpfunc=None, rankcheck=False):

579 """Returns (L, U, perm) where L is a lower triangular matrix with unit

580 diagonal, U is an upper triangular matrix, and perm is a list of row

581 swap index pairs. If A is the original matrix, then

582 ``A = (L*U).permuteBkwd(perm)``, and the row permutation matrix P such

583 that $P A = L U$ can be computed by ``P = eye(A.rows).permuteFwd(perm)``.

584

585 See documentation for LUCombined for details about the keyword argument

586 rankcheck, iszerofunc, and simpfunc.

587

588 Parameters

589 ==========

590

591 rankcheck : bool, optional

592 Determines if this function should detect the rank

593 deficiency of the matrixis and should raise a

594 ``ValueError``.

595

596 iszerofunc : function, optional

597 A function which determines if a given expression is zero.

598

599 The function should be a callable that takes a single

600 SymPy expression and returns a 3-valued boolean value

601 ``True``, ``False``, or ``None``.

602

603 It is internally used by the pivot searching algorithm.

604 See the notes section for a more information about the

605 pivot searching algorithm.

606

607 simpfunc : function or None, optional

608 A function that simplifies the input.

609

610 If this is specified as a function, this function should be

611 a callable that takes a single SymPy expression and returns

612 an another SymPy expression that is algebraically

613 equivalent.

614

615 If ``None``, it indicates that the pivot search algorithm

616 should not attempt to simplify any candidate pivots.

617

618 It is internally used by the pivot searching algorithm.

619 See the notes section for a more information about the

620 pivot searching algorithm.

621

622 Examples

623 ========

624

625 >>> from sympy import Matrix

626 >>> a = Matrix([[4, 3], [6, 3]])

627 >>> L, U, _ = a.LUdecomposition()

628 >>> L

629 Matrix([

630 [ 1, 0],

631 [3/2, 1]])

632 >>> U

633 Matrix([

634 [4, 3],

635 [0, -3/2]])

636

637 See Also

638 ========

639

640 sympy.matrices.dense.DenseMatrix.cholesky

641 sympy.matrices.dense.DenseMatrix.LDLdecomposition

642 QRdecomposition

643 LUdecomposition_Simple

644 LUdecompositionFF

645 LUsolve

646 """

647

648 combined, p = M.LUdecomposition_Simple(iszerofunc=iszerofunc,

649 simpfunc=simpfunc, rankcheck=rankcheck)

650

651 # L is lower triangular ``M.rows x M.rows``

652 # U is upper triangular ``M.rows x M.cols``

653 # L has unit diagonal. For each column in combined, the subcolumn

654 # below the diagonal of combined is shared by L.

655 # If L has more columns than combined, then the remaining subcolumns

656 # below the diagonal of L are zero.

657 # The upper triangular portion of L and combined are equal.

658 def entry_L(i, j):

659 if i < j:

660 # Super diagonal entry

661 return M.zero

662 elif i == j:

663 return M.one

664 elif j < combined.cols:

665 return combined[i, j]

666

667 # Subdiagonal entry of L with no corresponding

668 # entry in combined

669 return M.zero

670

671 def entry_U(i, j):

672 return M.zero if i > j else combined[i, j]

673

674 L = M._new(combined.rows, combined.rows, entry_L)

675 U = M._new(combined.rows, combined.cols, entry_U)

676

677 return L, U, p

678

679def _LUdecomposition_Simple(M, iszerofunc=_iszero, simpfunc=None,

680 rankcheck=False):

681 r"""Compute the PLU decomposition of the matrix.

682

683 Parameters

684 ==========

685

686 rankcheck : bool, optional

687 Determines if this function should detect the rank

688 deficiency of the matrixis and should raise a

689 ``ValueError``.

690

691 iszerofunc : function, optional

692 A function which determines if a given expression is zero.

693

694 The function should be a callable that takes a single

695 SymPy expression and returns a 3-valued boolean value

696 ``True``, ``False``, or ``None``.

697

698 It is internally used by the pivot searching algorithm.

699 See the notes section for a more information about the

700 pivot searching algorithm.

701

702 simpfunc : function or None, optional

703 A function that simplifies the input.

704

705 If this is specified as a function, this function should be

706 a callable that takes a single SymPy expression and returns

707 an another SymPy expression that is algebraically

708 equivalent.

709

710 If ``None``, it indicates that the pivot search algorithm

711 should not attempt to simplify any candidate pivots.

712

713 It is internally used by the pivot searching algorithm.

714 See the notes section for a more information about the

715 pivot searching algorithm.

716

717 Returns

718 =======

719

720 (lu, row_swaps) : (Matrix, list)

721 If the original matrix is a $m, n$ matrix:

722

723 *lu* is a $m, n$ matrix, which contains result of the

724 decomposition in a compressed form. See the notes section

725 to see how the matrix is compressed.

726

727 *row_swaps* is a $m$-element list where each element is a

728 pair of row exchange indices.

729

730 ``A = (L*U).permute_backward(perm)``, and the row

731 permutation matrix $P$ from the formula $P A = L U$ can be

732 computed by ``P=eye(A.row).permute_forward(perm)``.

733

734 Raises

735 ======

736

737 ValueError

738 Raised if ``rankcheck=True`` and the matrix is found to

739 be rank deficient during the computation.

740

741 Notes

742 =====

743

744 About the PLU decomposition:

745

746 PLU decomposition is a generalization of a LU decomposition

747 which can be extended for rank-deficient matrices.

748

749 It can further be generalized for non-square matrices, and this

750 is the notation that SymPy is using.

751

752 PLU decomposition is a decomposition of a $m, n$ matrix $A$ in

753 the form of $P A = L U$ where

754

755 * $L$ is a $m, m$ lower triangular matrix with unit diagonal

756 entries.

757 * $U$ is a $m, n$ upper triangular matrix.

758 * $P$ is a $m, m$ permutation matrix.

759

760 So, for a square matrix, the decomposition would look like:

761

762 .. math::

763 L = \begin{bmatrix}

764 1 & 0 & 0 & \cdots & 0 \\

765 L_{1, 0} & 1 & 0 & \cdots & 0 \\

766 L_{2, 0} & L_{2, 1} & 1 & \cdots & 0 \\

767 \vdots & \vdots & \vdots & \ddots & \vdots \\

768 L_{n-1, 0} & L_{n-1, 1} & L_{n-1, 2} & \cdots & 1

769 \end{bmatrix}

770

771 .. math::

772 U = \begin{bmatrix}

773 U_{0, 0} & U_{0, 1} & U_{0, 2} & \cdots & U_{0, n-1} \\

774 0 & U_{1, 1} & U_{1, 2} & \cdots & U_{1, n-1} \\

775 0 & 0 & U_{2, 2} & \cdots & U_{2, n-1} \\

776 \vdots & \vdots & \vdots & \ddots & \vdots \\

777 0 & 0 & 0 & \cdots & U_{n-1, n-1}

778 \end{bmatrix}

779

780 And for a matrix with more rows than the columns,

781 the decomposition would look like:

782

783 .. math::

784 L = \begin{bmatrix}

785 1 & 0 & 0 & \cdots & 0 & 0 & \cdots & 0 \\

786 L_{1, 0} & 1 & 0 & \cdots & 0 & 0 & \cdots & 0 \\

787 L_{2, 0} & L_{2, 1} & 1 & \cdots & 0 & 0 & \cdots & 0 \\

788 \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots

789 & \vdots \\

790 L_{n-1, 0} & L_{n-1, 1} & L_{n-1, 2} & \cdots & 1 & 0

791 & \cdots & 0 \\

792 L_{n, 0} & L_{n, 1} & L_{n, 2} & \cdots & L_{n, n-1} & 1

793 & \cdots & 0 \\

794 \vdots & \vdots & \vdots & \ddots & \vdots & \vdots

795 & \ddots & \vdots \\

796 L_{m-1, 0} & L_{m-1, 1} & L_{m-1, 2} & \cdots & L_{m-1, n-1}

797 & 0 & \cdots & 1 \\

798 \end{bmatrix}

799

800 .. math::

801 U = \begin{bmatrix}

802 U_{0, 0} & U_{0, 1} & U_{0, 2} & \cdots & U_{0, n-1} \\

803 0 & U_{1, 1} & U_{1, 2} & \cdots & U_{1, n-1} \\

804 0 & 0 & U_{2, 2} & \cdots & U_{2, n-1} \\

805 \vdots & \vdots & \vdots & \ddots & \vdots \\

806 0 & 0 & 0 & \cdots & U_{n-1, n-1} \\

807 0 & 0 & 0 & \cdots & 0 \\

808 \vdots & \vdots & \vdots & \ddots & \vdots \\

809 0 & 0 & 0 & \cdots & 0

810 \end{bmatrix}

811

812 Finally, for a matrix with more columns than the rows, the

813 decomposition would look like:

814

815 .. math::

816 L = \begin{bmatrix}

817 1 & 0 & 0 & \cdots & 0 \\

818 L_{1, 0} & 1 & 0 & \cdots & 0 \\

819 L_{2, 0} & L_{2, 1} & 1 & \cdots & 0 \\

820 \vdots & \vdots & \vdots & \ddots & \vdots \\

821 L_{m-1, 0} & L_{m-1, 1} & L_{m-1, 2} & \cdots & 1

822 \end{bmatrix}

823

824 .. math::

825 U = \begin{bmatrix}

826 U_{0, 0} & U_{0, 1} & U_{0, 2} & \cdots & U_{0, m-1}

827 & \cdots & U_{0, n-1} \\

828 0 & U_{1, 1} & U_{1, 2} & \cdots & U_{1, m-1}

829 & \cdots & U_{1, n-1} \\

830 0 & 0 & U_{2, 2} & \cdots & U_{2, m-1}

831 & \cdots & U_{2, n-1} \\

832 \vdots & \vdots & \vdots & \ddots & \vdots

833 & \cdots & \vdots \\

834 0 & 0 & 0 & \cdots & U_{m-1, m-1}

835 & \cdots & U_{m-1, n-1} \\

836 \end{bmatrix}

837

838 About the compressed LU storage:

839

840 The results of the decomposition are often stored in compressed

841 forms rather than returning $L$ and $U$ matrices individually.

842

843 It may be less intiuitive, but it is commonly used for a lot of

844 numeric libraries because of the efficiency.

845

846 The storage matrix is defined as following for this specific

847 method:

848

849 * The subdiagonal elements of $L$ are stored in the subdiagonal

850 portion of $LU$, that is $LU_{i, j} = L_{i, j}$ whenever

851 $i > j$.

852 * The elements on the diagonal of $L$ are all 1, and are not

853 explicitly stored.

854 * $U$ is stored in the upper triangular portion of $LU$, that is

855 $LU_{i, j} = U_{i, j}$ whenever $i <= j$.

856 * For a case of $m > n$, the right side of the $L$ matrix is

857 trivial to store.

858 * For a case of $m < n$, the below side of the $U$ matrix is

859 trivial to store.

860

861 So, for a square matrix, the compressed output matrix would be:

862

863 .. math::

864 LU = \begin{bmatrix}

865 U_{0, 0} & U_{0, 1} & U_{0, 2} & \cdots & U_{0, n-1} \\

866 L_{1, 0} & U_{1, 1} & U_{1, 2} & \cdots & U_{1, n-1} \\

867 L_{2, 0} & L_{2, 1} & U_{2, 2} & \cdots & U_{2, n-1} \\

868 \vdots & \vdots & \vdots & \ddots & \vdots \\

869 L_{n-1, 0} & L_{n-1, 1} & L_{n-1, 2} & \cdots & U_{n-1, n-1}

870 \end{bmatrix}

871

872 For a matrix with more rows than the columns, the compressed

873 output matrix would be:

874

875 .. math::

876 LU = \begin{bmatrix}

877 U_{0, 0} & U_{0, 1} & U_{0, 2} & \cdots & U_{0, n-1} \\

878 L_{1, 0} & U_{1, 1} & U_{1, 2} & \cdots & U_{1, n-1} \\

879 L_{2, 0} & L_{2, 1} & U_{2, 2} & \cdots & U_{2, n-1} \\

880 \vdots & \vdots & \vdots & \ddots & \vdots \\

881 L_{n-1, 0} & L_{n-1, 1} & L_{n-1, 2} & \cdots

882 & U_{n-1, n-1} \\

883 \vdots & \vdots & \vdots & \ddots & \vdots \\

884 L_{m-1, 0} & L_{m-1, 1} & L_{m-1, 2} & \cdots

885 & L_{m-1, n-1} \\

886 \end{bmatrix}

887

888 For a matrix with more columns than the rows, the compressed

889 output matrix would be:

890

891 .. math::

892 LU = \begin{bmatrix}

893 U_{0, 0} & U_{0, 1} & U_{0, 2} & \cdots & U_{0, m-1}

894 & \cdots & U_{0, n-1} \\

895 L_{1, 0} & U_{1, 1} & U_{1, 2} & \cdots & U_{1, m-1}

896 & \cdots & U_{1, n-1} \\

897 L_{2, 0} & L_{2, 1} & U_{2, 2} & \cdots & U_{2, m-1}

898 & \cdots & U_{2, n-1} \\

899 \vdots & \vdots & \vdots & \ddots & \vdots

900 & \cdots & \vdots \\

901 L_{m-1, 0} & L_{m-1, 1} & L_{m-1, 2} & \cdots & U_{m-1, m-1}

902 & \cdots & U_{m-1, n-1} \\

903 \end{bmatrix}

904

905 About the pivot searching algorithm:

906

907 When a matrix contains symbolic entries, the pivot search algorithm

908 differs from the case where every entry can be categorized as zero or

909 nonzero.

910 The algorithm searches column by column through the submatrix whose

911 top left entry coincides with the pivot position.

912 If it exists, the pivot is the first entry in the current search

913 column that iszerofunc guarantees is nonzero.

914 If no such candidate exists, then each candidate pivot is simplified

915 if simpfunc is not None.

916 The search is repeated, with the difference that a candidate may be

917 the pivot if ``iszerofunc()`` cannot guarantee that it is nonzero.

918 In the second search the pivot is the first candidate that

919 iszerofunc can guarantee is nonzero.

920 If no such candidate exists, then the pivot is the first candidate

921 for which iszerofunc returns None.

922 If no such candidate exists, then the search is repeated in the next

923 column to the right.

924 The pivot search algorithm differs from the one in ``rref()``, which

925 relies on ``_find_reasonable_pivot()``.

926 Future versions of ``LUdecomposition_simple()`` may use

927 ``_find_reasonable_pivot()``.

928

929 See Also

930 ========

931

932 sympy.matrices.matrices.MatrixBase.LUdecomposition

933 LUdecompositionFF

934 LUsolve

935 """

936

937 if rankcheck:

938 # https://github.com/sympy/sympy/issues/9796

939 pass

940

941 if S.Zero in M.shape:

942 # Define LU decomposition of a matrix with no entries as a matrix

943 # of the same dimensions with all zero entries.

944 return M.zeros(M.rows, M.cols), []

945

946 dps = _get_intermediate_simp()

947 lu = M.as_mutable()

948 row_swaps = []

949

950 pivot_col = 0

951

952 for pivot_row in range(0, lu.rows - 1):

953 # Search for pivot. Prefer entry that iszeropivot determines

954 # is nonzero, over entry that iszeropivot cannot guarantee

955 # is zero.

956 # XXX ``_find_reasonable_pivot`` uses slow zero testing. Blocked by bug #10279

957 # Future versions of LUdecomposition_simple can pass iszerofunc and simpfunc

958 # to _find_reasonable_pivot().

959 # In pass 3 of _find_reasonable_pivot(), the predicate in ``if x.equals(S.Zero):``

960 # calls sympy.simplify(), and not the simplification function passed in via

961 # the keyword argument simpfunc.

962 iszeropivot = True

963

964 while pivot_col != M.cols and iszeropivot:

965 sub_col = (lu[r, pivot_col] for r in range(pivot_row, M.rows))

966

967 pivot_row_offset, pivot_value, is_assumed_non_zero, ind_simplified_pairs =\

968 _find_reasonable_pivot_naive(sub_col, iszerofunc, simpfunc)

969

970 iszeropivot = pivot_value is None

971

972 if iszeropivot:

973 # All candidate pivots in this column are zero.

974 # Proceed to next column.

975 pivot_col += 1

976

977 if rankcheck and pivot_col != pivot_row:

978 # All entries including and below the pivot position are

979 # zero, which indicates that the rank of the matrix is

980 # strictly less than min(num rows, num cols)

981 # Mimic behavior of previous implementation, by throwing a

982 # ValueError.

983 raise ValueError("Rank of matrix is strictly less than"

984 " number of rows or columns."

985 " Pass keyword argument"

986 " rankcheck=False to compute"

987 " the LU decomposition of this matrix.")

988

989 candidate_pivot_row = None if pivot_row_offset is None else pivot_row + pivot_row_offset

990

991 if candidate_pivot_row is None and iszeropivot:

992 # If candidate_pivot_row is None and iszeropivot is True

993 # after pivot search has completed, then the submatrix

994 # below and to the right of (pivot_row, pivot_col) is

995 # all zeros, indicating that Gaussian elimination is

996 # complete.

997 return lu, row_swaps

998

999 # Update entries simplified during pivot search.

1000 for offset, val in ind_simplified_pairs:

1001 lu[pivot_row + offset, pivot_col] = val

1002

1003 if pivot_row != candidate_pivot_row:

1004 # Row swap book keeping:

1005 # Record which rows were swapped.

1006 # Update stored portion of L factor by multiplying L on the

1007 # left and right with the current permutation.

1008 # Swap rows of U.

1009 row_swaps.append([pivot_row, candidate_pivot_row])

1010

1011 # Update L.

1012 lu[pivot_row, 0:pivot_row], lu[candidate_pivot_row, 0:pivot_row] = \

1013 lu[candidate_pivot_row, 0:pivot_row], lu[pivot_row, 0:pivot_row]

1014

1015 # Swap pivot row of U with candidate pivot row.

1016 lu[pivot_row, pivot_col:lu.cols], lu[candidate_pivot_row, pivot_col:lu.cols] = \

1017 lu[candidate_pivot_row, pivot_col:lu.cols], lu[pivot_row, pivot_col:lu.cols]

1018

1019 # Introduce zeros below the pivot by adding a multiple of the

1020 # pivot row to a row under it, and store the result in the

1021 # row under it.

1022 # Only entries in the target row whose index is greater than

1023 # start_col may be nonzero.

1024 start_col = pivot_col + 1

1025

1026 for row in range(pivot_row + 1, lu.rows):

1027 # Store factors of L in the subcolumn below

1028 # (pivot_row, pivot_row).

1029 lu[row, pivot_row] = \

1030 dps(lu[row, pivot_col]/lu[pivot_row, pivot_col])

1031

1032 # Form the linear combination of the pivot row and the current

1033 # row below the pivot row that zeros the entries below the pivot.

1034 # Employing slicing instead of a loop here raises

1035 # NotImplementedError: Cannot add Zero to MutableSparseMatrix

1036 # in sympy/matrices/tests/test_sparse.py.

1037 # c = pivot_row + 1 if pivot_row == pivot_col else pivot_col

1038 for c in range(start_col, lu.cols):

1039 lu[row, c] = dps(lu[row, c] - lu[row, pivot_row]*lu[pivot_row, c])

1040

1041 if pivot_row != pivot_col:

1042 # matrix rank < min(num rows, num cols),

1043 # so factors of L are not stored directly below the pivot.

1044 # These entries are zero by construction, so don't bother

1045 # computing them.

1046 for row in range(pivot_row + 1, lu.rows):

1047 lu[row, pivot_col] = M.zero

1048

1049 pivot_col += 1

1050

1051 if pivot_col == lu.cols:

1052 # All candidate pivots are zero implies that Gaussian

1053 # elimination is complete.

1054 return lu, row_swaps

1055

1056 if rankcheck:

1057 if iszerofunc(

1058 lu[Min(lu.rows, lu.cols) - 1, Min(lu.rows, lu.cols) - 1]):

1059 raise ValueError("Rank of matrix is strictly less than"

1060 " number of rows or columns."

1061 " Pass keyword argument"

1062 " rankcheck=False to compute"

1063 " the LU decomposition of this matrix.")

1064

1065 return lu, row_swaps

1066

1067def _LUdecompositionFF(M):

1068 """Compute a fraction-free LU decomposition.

1069

1070 Returns 4 matrices P, L, D, U such that PA = L D**-1 U.

1071 If the elements of the matrix belong to some integral domain I, then all

1072 elements of L, D and U are guaranteed to belong to I.

1073

1074 See Also

1075 ========

1076

1077 sympy.matrices.matrices.MatrixBase.LUdecomposition

1078 LUdecomposition_Simple

1079 LUsolve

1080

1081 References

1082 ==========

1083

1084 .. [1] W. Zhou & D.J. Jeffrey, "Fraction-free matrix factors: new forms

1085 for LU and QR factors". Frontiers in Computer Science in China,

1086 Vol 2, no. 1, pp. 67-80, 2008.

1087 """

1088

1089 from sympy.matrices import SparseMatrix

1090

1091 zeros = SparseMatrix.zeros

1092 eye = SparseMatrix.eye

1093 n, m = M.rows, M.cols

1094 U, L, P = M.as_mutable(), eye(n), eye(n)

1095 DD = zeros(n, n)

1096 oldpivot = 1

1097

1098 for k in range(n - 1):

1099 if U[k, k] == 0:

1100 for kpivot in range(k + 1, n):

1101 if U[kpivot, k]:

1102 break

1103 else:

1104 raise ValueError("Matrix is not full rank")

1105

1106 U[k, k:], U[kpivot, k:] = U[kpivot, k:], U[k, k:]

1107 L[k, :k], L[kpivot, :k] = L[kpivot, :k], L[k, :k]

1108 P[k, :], P[kpivot, :] = P[kpivot, :], P[k, :]

1109

1110 L [k, k] = Ukk = U[k, k]

1111 DD[k, k] = oldpivot * Ukk

1112

1113 for i in range(k + 1, n):

1114 L[i, k] = Uik = U[i, k]

1115

1116 for j in range(k + 1, m):

1117 U[i, j] = (Ukk * U[i, j] - U[k, j] * Uik) / oldpivot

1118

1119 U[i, k] = 0

1120

1121 oldpivot = Ukk

1122

1123 DD[n - 1, n - 1] = oldpivot

1124

1125 return P, L, DD, U

1126

1127def _singular_value_decomposition(A):

1128 r"""Returns a Condensed Singular Value decomposition.

1129

1130 Explanation

1131 ===========

1132

1133 A Singular Value decomposition is a decomposition in the form $A = U \Sigma V$

1134 where

1135

1136 - $U, V$ are column orthogonal matrix.

1137 - $\Sigma$ is a diagonal matrix, where the main diagonal contains singular

1138 values of matrix A.

1139

1140 A column orthogonal matrix satisfies

1141 $\mathbb{I} = U^H U$ while a full orthogonal matrix satisfies

1142 relation $\mathbb{I} = U U^H = U^H U$ where $\mathbb{I}$ is an identity

1143 matrix with matching dimensions.

1144

1145 For matrices which are not square or are rank-deficient, it is

1146 sufficient to return a column orthogonal matrix because augmenting

1147 them may introduce redundant computations.

1148 In condensed Singular Value Decomposition we only return column orthogonal

1149 matrices because of this reason

1150

1151 If you want to augment the results to return a full orthogonal

1152 decomposition, you should use the following procedures.

1153

1154 - Augment the $U, V$ matrices with columns that are orthogonal to every

1155 other columns and make it square.

1156 - Augment the $\Sigma$ matrix with zero rows to make it have the same

1157 shape as the original matrix.

1158

1159 The procedure will be illustrated in the examples section.

1160

1161 Examples

1162 ========

1163

1164 we take a full rank matrix first:

1165

1166 >>> from sympy import Matrix

1167 >>> A = Matrix([[1, 2],[2,1]])

1168 >>> U, S, V = A.singular_value_decomposition()

1169 >>> U

1170 Matrix([

1171 [ sqrt(2)/2, sqrt(2)/2],

1172 [-sqrt(2)/2, sqrt(2)/2]])

1173 >>> S

1174 Matrix([

1175 [1, 0],

1176 [0, 3]])

1177 >>> V

1178 Matrix([

1179 [-sqrt(2)/2, sqrt(2)/2],

1180 [ sqrt(2)/2, sqrt(2)/2]])

1181

1182 If a matrix if square and full rank both U, V

1183 are orthogonal in both directions

1184

1185 >>> U * U.H

1186 Matrix([

1187 [1, 0],

1188 [0, 1]])

1189 >>> U.H * U

1190 Matrix([

1191 [1, 0],

1192 [0, 1]])

1193

1194 >>> V * V.H

1195 Matrix([

1196 [1, 0],

1197 [0, 1]])

1198 >>> V.H * V

1199 Matrix([

1200 [1, 0],

1201 [0, 1]])

1202 >>> A == U * S * V.H

1203 True

1204

1205 >>> C = Matrix([

1206 ... [1, 0, 0, 0, 2],

1207 ... [0, 0, 3, 0, 0],

1208 ... [0, 0, 0, 0, 0],

1209 ... [0, 2, 0, 0, 0],

1210 ... ])

1211 >>> U, S, V = C.singular_value_decomposition()

1212

1213 >>> V.H * V

1214 Matrix([

1215 [1, 0, 0],

1216 [0, 1, 0],

1217 [0, 0, 1]])

1218 >>> V * V.H

1219 Matrix([

1220 [1/5, 0, 0, 0, 2/5],

1221 [ 0, 1, 0, 0, 0],

1222 [ 0, 0, 1, 0, 0],

1223 [ 0, 0, 0, 0, 0],

1224 [2/5, 0, 0, 0, 4/5]])

1225

1226 If you want to augment the results to be a full orthogonal

1227 decomposition, you should augment $V$ with an another orthogonal

1228 column.

1229

1230 You are able to append an arbitrary standard basis that are linearly

1231 independent to every other columns and you can run the Gram-Schmidt

1232 process to make them augmented as orthogonal basis.

1233

1234 >>> V_aug = V.row_join(Matrix([[0,0,0,0,1],

1235 ... [0,0,0,1,0]]).H)

1236 >>> V_aug = V_aug.QRdecomposition()[0]

1237 >>> V_aug

1238 Matrix([

1239 [0, sqrt(5)/5, 0, -2*sqrt(5)/5, 0],

1240 [1, 0, 0, 0, 0],

1241 [0, 0, 1, 0, 0],

1242 [0, 0, 0, 0, 1],

1243 [0, 2*sqrt(5)/5, 0, sqrt(5)/5, 0]])

1244 >>> V_aug.H * V_aug

1245 Matrix([

1246 [1, 0, 0, 0, 0],

1247 [0, 1, 0, 0, 0],

1248 [0, 0, 1, 0, 0],

1249 [0, 0, 0, 1, 0],

1250 [0, 0, 0, 0, 1]])

1251 >>> V_aug * V_aug.H

1252 Matrix([

1253 [1, 0, 0, 0, 0],

1254 [0, 1, 0, 0, 0],

1255 [0, 0, 1, 0, 0],

1256 [0, 0, 0, 1, 0],

1257 [0, 0, 0, 0, 1]])

1258

1259 Similarly we augment U

1260

1261 >>> U_aug = U.row_join(Matrix([0,0,1,0]))

1262 >>> U_aug = U_aug.QRdecomposition()[0]

1263 >>> U_aug

1264 Matrix([

1265 [0, 1, 0, 0],

1266 [0, 0, 1, 0],

1267 [0, 0, 0, 1],

1268 [1, 0, 0, 0]])

1269

1270 >>> U_aug.H * U_aug

1271 Matrix([

1272 [1, 0, 0, 0],

1273 [0, 1, 0, 0],

1274 [0, 0, 1, 0],

1275 [0, 0, 0, 1]])

1276 >>> U_aug * U_aug.H

1277 Matrix([

1278 [1, 0, 0, 0],

1279 [0, 1, 0, 0],

1280 [0, 0, 1, 0],

1281 [0, 0, 0, 1]])

1282

1283 We add 2 zero columns and one row to S

1284

1285 >>> S_aug = S.col_join(Matrix([[0,0,0]]))

1286 >>> S_aug = S_aug.row_join(Matrix([[0,0,0,0],

1287 ... [0,0,0,0]]).H)

1288 >>> S_aug

1289 Matrix([

1290 [2, 0, 0, 0, 0],

1291 [0, sqrt(5), 0, 0, 0],

1292 [0, 0, 3, 0, 0],

1293 [0, 0, 0, 0, 0]])

1294

1295

1296

1297 >>> U_aug * S_aug * V_aug.H == C

1298 True

1299

1300 """

1301

1302 AH = A.H

1303 m, n = A.shape

1304 if m >= n:

1305 V, S = (AH * A).diagonalize()

1306

1307 ranked = []

1308 for i, x in enumerate(S.diagonal()):

1309 if not x.is_zero:

1310 ranked.append(i)

1311

1312 V = V[:, ranked]

1313

1314 Singular_vals = [sqrt(S[i, i]) for i in range(S.rows) if i in ranked]

1315

1316 S = S.zeros(len(Singular_vals))

1317

1318 for i, sv in enumerate(Singular_vals):

1319 S[i, i] = sv

1320

1321 V, _ = V.QRdecomposition()

1322 U = A * V * S.inv()

1323 else:

1324 U, S = (A * AH).diagonalize()

1325

1326 ranked = []

1327 for i, x in enumerate(S.diagonal()):

1328 if not x.is_zero:

1329 ranked.append(i)

1330

1331 U = U[:, ranked]

1332 Singular_vals = [sqrt(S[i, i]) for i in range(S.rows) if i in ranked]

1333

1334 S = S.zeros(len(Singular_vals))

1335

1336 for i, sv in enumerate(Singular_vals):

1337 S[i, i] = sv

1338

1339 U, _ = U.QRdecomposition()

1340 V = AH * U * S.inv()

1341

1342 return U, S, V

1343

1344def _QRdecomposition_optional(M, normalize=True):

1345 def dot(u, v):

1346 return u.dot(v, hermitian=True)

1347

1348 dps = _get_intermediate_simp(expand_mul, expand_mul)

1349

1350 A = M.as_mutable()

1351 ranked = []

1352

1353 Q = A

1354 R = A.zeros(A.cols)

1355

1356 for j in range(A.cols):

1357 for i in range(j):

1358 if Q[:, i].is_zero_matrix:

1359 continue

1360

1361 R[i, j] = dot(Q[:, i], Q[:, j]) / dot(Q[:, i], Q[:, i])

1362 R[i, j] = dps(R[i, j])

1363 Q[:, j] -= Q[:, i] * R[i, j]

1364

1365 Q[:, j] = dps(Q[:, j])

1366 if Q[:, j].is_zero_matrix is not True:

1367 ranked.append(j)

1368 R[j, j] = M.one

1369

1370 Q = Q.extract(range(Q.rows), ranked)

1371 R = R.extract(ranked, range(R.cols))

1372

1373 if normalize:

1374 # Normalization

1375 for i in range(Q.cols):

1376 norm = Q[:, i].norm()

1377 Q[:, i] /= norm

1378 R[i, :] *= norm

1379

1380 return M.__class__(Q), M.__class__(R)

1381

1382

1383def _QRdecomposition(M):

1384 r"""Returns a QR decomposition.

1385

1386 Explanation

1387 ===========

1388

1389 A QR decomposition is a decomposition in the form $A = Q R$

1390 where

1391

1392 - $Q$ is a column orthogonal matrix.

1393 - $R$ is a upper triangular (trapezoidal) matrix.

1394

1395 A column orthogonal matrix satisfies

1396 $\mathbb{I} = Q^H Q$ while a full orthogonal matrix satisfies

1397 relation $\mathbb{I} = Q Q^H = Q^H Q$ where $I$ is an identity

1398 matrix with matching dimensions.

1399

1400 For matrices which are not square or are rank-deficient, it is

1401 sufficient to return a column orthogonal matrix because augmenting

1402 them may introduce redundant computations.

1403 And an another advantage of this is that you can easily inspect the

1404 matrix rank by counting the number of columns of $Q$.

1405

1406 If you want to augment the results to return a full orthogonal

1407 decomposition, you should use the following procedures.

1408

1409 - Augment the $Q$ matrix with columns that are orthogonal to every

1410 other columns and make it square.

1411 - Augment the $R$ matrix with zero rows to make it have the same

1412 shape as the original matrix.

1413

1414 The procedure will be illustrated in the examples section.

1415

1416 Examples

1417 ========

1418

1419 A full rank matrix example:

1420

1421 >>> from sympy import Matrix

1422 >>> A = Matrix([[12, -51, 4], [6, 167, -68], [-4, 24, -41]])

1423 >>> Q, R = A.QRdecomposition()

1424 >>> Q

1425 Matrix([

1426 [ 6/7, -69/175, -58/175],

1427 [ 3/7, 158/175, 6/175],

1428 [-2/7, 6/35, -33/35]])

1429 >>> R

1430 Matrix([

1431 [14, 21, -14],

1432 [ 0, 175, -70],

1433 [ 0, 0, 35]])

1434

1435 If the matrix is square and full rank, the $Q$ matrix becomes

1436 orthogonal in both directions, and needs no augmentation.

1437

1438 >>> Q * Q.H

1439 Matrix([

1440 [1, 0, 0],

1441 [0, 1, 0],

1442 [0, 0, 1]])

1443 >>> Q.H * Q

1444 Matrix([

1445 [1, 0, 0],

1446 [0, 1, 0],

1447 [0, 0, 1]])

1448

1449 >>> A == Q*R

1450 True

1451

1452 A rank deficient matrix example:

1453

1454 >>> A = Matrix([[12, -51, 0], [6, 167, 0], [-4, 24, 0]])

1455 >>> Q, R = A.QRdecomposition()

1456 >>> Q

1457 Matrix([

1458 [ 6/7, -69/175],

1459 [ 3/7, 158/175],

1460 [-2/7, 6/35]])

1461 >>> R

1462 Matrix([

1463 [14, 21, 0],

1464 [ 0, 175, 0]])

1465

1466 QRdecomposition might return a matrix Q that is rectangular.

1467 In this case the orthogonality condition might be satisfied as

1468 $\mathbb{I} = Q.H*Q$ but not in the reversed product

1469 $\mathbb{I} = Q * Q.H$.

1470

1471 >>> Q.H * Q

1472 Matrix([

1473 [1, 0],

1474 [0, 1]])

1475 >>> Q * Q.H

1476 Matrix([

1477 [27261/30625, 348/30625, -1914/6125],

1478 [ 348/30625, 30589/30625, 198/6125],

1479 [ -1914/6125, 198/6125, 136/1225]])

1480

1481 If you want to augment the results to be a full orthogonal

1482 decomposition, you should augment $Q$ with an another orthogonal

1483 column.

1484

1485 You are able to append an arbitrary standard basis that are linearly

1486 independent to every other columns and you can run the Gram-Schmidt

1487 process to make them augmented as orthogonal basis.

1488

1489 >>> Q_aug = Q.row_join(Matrix([0, 0, 1]))

1490 >>> Q_aug = Q_aug.QRdecomposition()[0]

1491 >>> Q_aug

1492 Matrix([

1493 [ 6/7, -69/175, 58/175],

1494 [ 3/7, 158/175, -6/175],

1495 [-2/7, 6/35, 33/35]])

1496 >>> Q_aug.H * Q_aug

1497 Matrix([

1498 [1, 0, 0],

1499 [0, 1, 0],

1500 [0, 0, 1]])

1501 >>> Q_aug * Q_aug.H

1502 Matrix([

1503 [1, 0, 0],

1504 [0, 1, 0],

1505 [0, 0, 1]])

1506

1507 Augmenting the $R$ matrix with zero row is straightforward.

1508

1509 >>> R_aug = R.col_join(Matrix([[0, 0, 0]]))

1510 >>> R_aug

1511 Matrix([

1512 [14, 21, 0],

1513 [ 0, 175, 0],

1514 [ 0, 0, 0]])

1515 >>> Q_aug * R_aug == A

1516 True

1517

1518 A zero matrix example:

1519

1520 >>> from sympy import Matrix

1521 >>> A = Matrix.zeros(3, 4)

1522 >>> Q, R = A.QRdecomposition()

1523

1524 They may return matrices with zero rows and columns.

1525

1526 >>> Q

1527 Matrix(3, 0, [])

1528 >>> R

1529 Matrix(0, 4, [])

1530 >>> Q*R

1531 Matrix([

1532 [0, 0, 0, 0],

1533 [0, 0, 0, 0],

1534 [0, 0, 0, 0]])

1535

1536 As the same augmentation rule described above, $Q$ can be augmented

1537 with columns of an identity matrix and $R$ can be augmented with

1538 rows of a zero matrix.

1539

1540 >>> Q_aug = Q.row_join(Matrix.eye(3))

1541 >>> R_aug = R.col_join(Matrix.zeros(3, 4))

1542 >>> Q_aug * Q_aug.T

1543 Matrix([

1544 [1, 0, 0],

1545 [0, 1, 0],

1546 [0, 0, 1]])

1547 >>> R_aug

1548 Matrix([

1549 [0, 0, 0, 0],

1550 [0, 0, 0, 0],

1551 [0, 0, 0, 0]])

1552 >>> Q_aug * R_aug == A

1553 True

1554

1555 See Also

1556 ========

1557

1558 sympy.matrices.dense.DenseMatrix.cholesky

1559 sympy.matrices.dense.DenseMatrix.LDLdecomposition

1560 sympy.matrices.matrices.MatrixBase.LUdecomposition

1561 QRsolve

1562 """

1563 return _QRdecomposition_optional(M, normalize=True)

1564

1565def _upper_hessenberg_decomposition(A):

1566 """Converts a matrix into Hessenberg matrix H.

1567

1568 Returns 2 matrices H, P s.t.

1569 $P H P^{T} = A$, where H is an upper hessenberg matrix

1570 and P is an orthogonal matrix

1571

1572 Examples

1573 ========

1574

1575 >>> from sympy import Matrix

1576 >>> A = Matrix([

1577 ... [1,2,3],

1578 ... [-3,5,6],

1579 ... [4,-8,9],

1580 ... ])

1581 >>> H, P = A.upper_hessenberg_decomposition()

1582 >>> H

1583 Matrix([

1584 [1, 6/5, 17/5],

1585 [5, 213/25, -134/25],

1586 [0, 216/25, 137/25]])

1587 >>> P

1588 Matrix([

1589 [1, 0, 0],

1590 [0, -3/5, 4/5],

1591 [0, 4/5, 3/5]])

1592 >>> P * H * P.H == A

1593 True

1594

1595

1596 References

1597 ==========

1598

1599 .. [#] https://mathworld.wolfram.com/HessenbergDecomposition.html

1600 """

1601

1602 M = A.as_mutable()

1603

1604 if not M.is_square:

1605 raise NonSquareMatrixError("Matrix must be square.")

1606

1607 n = M.cols

1608 P = M.eye(n)

1609 H = M

1610

1611 for j in range(n - 2):

1612

1613 u = H[j + 1:, j]

1614

1615 if u[1:, :].is_zero_matrix:

1616 continue

1617

1618 if sign(u[0]) != 0:

1619 u[0] = u[0] + sign(u[0]) * u.norm()

1620 else:

1621 u[0] = u[0] + u.norm()

1622

1623 v = u / u.norm()

1624

1625 H[j + 1:, :] = H[j + 1:, :] - 2 * v * (v.H * H[j + 1:, :])

1626 H[:, j + 1:] = H[:, j + 1:] - (H[:, j + 1:] * (2 * v)) * v.H

1627 P[:, j + 1:] = P[:, j + 1:] - (P[:, j + 1:] * (2 * v)) * v.H

1628

1629 return H, P