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

[Sturla Molden]

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