An optimized matrix multiplication library in Python (similar to Matlab) but is NOT numpy

Question

According to the NumPy documentation they may deprecate their np.matrix class. And while arrays do have their multitude of use cases, they cannot do everything. Specifically, they will "break" when doing pretty basic linear algebra operations (you can read more about it here).

Building my own matrix multiplication module in python is not too difficult, but it would not be optimized at all. I am looking for another library that has full linear algebra support which is optimized upon BLAS (Basic Linear Algebra Subprograms). Or at the least, is there any documents on how to DIY integrate a BLAS to python.

Edit: So some are suggesting the @ operator, which is like pushing a mole down a hole and having him pop up immediately in the neighbouring one. In essence, what is happening is a debuggers nightmare:

W*x == w*x.T
W@x == W@x.T

You would hope that an error is raised here letting you know that you made a mistake in defining your matrices. But since arrays don't store 2D information if they are along one axis, I am not sure that the issue can ever be solved via np.array. (These problems don't exist with np.matrix but for some reason the developers seem insistent on removing it).

Those 'faults' superficial. `np.matrix` can only be 2d, like MATLAB of the 1990s. It doesn't offer any more linear algebra power than `ndarray`, just some different syntax. The `@` operator makes matrix product just as convenient (it already makes good use of BLAS). — hpaulj, Jun 26 '20 at 23:15
@hpaulj check the edit to my post, the `np.matrix` does give more "linear algebra power" than `ndarray` — Al-Baraa El-Hag, Jun 26 '20 at 23:48
In `numpy` we often deal with more than 2 dimensions. For example I just answered a question about doing an `inner` product on the last dimension of 3d arrays: https://stackoverflow.com/questions/62603377/fastest-way-to-find-the-dot-product-of-a-large-matrix-of-vectors. `@` was introduced to handle batches of matrix products. In contrast MATLAB says: "Matrix operations follow the rules of linear algebra and are not compatible with multidimensional arrays." — hpaulj, Jun 27 '20 at 00:51
`np.einsum` is another extension to the sum-of-products task, using Einstein notation to handle a wider variety of dimension combinations. That's linear algebra applied to physics. https://en.wikipedia.org/wiki/Einstein_notation — hpaulj, Jun 27 '20 at 00:55
@hpaulj Oh wow! I didn't know that you could do that, in Matlab I just gave up and used loops when I had +2D matrices. So it seems that `ndarray` is better at 3D matrices while `np.matrix` is better at "1D" matrices with them both being equal at the 2D level (though `np.matrix` syntax is more _pythonic/intuitive_ imo). — Al-Baraa El-Hag, Jun 27 '20 at 01:07
You can also use Tensorflow https://www.tensorflow.org/api_docs/python/tf/linalg . It has the additional benefit of allowing you to use a gpu to perform the calculations. — Diego, Jun 27 '20 at 01:44

score 2 · Answer 1 · answered Jun 27 '20 at 12:56

If you insist on the distinction between column and row vectors, you can do that.

>>> x = np.array([1, 2, 3]).reshape(-1, 1)
>>> W = np.arange(15).reshape(5, 3)
>>> x
array([[1],
       [2],
       [3]])

>>> W
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11],
       [12, 13, 14]])
>>> W @ x
array([[ 8],
       [26],
       [44],
       [62],
       [80]])
>>> W @ x.T
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0,
with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 1 is different from 3)

You could create helper functions to create column and row vectors:

def rowvec(x):
    return np.array(x).reshape(1, -1)

def colvec(x):
    return np.array(x).reshape(-1, 1)

>>> rowvec([1, 2, 3])
array([[1, 2, 3]])

>>> colvec([1, 2, 3])
array([[1],
       [2],
       [3]])

I would recommend that you only use this type of constructs when you're porting existing Matlab code. You'll have trouble reading numpy code written by others and many library functions expect 1D arrays as inputs, not (1, n)-shaped arrays.

Pedro · Answer 2 · 2020-06-27T00:22:51.840

1

Actually, numpy offers BLAS-powered matrix mutiplication through the matmul operator @. This invokes the __matmul__ magic method for a given class.

All you have to do in the above example is W @ x.

Other linear algebra stuff can be found on the np.linalg module.

Edit: I guess your problem is way more about the language's style than any technical issues. I found this answer very elucidative: Transposing a NumPy array

Also, I find it very improbable that you will find something that is NOT numpy since most of the major machine learning/data science frameworks rely on it.

edited Jun 27 '20 at 00:22

answered Jun 26 '20 at 23:09

Pedro

1,121
7
16

1

So while the method you mentioned does indeed solve the immediate issue at hand, it raises a problem. W@x == W@x.T (which it shouldn't, it should raise an error like W.T@x). This will make debugging and finding this error a nightmare – Al-Baraa El-Hag Jun 26 '20 at 23:19
2

The idea that 1d arrays don't have a (distinct) transpose does take getting used too, but it's not impossible. Even when I worked in MATLAB I found that keeping track of dimensions was 80% of the debugging battle. – hpaulj Jun 26 '20 at 23:49
@hpaulj Fear is that after setting up a cascade of matrix operations (eg. `W1*W2*W3*W4`) I discover that one of them has been rendered into a 1D blob incapable of being used for further matrix math. Could I just "remember" the dimensions of my matrices while coding, inserting the needed `*` or `@` operator appropriately, most likely. As you said, I would have had to maintain them anyway with Matlab. But python throwing up an error when I make a mistake is really comforting and I would rather have it than not. – Al-Baraa El-Hag Jun 27 '20 at 00:16

An optimized matrix multiplication library in Python (similar to Matlab) but is NOT numpy

2 Answers2