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