Adding a linear system type to NumPy?
The problem I am thinking we might try to might fix is that programmers with less numerical competence is unaware that the matrix expression (X**-1) * Y should be written as np.linalg.solve(X,Y) I've seen numerous times matrix expressions being typed exactly as written in linear algebra text books. Matlab does this with the backslash operator, though it is not that intuitive. Also, there seems to be general agreement against a backslash division operator in Python. So I suggest inverting a NumPy matrix could result in an unsolved linear system type, instead of actually computing the matrix inverse and returning a new matrix. The linear system type would store two arrays, A and X, symbolically representing A * (X**-1). Initially, A would be set to the indenty matrix, I. A matrix expression Y * (X**-1) would result in (1) creation of a LinearSystem object for the iversion of X, and (2) matrix multiplication of Y by A, returning a new LinearSystem object with A updated by A = Y * A. The matrix expression (X**-1) * Y Would result in (1) creation of a LinearSystem object for the iversion of X, and (2) solution of the linear system by calling np.linalg.solve, i.e. np np.linalg.solve(X,Y) The matrix expression would Z * (X**-1) * Y would form a linear system type for X, initialize A to Z, and then evaluate np.dot(A, np np.linalg.solve(X,Y)) cf. Python's evaluation order is left to right. Any other operation on a linear system (e.g. slicing) would result in formation of the inverse, by solving it against the identity matrix, set a flag that the system is solved, and then just behave as an ordinary np.matrix. Thus, (X**-1) * Y would behave as before, but do the "correct" math (instead of explicitely forming the inverse and then multiplying). Consequently this would be the same as Matlab's backslash operator, only more intuitive, as the syntax would be the same as textbook linear algebra notation. A for naming, it could e.g. be np.linear_system. I am just thinking out loudly, so forgive me for spamming the list :-) Sturla
On 7/17/2011 1:57 PM, Sturla Molden wrote:
I suggest inverting a NumPy matrix could result in an unsolved linear system type, instead of actually computing the matrix inverse and returning a new matrix.
1. Too implicit. 2. Too confusing for new users. 2a. Too confusing for students. However a "project" method might be nice, where X.project(Y) would do an orthogonal projection of Y onto X. (Then the underlying computation becomes an implementation detail.) fwiw, Alan Isaac
Something related: This autumn I expect to invest a significant amount of time (more than four weeks full-time) in a package for lazily evaluated, polymorphic linear algebra. Matrix = linear operator, a type seperate from arrays -- arrays are treated as vectors/stacked vectors Matrices can be of a variety of storage formats (diagonal, dense, the sparse formats, block-diagonal, a fortran routine that you promise acts linearly on a vector linear, and so on). The point is allowing to write linear algebra code without caring about the storage formats of the inputs. Use * for matmul, and A.I for lazy inversion i.e. (A.I * u) does a LU. In summary, a) I add my vote to this being outside the scope of numpy, b) I hope to do something about this outside of numpy. (I'll only do what is actually relevant to my research of course... But I think that will be enough for an interesting prototype of a full-fledged system for object oriented/polymorphic linear algebra) -- Sent from my Android phone with K-9 Mail. Please excuse my brevity. Sturla Molden <sturla@molden.no> wrote: The problem I am thinking we might try to might fix is that programmers with less numerical competence is unaware that the matrix expression (X**-1) * Y should be written as np.linalg.solve(X,Y) I've seen numerous times matrix expressions being typed exactly as written in linear algebra text books. Matlab does this with the backslash operator, though it is not that intuitive. Also, there seems to be general agreement against a backslash division operator in Python. So I suggest inverting a NumPy matrix could result in an unsolved linear system type, instead of actually computing the matrix inverse and returning a new matrix. The linear system type would store two arrays, A and X, symbolically representing A * (X**-1). Initially, A would be set to the indenty matrix, I. A matrix expression Y * (X**-1) would result in (1) creation of a LinearSystem object for the iversion of X, and (2) matrix multiplication of Y by A, returning a new LinearSystem object with A updated by A = Y * A. The matrix expression (X**-1) * Y Would result in (1) creation of a LinearSystem object for the iversion of X, and (2) solution of the linear system by calling np.linalg.solve, i.e. np np.linalg.solve(X,Y) The matrix expression would Z * (X**-1) * Y would form a linear system type for X, initialize A to Z, and then evaluate np.dot(A, np np.linalg.solve(X,Y)) cf. Python's evaluation order is left to right. Any other operation on a linear system (e.g. slicing) would result in formation of the inverse, by solving it against the identity matrix, set a flag that the system is solved, and then just behave as an ordinary np.matrix. Thus, (X**-1) * Y would behave as before, but do the "correct" math (instead of explicitely forming the inverse and then multiplying). Consequently this would be the same as Matlab's backslash operator, only more intuitive, as the syntax would be the same as textbook linear algebra notation. A for naming, it could e.g. be np.linear_system. I am just thinking out loudly, so forgive me for spamming the list :-) Sturla_____________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
More concrete feedback about Sturla's proposal: The problem I have is if you do A = B**-1 Then, A is some 'magic' object, not a NumPy array. That means that it is very different from Matlab's \, which restricts the context, you simply can't do A = B \ I think A.solve(u) is a lot clearer in the case of numpy arrays. Dag -- Sent from my Android phone with K-9 Mail. Please excuse my brevity. Dag Sverre Seljebotn <d.s.seljebotn@astro.uio.no> wrote: Something related: This autumn I expect to invest a significant amount of time (more than four weeks full-time) in a package for lazily evaluated, polymorphic linear algebra. Matrix = linear operator, a type seperate from arrays -- arrays are treated as vectors/stacked vectors Matrices can be of a variety of storage formats (diagonal, dense, the sparse formats, block-diagonal, a fortran routine that you promise acts linearly on a vector linear, and so on). The point is allowing to write linear algebra code without caring about the storage formats of the inputs. Use * for matmul, and A.I for lazy inversion i.e. (A.I * u) does a LU. In summary, a) I add my vote to this being outside the scope of numpy, b) I hope to do something about this outside of numpy. (I'll only do what is actually relevant to my research of course... But I think that will be enough for an interesting prototype of a full-fledged system for object oriented/polymorphic linear algebra) -- Sent from my Android phone with K-9 Mail. Please excuse my brevity. Sturla Molden <sturla@molden.no> wrote: The problem I am thinking we might try to might fix is that programmers with less numerical competence is unaware that the matrix expression (X**-1) * Y should be written as np.linalg.solve(X,Y) I've seen numerous times matrix expressions being typed exactly as written in linear algebra text books. Matlab does this with the backslash operator, though it is not that intuitive. Also, there seems to be general agreement against a backslash division operator in Python. So I suggest inverting a NumPy matrix could result in an unsolved linear system type, instead of actually computing the matrix inverse and returning a new matrix. The linear system type would store two arrays, A and X, symbolically representing A * (X**-1). Initially, A would be set to the indenty matrix, I. A matrix expression Y * (X**-1) would result in (1) creation of a LinearSystem object for the iversion of X, and (2) matrix multiplication of Y by A, returning a new LinearSystem object with A updated by A = Y * A. The matrix expression (X**-1) * Y Would result in (1) creation of a LinearSystem object for the iversion of X, and (2) solution of the linear system by calling np.linalg.solve, i.e. np np.linalg.solve(X,Y) The matrix expression would Z * (X**-1) * Y would form a linear system type for X, initialize A to Z, and then evaluate np.dot(A, np np.linalg.solve(X,Y)) cf. Python's evaluation order is left to right. Any other operation on a linear system (e.g. slicing) would result in formation of the inverse, by solving it against the identity matrix, set a flag that the system is solved, and then just behave as an ordinary np.matrix. Thus, (X**-1) * Y would behave as before, but do the "correct" math (instead of explicitely forming the inverse and then multiplying). Consequently this would be the same as Matlab's backslash operator, only more intuitive, as the syntax would be the same as textbook linear algebra notation. A for naming, it could e.g. be np.linear_system. I am just thinking out loudly, so forgive me for spamming the list :-) Sturla_____________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
participants (3)
-
Alan G Isaac -
Dag Sverre Seljebotn -
Sturla Molden