[Numpy-discussion] Adding a linear system type to NumPy?

[Dag Sverre Seljebotn]

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 at 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 at scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20110717/3d811f9d/attachment.html>