Here is an improved patch. Vinzent
2009/8/13 Ondrej Certik <[email protected]> > > thanks! > > On Thu, Aug 13, 2009 at 6:15 AM, Vinzent > Steinberg<[email protected]> wrote: > > Jochen, any progress on this? Your patch is ready to merge in my opinion, > if > > you don't have have time to fix the docstrings, I'll do this trivial > stuff. > > :) > > > > Vinzent > > > > 2009/5/20 Ondrej Certik <[email protected]> > >> > >> On Fri, May 15, 2009 at 4:15 AM, Vinzent Steinberg > >> <[email protected]> wrote: > >> > 2009/5/15 Jochen Voß <[email protected]> > >> >> > >> >> Hi, > >> >> > >> >> I'll be travelling for a few days and will be able to send a new > patch > >> >> only sometimes next week. What I suggest to do is: > >> >> > >> >> 1) Remove the "more stable" wording from the QRsolve docstring. > >> >> 2) Add warnings to the both the LUsolve and the QRsolve docstrings, > >> >> stating that these functions should not be used for real matrics (and > >> >> point to mpmath instead). > >> >> 3) Add one or more examples with symbolic matrices to the tests. > >> >> > >> >> Does this sound ok for everybody? > >> > > >> > > >> > Yeah, this sounds very reasonable, let's do it. > >> > > >> > Have a nice journey! > >> > >> Please ping me if you have something ready for review/merge. > >> > >> Thanks, > >> Ondrej > >> > >> > > > > > > > > > > > > > --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy-patches" group. To post to this group, send email to [email protected] To unsubscribe from this group, send email to [email protected] For more options, visit this group at http://groups.google.com/group/sympy-patches?hl=en -~----------~----~----~----~------~----~------~--~---
From f9a2eb9f5a1b7547d75e60c5679d776c2de75edb Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg <[email protected]> Date: Sat, 7 Nov 2009 23:47:25 +0100 Subject: [PATCH] Add Matrix.QRsolve method This method is for educational purposes, it has not advantage over LUsolve(), except for floating-point arithmetics, for which mpmath's optimized functions should be used anyway. Also improved some matrix docstrings. Most of this was implemented by Jochen Voss, thanks! Signed-off-by: Vinzent Steinberg <[email protected]> --- sympy/matrices/matrices.py | 41 ++++++++++++++++++++++++++++++-- sympy/matrices/tests/test_matrices.py | 27 ++++++++++++++++++++- 2 files changed, 63 insertions(+), 5 deletions(-) diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py index 7735e75..b5a5811 100644 --- a/sympy/matrices/matrices.py +++ b/sympy/matrices/matrices.py @@ -802,8 +802,11 @@ def print_nonzero (self, symb="X"): def LUsolve(self, rhs, iszerofunc=_iszero): """ - Solve the linear system Ax = b. + Solve the linear system Ax = b for x. self is the coefficient matrix A and rhs is the right side b. + + This is for symbolic matrices, for real or complex ones use + sympy.mpmath.lu_solve or sympy.mpmath.qr_solve. """ assert rhs.rows == self.rows A, perm = self.LUdecomposition_Simple(iszerofunc=_iszero) @@ -985,7 +988,7 @@ def jacobian(self, X): def QRdecomposition(self): """ - Return Q*R where Q is orthogonal and R is upper triangular. + Return Q,R where A = Q*R, Q is orthogonal and R is upper triangular. Assumes full-rank square (for now). """ @@ -1006,10 +1009,42 @@ def QRdecomposition(self): R[i,j] = Q[:,i].dot(self[:,j]) return Q,R - # TODO: QRsolve + def QRsolve(self, b): + """ + Solve the linear system 'Ax = b'. + + 'self' is the matrix 'A', the method argument is the vector + 'b'. The method returns the solution vector 'x'. If 'b' is a + matrix, the system is solved for each column of 'b' and the + return value is a matrix of the same shape as 'b'. + + This method is slower (approximately by a factor of 2) but + more stable for floating-point arithmetic than the LUsolve method. + However, LUsolve usually uses an exact arithmetic, so you don't need + to use QRsolve. + + This is mainly for educational purposes and symbolic matrices, for real + (or complex) matrices use sympy.mpmath.qr_solve. + """ + + Q, R = self.QRdecomposition() + y = Q.T * b + + # back substitution to solve R*x = y: + # We build up the result "backwards" in the vector 'x' and reverse it + # only in the end. + x = [] + n = R.rows + for j in range(n-1, -1, -1): + tmp = y[j,:] + for k in range(j+1, n): + tmp -= R[j,k] * x[n-1-k] + x.append(tmp/R[j,j]) + return Matrix([row.mat for row in reversed(x)]) # Utility functions def simplify(self): + """Simplify the elements of a matrix in place.""" for i in xrange(len(self.mat)): self.mat[i] = simplify(self.mat[i]) diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py index aac4f30..3c0cf02 100644 --- a/sympy/matrices/tests/test_matrices.py +++ b/sympy/matrices/tests/test_matrices.py @@ -282,8 +282,6 @@ def test_LUdecomp(): P, L, Dee, U = M.LUdecompositionFF() assert P*M == L*Dee.inv()*U - - def test_LUsolve(): A = Matrix([[2,3,5], [3,6,2], @@ -300,6 +298,31 @@ def test_LUsolve(): soln = A.LUsolve(b) assert soln == x +def test_QRsolve(): + A = Matrix([[2,3,5], + [3,6,2], + [8,3,6]]) + x = Matrix(3,1,[3,7,5]) + b = A*x + soln = A.QRsolve(b) + assert soln == x + x = Matrix([[1,2],[3,4],[5,6]]) + b = A*x + soln = A.QRsolve(b) + assert soln == x + + A = Matrix([[0,-1,2], + [5,10,7], + [8,3,4]]) + x = Matrix(3,1,[-1,2,5]) + b = A*x + soln = A.QRsolve(b) + assert soln == x + x = Matrix([[7,8],[9,10],[11,12]]) + b = A*x + soln = A.QRsolve(b) + assert soln == x + def test_inverse(): A = eye(4) assert A.inv() == eye(4) -- 1.6.5.2
