Coverage for /usr/lib/python3/dist-packages/sympy/matrices/solvers.py: 8%
231 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core.function import expand_mul
2from sympy.core.symbol import Dummy, uniquely_named_symbol, symbols
3from sympy.utilities.iterables import numbered_symbols
5from .common import ShapeError, NonSquareMatrixError, NonInvertibleMatrixError
6from .eigen import _fuzzy_positive_definite
7from .utilities import _get_intermediate_simp, _iszero
10def _diagonal_solve(M, rhs):
11 """Solves ``Ax = B`` efficiently, where A is a diagonal Matrix,
12 with non-zero diagonal entries.
14 Examples
15 ========
17 >>> from sympy import Matrix, eye
18 >>> A = eye(2)*2
19 >>> B = Matrix([[1, 2], [3, 4]])
20 >>> A.diagonal_solve(B) == B/2
21 True
23 See Also
24 ========
26 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
27 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
28 gauss_jordan_solve
29 cholesky_solve
30 LDLsolve
31 LUsolve
32 QRsolve
33 pinv_solve
34 """
36 if not M.is_diagonal():
37 raise TypeError("Matrix should be diagonal")
38 if rhs.rows != M.rows:
39 raise TypeError("Size mismatch")
41 return M._new(
42 rhs.rows, rhs.cols, lambda i, j: rhs[i, j] / M[i, i])
45def _lower_triangular_solve(M, rhs):
46 """Solves ``Ax = B``, where A is a lower triangular matrix.
48 See Also
49 ========
51 upper_triangular_solve
52 gauss_jordan_solve
53 cholesky_solve
54 diagonal_solve
55 LDLsolve
56 LUsolve
57 QRsolve
58 pinv_solve
59 """
61 from .dense import MutableDenseMatrix
63 if not M.is_square:
64 raise NonSquareMatrixError("Matrix must be square.")
65 if rhs.rows != M.rows:
66 raise ShapeError("Matrices size mismatch.")
67 if not M.is_lower:
68 raise ValueError("Matrix must be lower triangular.")
70 dps = _get_intermediate_simp()
71 X = MutableDenseMatrix.zeros(M.rows, rhs.cols)
73 for j in range(rhs.cols):
74 for i in range(M.rows):
75 if M[i, i] == 0:
76 raise TypeError("Matrix must be non-singular.")
78 X[i, j] = dps((rhs[i, j] - sum(M[i, k]*X[k, j]
79 for k in range(i))) / M[i, i])
81 return M._new(X)
83def _lower_triangular_solve_sparse(M, rhs):
84 """Solves ``Ax = B``, where A is a lower triangular matrix.
86 See Also
87 ========
89 upper_triangular_solve
90 gauss_jordan_solve
91 cholesky_solve
92 diagonal_solve
93 LDLsolve
94 LUsolve
95 QRsolve
96 pinv_solve
97 """
99 if not M.is_square:
100 raise NonSquareMatrixError("Matrix must be square.")
101 if rhs.rows != M.rows:
102 raise ShapeError("Matrices size mismatch.")
103 if not M.is_lower:
104 raise ValueError("Matrix must be lower triangular.")
106 dps = _get_intermediate_simp()
107 rows = [[] for i in range(M.rows)]
109 for i, j, v in M.row_list():
110 if i > j:
111 rows[i].append((j, v))
113 X = rhs.as_mutable()
115 for j in range(rhs.cols):
116 for i in range(rhs.rows):
117 for u, v in rows[i]:
118 X[i, j] -= v*X[u, j]
120 X[i, j] = dps(X[i, j] / M[i, i])
122 return M._new(X)
125def _upper_triangular_solve(M, rhs):
126 """Solves ``Ax = B``, where A is an upper triangular matrix.
128 See Also
129 ========
131 lower_triangular_solve
132 gauss_jordan_solve
133 cholesky_solve
134 diagonal_solve
135 LDLsolve
136 LUsolve
137 QRsolve
138 pinv_solve
139 """
141 from .dense import MutableDenseMatrix
143 if not M.is_square:
144 raise NonSquareMatrixError("Matrix must be square.")
145 if rhs.rows != M.rows:
146 raise ShapeError("Matrix size mismatch.")
147 if not M.is_upper:
148 raise TypeError("Matrix is not upper triangular.")
150 dps = _get_intermediate_simp()
151 X = MutableDenseMatrix.zeros(M.rows, rhs.cols)
153 for j in range(rhs.cols):
154 for i in reversed(range(M.rows)):
155 if M[i, i] == 0:
156 raise ValueError("Matrix must be non-singular.")
158 X[i, j] = dps((rhs[i, j] - sum(M[i, k]*X[k, j]
159 for k in range(i + 1, M.rows))) / M[i, i])
161 return M._new(X)
163def _upper_triangular_solve_sparse(M, rhs):
164 """Solves ``Ax = B``, where A is an upper triangular matrix.
166 See Also
167 ========
169 lower_triangular_solve
170 gauss_jordan_solve
171 cholesky_solve
172 diagonal_solve
173 LDLsolve
174 LUsolve
175 QRsolve
176 pinv_solve
177 """
179 if not M.is_square:
180 raise NonSquareMatrixError("Matrix must be square.")
181 if rhs.rows != M.rows:
182 raise ShapeError("Matrix size mismatch.")
183 if not M.is_upper:
184 raise TypeError("Matrix is not upper triangular.")
186 dps = _get_intermediate_simp()
187 rows = [[] for i in range(M.rows)]
189 for i, j, v in M.row_list():
190 if i < j:
191 rows[i].append((j, v))
193 X = rhs.as_mutable()
195 for j in range(rhs.cols):
196 for i in reversed(range(rhs.rows)):
197 for u, v in reversed(rows[i]):
198 X[i, j] -= v*X[u, j]
200 X[i, j] = dps(X[i, j] / M[i, i])
202 return M._new(X)
205def _cholesky_solve(M, rhs):
206 """Solves ``Ax = B`` using Cholesky decomposition,
207 for a general square non-singular matrix.
208 For a non-square matrix with rows > cols,
209 the least squares solution is returned.
211 See Also
212 ========
214 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
215 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
216 gauss_jordan_solve
217 diagonal_solve
218 LDLsolve
219 LUsolve
220 QRsolve
221 pinv_solve
222 """
224 if M.rows < M.cols:
225 raise NotImplementedError(
226 'Under-determined System. Try M.gauss_jordan_solve(rhs)')
228 hermitian = True
229 reform = False
231 if M.is_symmetric():
232 hermitian = False
233 elif not M.is_hermitian:
234 reform = True
236 if reform or _fuzzy_positive_definite(M) is False:
237 H = M.H
238 M = H.multiply(M)
239 rhs = H.multiply(rhs)
240 hermitian = not M.is_symmetric()
242 L = M.cholesky(hermitian=hermitian)
243 Y = L.lower_triangular_solve(rhs)
245 if hermitian:
246 return (L.H).upper_triangular_solve(Y)
247 else:
248 return (L.T).upper_triangular_solve(Y)
251def _LDLsolve(M, rhs):
252 """Solves ``Ax = B`` using LDL decomposition,
253 for a general square and non-singular matrix.
255 For a non-square matrix with rows > cols,
256 the least squares solution is returned.
258 Examples
259 ========
261 >>> from sympy import Matrix, eye
262 >>> A = eye(2)*2
263 >>> B = Matrix([[1, 2], [3, 4]])
264 >>> A.LDLsolve(B) == B/2
265 True
267 See Also
268 ========
270 sympy.matrices.dense.DenseMatrix.LDLdecomposition
271 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
272 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
273 gauss_jordan_solve
274 cholesky_solve
275 diagonal_solve
276 LUsolve
277 QRsolve
278 pinv_solve
279 """
281 if M.rows < M.cols:
282 raise NotImplementedError(
283 'Under-determined System. Try M.gauss_jordan_solve(rhs)')
285 hermitian = True
286 reform = False
288 if M.is_symmetric():
289 hermitian = False
290 elif not M.is_hermitian:
291 reform = True
293 if reform or _fuzzy_positive_definite(M) is False:
294 H = M.H
295 M = H.multiply(M)
296 rhs = H.multiply(rhs)
297 hermitian = not M.is_symmetric()
299 L, D = M.LDLdecomposition(hermitian=hermitian)
300 Y = L.lower_triangular_solve(rhs)
301 Z = D.diagonal_solve(Y)
303 if hermitian:
304 return (L.H).upper_triangular_solve(Z)
305 else:
306 return (L.T).upper_triangular_solve(Z)
309def _LUsolve(M, rhs, iszerofunc=_iszero):
310 """Solve the linear system ``Ax = rhs`` for ``x`` where ``A = M``.
312 This is for symbolic matrices, for real or complex ones use
313 mpmath.lu_solve or mpmath.qr_solve.
315 See Also
316 ========
318 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
319 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
320 gauss_jordan_solve
321 cholesky_solve
322 diagonal_solve
323 LDLsolve
324 QRsolve
325 pinv_solve
326 LUdecomposition
327 """
329 if rhs.rows != M.rows:
330 raise ShapeError(
331 "``M`` and ``rhs`` must have the same number of rows.")
333 m = M.rows
334 n = M.cols
336 if m < n:
337 raise NotImplementedError("Underdetermined systems not supported.")
339 try:
340 A, perm = M.LUdecomposition_Simple(
341 iszerofunc=_iszero, rankcheck=True)
342 except ValueError:
343 raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
345 dps = _get_intermediate_simp()
346 b = rhs.permute_rows(perm).as_mutable()
348 # forward substitution, all diag entries are scaled to 1
349 for i in range(m):
350 for j in range(min(i, n)):
351 scale = A[i, j]
352 b.zip_row_op(i, j, lambda x, y: dps(x - y * scale))
354 # consistency check for overdetermined systems
355 if m > n:
356 for i in range(n, m):
357 for j in range(b.cols):
358 if not iszerofunc(b[i, j]):
359 raise ValueError("The system is inconsistent.")
361 b = b[0:n, :] # truncate zero rows if consistent
363 # backward substitution
364 for i in range(n - 1, -1, -1):
365 for j in range(i + 1, n):
366 scale = A[i, j]
367 b.zip_row_op(i, j, lambda x, y: dps(x - y * scale))
369 scale = A[i, i]
370 b.row_op(i, lambda x, _: dps(x / scale))
372 return rhs.__class__(b)
375def _QRsolve(M, b):
376 """Solve the linear system ``Ax = b``.
378 ``M`` is the matrix ``A``, the method argument is the vector
379 ``b``. The method returns the solution vector ``x``. If ``b`` is a
380 matrix, the system is solved for each column of ``b`` and the
381 return value is a matrix of the same shape as ``b``.
383 This method is slower (approximately by a factor of 2) but
384 more stable for floating-point arithmetic than the LUsolve method.
385 However, LUsolve usually uses an exact arithmetic, so you do not need
386 to use QRsolve.
388 This is mainly for educational purposes and symbolic matrices, for real
389 (or complex) matrices use mpmath.qr_solve.
391 See Also
392 ========
394 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
395 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
396 gauss_jordan_solve
397 cholesky_solve
398 diagonal_solve
399 LDLsolve
400 LUsolve
401 pinv_solve
402 QRdecomposition
403 """
405 dps = _get_intermediate_simp(expand_mul, expand_mul)
406 Q, R = M.QRdecomposition()
407 y = Q.T * b
409 # back substitution to solve R*x = y:
410 # We build up the result "backwards" in the vector 'x' and reverse it
411 # only in the end.
412 x = []
413 n = R.rows
415 for j in range(n - 1, -1, -1):
416 tmp = y[j, :]
418 for k in range(j + 1, n):
419 tmp -= R[j, k] * x[n - 1 - k]
421 tmp = dps(tmp)
423 x.append(tmp / R[j, j])
425 return M.vstack(*x[::-1])
428def _gauss_jordan_solve(M, B, freevar=False):
429 """
430 Solves ``Ax = B`` using Gauss Jordan elimination.
432 There may be zero, one, or infinite solutions. If one solution
433 exists, it will be returned. If infinite solutions exist, it will
434 be returned parametrically. If no solutions exist, It will throw
435 ValueError.
437 Parameters
438 ==========
440 B : Matrix
441 The right hand side of the equation to be solved for. Must have
442 the same number of rows as matrix A.
444 freevar : boolean, optional
445 Flag, when set to `True` will return the indices of the free
446 variables in the solutions (column Matrix), for a system that is
447 undetermined (e.g. A has more columns than rows), for which
448 infinite solutions are possible, in terms of arbitrary
449 values of free variables. Default `False`.
451 Returns
452 =======
454 x : Matrix
455 The matrix that will satisfy ``Ax = B``. Will have as many rows as
456 matrix A has columns, and as many columns as matrix B.
458 params : Matrix
459 If the system is underdetermined (e.g. A has more columns than
460 rows), infinite solutions are possible, in terms of arbitrary
461 parameters. These arbitrary parameters are returned as params
462 Matrix.
464 free_var_index : List, optional
465 If the system is underdetermined (e.g. A has more columns than
466 rows), infinite solutions are possible, in terms of arbitrary
467 values of free variables. Then the indices of the free variables
468 in the solutions (column Matrix) are returned by free_var_index,
469 if the flag `freevar` is set to `True`.
471 Examples
472 ========
474 >>> from sympy import Matrix
475 >>> A = Matrix([[1, 2, 1, 1], [1, 2, 2, -1], [2, 4, 0, 6]])
476 >>> B = Matrix([7, 12, 4])
477 >>> sol, params = A.gauss_jordan_solve(B)
478 >>> sol
479 Matrix([
480 [-2*tau0 - 3*tau1 + 2],
481 [ tau0],
482 [ 2*tau1 + 5],
483 [ tau1]])
484 >>> params
485 Matrix([
486 [tau0],
487 [tau1]])
488 >>> taus_zeroes = { tau:0 for tau in params }
489 >>> sol_unique = sol.xreplace(taus_zeroes)
490 >>> sol_unique
491 Matrix([
492 [2],
493 [0],
494 [5],
495 [0]])
498 >>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
499 >>> B = Matrix([3, 6, 9])
500 >>> sol, params = A.gauss_jordan_solve(B)
501 >>> sol
502 Matrix([
503 [-1],
504 [ 2],
505 [ 0]])
506 >>> params
507 Matrix(0, 1, [])
509 >>> A = Matrix([[2, -7], [-1, 4]])
510 >>> B = Matrix([[-21, 3], [12, -2]])
511 >>> sol, params = A.gauss_jordan_solve(B)
512 >>> sol
513 Matrix([
514 [0, -2],
515 [3, -1]])
516 >>> params
517 Matrix(0, 2, [])
520 >>> from sympy import Matrix
521 >>> A = Matrix([[1, 2, 1, 1], [1, 2, 2, -1], [2, 4, 0, 6]])
522 >>> B = Matrix([7, 12, 4])
523 >>> sol, params, freevars = A.gauss_jordan_solve(B, freevar=True)
524 >>> sol
525 Matrix([
526 [-2*tau0 - 3*tau1 + 2],
527 [ tau0],
528 [ 2*tau1 + 5],
529 [ tau1]])
530 >>> params
531 Matrix([
532 [tau0],
533 [tau1]])
534 >>> freevars
535 [1, 3]
538 See Also
539 ========
541 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
542 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
543 cholesky_solve
544 diagonal_solve
545 LDLsolve
546 LUsolve
547 QRsolve
548 pinv
550 References
551 ==========
553 .. [1] https://en.wikipedia.org/wiki/Gaussian_elimination
555 """
557 from sympy.matrices import Matrix, zeros
559 cls = M.__class__
560 aug = M.hstack(M.copy(), B.copy())
561 B_cols = B.cols
562 row, col = aug[:, :-B_cols].shape
564 # solve by reduced row echelon form
565 A, pivots = aug.rref(simplify=True)
566 A, v = A[:, :-B_cols], A[:, -B_cols:]
567 pivots = list(filter(lambda p: p < col, pivots))
568 rank = len(pivots)
570 # Get index of free symbols (free parameters)
571 # non-pivots columns are free variables
572 free_var_index = [c for c in range(A.cols) if c not in pivots]
574 # Bring to block form
575 permutation = Matrix(pivots + free_var_index).T
577 # check for existence of solutions
578 # rank of aug Matrix should be equal to rank of coefficient matrix
579 if not v[rank:, :].is_zero_matrix:
580 raise ValueError("Linear system has no solution")
582 # Free parameters
583 # what are current unnumbered free symbol names?
584 name = uniquely_named_symbol('tau', aug,
585 compare=lambda i: str(i).rstrip('1234567890'),
586 modify=lambda s: '_' + s).name
587 gen = numbered_symbols(name)
588 tau = Matrix([next(gen) for k in range((col - rank)*B_cols)]).reshape(
589 col - rank, B_cols)
591 # Full parametric solution
592 V = A[:rank, free_var_index]
593 vt = v[:rank, :]
594 free_sol = tau.vstack(vt - V * tau, tau)
596 # Undo permutation
597 sol = zeros(col, B_cols)
599 for k in range(col):
600 sol[permutation[k], :] = free_sol[k,:]
602 sol, tau = cls(sol), cls(tau)
604 if freevar:
605 return sol, tau, free_var_index
606 else:
607 return sol, tau
610def _pinv_solve(M, B, arbitrary_matrix=None):
611 """Solve ``Ax = B`` using the Moore-Penrose pseudoinverse.
613 There may be zero, one, or infinite solutions. If one solution
614 exists, it will be returned. If infinite solutions exist, one will
615 be returned based on the value of arbitrary_matrix. If no solutions
616 exist, the least-squares solution is returned.
618 Parameters
619 ==========
621 B : Matrix
622 The right hand side of the equation to be solved for. Must have
623 the same number of rows as matrix A.
624 arbitrary_matrix : Matrix
625 If the system is underdetermined (e.g. A has more columns than
626 rows), infinite solutions are possible, in terms of an arbitrary
627 matrix. This parameter may be set to a specific matrix to use
628 for that purpose; if so, it must be the same shape as x, with as
629 many rows as matrix A has columns, and as many columns as matrix
630 B. If left as None, an appropriate matrix containing dummy
631 symbols in the form of ``wn_m`` will be used, with n and m being
632 row and column position of each symbol.
634 Returns
635 =======
637 x : Matrix
638 The matrix that will satisfy ``Ax = B``. Will have as many rows as
639 matrix A has columns, and as many columns as matrix B.
641 Examples
642 ========
644 >>> from sympy import Matrix
645 >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
646 >>> B = Matrix([7, 8])
647 >>> A.pinv_solve(B)
648 Matrix([
649 [ _w0_0/6 - _w1_0/3 + _w2_0/6 - 55/18],
650 [-_w0_0/3 + 2*_w1_0/3 - _w2_0/3 + 1/9],
651 [ _w0_0/6 - _w1_0/3 + _w2_0/6 + 59/18]])
652 >>> A.pinv_solve(B, arbitrary_matrix=Matrix([0, 0, 0]))
653 Matrix([
654 [-55/18],
655 [ 1/9],
656 [ 59/18]])
658 See Also
659 ========
661 sympy.matrices.dense.DenseMatrix.lower_triangular_solve
662 sympy.matrices.dense.DenseMatrix.upper_triangular_solve
663 gauss_jordan_solve
664 cholesky_solve
665 diagonal_solve
666 LDLsolve
667 LUsolve
668 QRsolve
669 pinv
671 Notes
672 =====
674 This may return either exact solutions or least squares solutions.
675 To determine which, check ``A * A.pinv() * B == B``. It will be
676 True if exact solutions exist, and False if only a least-squares
677 solution exists. Be aware that the left hand side of that equation
678 may need to be simplified to correctly compare to the right hand
679 side.
681 References
682 ==========
684 .. [1] https://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#Obtaining_all_solutions_of_a_linear_system
686 """
688 from sympy.matrices import eye
690 A = M
691 A_pinv = M.pinv()
693 if arbitrary_matrix is None:
694 rows, cols = A.cols, B.cols
695 w = symbols('w:{}_:{}'.format(rows, cols), cls=Dummy)
696 arbitrary_matrix = M.__class__(cols, rows, w).T
698 return A_pinv.multiply(B) + (eye(A.cols) -
699 A_pinv.multiply(A)).multiply(arbitrary_matrix)
702def _solve(M, rhs, method='GJ'):
703 """Solves linear equation where the unique solution exists.
705 Parameters
706 ==========
708 rhs : Matrix
709 Vector representing the right hand side of the linear equation.
711 method : string, optional
712 If set to ``'GJ'`` or ``'GE'``, the Gauss-Jordan elimination will be
713 used, which is implemented in the routine ``gauss_jordan_solve``.
715 If set to ``'LU'``, ``LUsolve`` routine will be used.
717 If set to ``'QR'``, ``QRsolve`` routine will be used.
719 If set to ``'PINV'``, ``pinv_solve`` routine will be used.
721 It also supports the methods available for special linear systems
723 For positive definite systems:
725 If set to ``'CH'``, ``cholesky_solve`` routine will be used.
727 If set to ``'LDL'``, ``LDLsolve`` routine will be used.
729 To use a different method and to compute the solution via the
730 inverse, use a method defined in the .inv() docstring.
732 Returns
733 =======
735 solutions : Matrix
736 Vector representing the solution.
738 Raises
739 ======
741 ValueError
742 If there is not a unique solution then a ``ValueError`` will be
743 raised.
745 If ``M`` is not square, a ``ValueError`` and a different routine
746 for solving the system will be suggested.
747 """
749 if method in ('GJ', 'GE'):
750 try:
751 soln, param = M.gauss_jordan_solve(rhs)
753 if param:
754 raise NonInvertibleMatrixError("Matrix det == 0; not invertible. "
755 "Try ``M.gauss_jordan_solve(rhs)`` to obtain a parametric solution.")
757 except ValueError:
758 raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
760 return soln
762 elif method == 'LU':
763 return M.LUsolve(rhs)
764 elif method == 'CH':
765 return M.cholesky_solve(rhs)
766 elif method == 'QR':
767 return M.QRsolve(rhs)
768 elif method == 'LDL':
769 return M.LDLsolve(rhs)
770 elif method == 'PINV':
771 return M.pinv_solve(rhs)
772 else:
773 return M.inv(method=method).multiply(rhs)
776def _solve_least_squares(M, rhs, method='CH'):
777 """Return the least-square fit to the data.
779 Parameters
780 ==========
782 rhs : Matrix
783 Vector representing the right hand side of the linear equation.
785 method : string or boolean, optional
786 If set to ``'CH'``, ``cholesky_solve`` routine will be used.
788 If set to ``'LDL'``, ``LDLsolve`` routine will be used.
790 If set to ``'QR'``, ``QRsolve`` routine will be used.
792 If set to ``'PINV'``, ``pinv_solve`` routine will be used.
794 Otherwise, the conjugate of ``M`` will be used to create a system
795 of equations that is passed to ``solve`` along with the hint
796 defined by ``method``.
798 Returns
799 =======
801 solutions : Matrix
802 Vector representing the solution.
804 Examples
805 ========
807 >>> from sympy import Matrix, ones
808 >>> A = Matrix([1, 2, 3])
809 >>> B = Matrix([2, 3, 4])
810 >>> S = Matrix(A.row_join(B))
811 >>> S
812 Matrix([
813 [1, 2],
814 [2, 3],
815 [3, 4]])
817 If each line of S represent coefficients of Ax + By
818 and x and y are [2, 3] then S*xy is:
820 >>> r = S*Matrix([2, 3]); r
821 Matrix([
822 [ 8],
823 [13],
824 [18]])
826 But let's add 1 to the middle value and then solve for the
827 least-squares value of xy:
829 >>> xy = S.solve_least_squares(Matrix([8, 14, 18])); xy
830 Matrix([
831 [ 5/3],
832 [10/3]])
834 The error is given by S*xy - r:
836 >>> S*xy - r
837 Matrix([
838 [1/3],
839 [1/3],
840 [1/3]])
841 >>> _.norm().n(2)
842 0.58
844 If a different xy is used, the norm will be higher:
846 >>> xy += ones(2, 1)/10
847 >>> (S*xy - r).norm().n(2)
848 1.5
850 """
852 if method == 'CH':
853 return M.cholesky_solve(rhs)
854 elif method == 'QR':
855 return M.QRsolve(rhs)
856 elif method == 'LDL':
857 return M.LDLsolve(rhs)
858 elif method == 'PINV':
859 return M.pinv_solve(rhs)
860 else:
861 t = M.H
862 return (t * M).solve(t * rhs, method=method)