What is the status of the
matrix class in numpy? I keep being told that I should use the
ndarray class instead. Is it worth/safe using the
matrix class in new code I write? I don't understand why I should use
tl; dr: the
numpy.matrix class is getting deprecated. There are some high-profile libraries that depend on the class as a dependency (the largest one being
scipy.sparse) which hinders proper short-term deprecation of the class, but users are strongly encouraged to use the
ndarray class instead. With the introduction of the
@ operator for matrix multiplication a lot of the relative advantages of matrices have been removed.
Why (not) the matrix class?
numpy.matrix is a subclass of
numpy.ndarray. It was originally meant for convenient use in computations involving linear algebra, but there are both limitations and surprising differences in how they behave compared to instances of the more general array class. Examples for fundamental differences in behaviour:
- Shapes: arrays can have an arbitrary number of dimensions ranging from 0 infinity (or 32). Matrices are always two-dimensional. Oddly enough, while a matrix can't be created with more dimensions, it's possible to inject singleton dimensions into a matrix to end up with technically a multidimensional matrix:
np.matrix(np.random.rand(2,3))[None,...,None].shape == (1,2,3,1)(not that this is of any practical importance).
- Indexing: indexing arrays can give you arrays of any size depending on how you are indexing it. Indexing expressions on matrices will always give you a matrix. This means that both
arr[0,:]for a 2d array gives you a 1d
(1,M)in case of a
- Arithmetic operations: the main reason for using matrices in the old days was that arithmetic operations (in particular, multiplication and power) on matrices performs matrix operations (matrix multiplication and matrix power). The same for arrays results in elementwise multiplication and power. Consequently
mat1 * mat2is valid if
mat1.shape == mat2.shape, but
arr1 * arr2is valid if
arr1.shape == arr2.shape(and of course the result means something completely different). Also, surprisingly,
mat1 / mat2performs elementwise division of two matrices. This behaviour is probably inherited from
ndarraybut makes no sense for matrices, especially in light of the meaning of
- Special attributes: matrices have a few handy attributes in addition to what arrays have:
mat.A1are array views with the same value as
mat.Hare the transpose and conjugate transpose (adjoint) of the matrix;
arr.Tis the only such attribute that exists for the
mat.Iis the inverse matrix of
It's easy enough writing code that works either for ndarrays or for matrices. But when there's a chance that the two classes have to interact in code, things start to become difficult. In particular, a lot of code could work naturally for subclasses of
matrix is an ill-behaved subclass that can easily break code that tries to rely on duck typing. Consider the following example using arrays and matrices of shape
import numpy as np shape = (3, 4) arr = np.arange(np.prod(shape)).reshape(shape) # ndarray mat = np.matrix(arr) # same data in a matrix print((arr + mat).shape) # (3, 4), makes sense print((arr[0,:] + mat[0,:]).shape) # (1, 4), makes sense print((arr[:,0] + mat[:,0]).shape) # (3, 3), surprising
Adding slices of the two objects is catastrophically different depending on the dimension along which we slice. Addition on both matrices and arrays happens elementwise when the shapes are the same. The first two cases in the above are intuitive: we add two arrays (matrices), then we add two rows from each. The last case is really surprising: we probably meant to add two columns and ended up with a matrix. The reason of course is that
arr[:,0] has shape
(3,) which is compatible with shape
mat[:.0] has shape
(3,1). The two are broadcast together to shape
Finally, the largest advantage of the matrix class (i.e. the possibility to concisely formulate complicated matrix expressions involving a lot of matrix products) was removed when the
@ matmul operator was introduced in python 3.5, first implemented in numpy 1.10. Compare the computation of a simple quadratic form:
v = np.random.rand(3); v_row = np.matrix(v) arr = np.random.rand(3,3); mat = np.matrix(arr) print(v.dot(arr.dot(v))) # pre-matmul style # 0.713447037658556, yours will vary print(v_row * mat * v_row.T) # pre-matmul matrix style # [[0.71344704]] print(v @ arr @ v) # matmul style # 0.713447037658556
Looking at the above it's clear why the matrix class was widely preferred for working with linear algebra: the infix
* operator made the expressions much less verbose and much easier to read. However, we get the same readability with the
@ operator using modern python and numpy. Furthermore, note that the matrix case gives us a matrix of shape
(1,1) which should technically be a scalar. This also implies that we can't multiply a column vector with this "scalar":
(v_row * mat * v_row.T) * v_row.T in the above example raises an error because matrices with shape
(3,1) can't be multiplied in this order.
For completeness' sake it should be noted that while the matmul operator fixes the most common scenario in which ndarrays are suboptimal compared to matrices, there are still a few shortcomings in handling linear algebra elegantly using ndarrays (although people still tend to believe that overall it's preferable to stick to the latter). One such example is matrix power:
mat ** 3 is the proper third matrix power of a matrix (whereas it's the elementwise cube of an ndarray). Unfortunately
numpy.linalg.matrix_power is quite more verbose. Furthermore, in-place matrix multiplication only works fine for the matrix class. In contrast, while both PEP 465 and the python grammar allow
@= as an augmented assignment with matmul, this is not implemented for ndarrays as of numpy 1.15.
Considering the above complications concerning the
matrix class there have been recurring discussions of its possible deprecation for a long time. The introduction of the
@ infix operator which was a huge prerequisite for this process happened in September 2015. Unfortunately the advantages of the matrix class in earlier days meant that its use spread wide. There are libraries that depend on the matrix class (one of the most important dependent is
scipy.sparse which uses both
numpy.matrix semantics and often returns matrices when densifying), so fully deprecating them has always been problematic.
Already in a numpy mailing list thread from 2009 I found remarks such as
numpy was designed for general purpose computational needs, not any one branch of math. nd-arrays are very useful for lots of things. In contrast, Matlab, for instance, was originally designed to be an easy front-end to linear algebra package. Personally, when I used Matlab, I found that very awkward -- I was usually writing 100s of lines of code that had nothing to do with linear algebra, for every few lines that actually did matrix math. So I much prefer numpy's way -- the linear algebra lines of code are longer an more awkward, but the rest is much better.
The Matrix class is the exception to this: is was written to provide a natural way to express linear algebra. However, things get a bit tricky when you mix matrices and arrays, and even when sticking with matrices there are confusions and limitations -- how do you express a row vs a column vector? what do you get when you iterate over a matrix? etc.
There has been a bunch of discussion about these issues, a lot of good ideas, a little bit of consensus about how to improve it, but no one with the skill to do it has enough motivation to do it.
These reflect the benefits and difficulties arising from the matrix class. The earliest suggestion for deprecation I could find is from 2008, although partly motivated by unintuitive behaviour that has changed since (in particular, slicing and iterating over a matrix will result in (row) matrices as one would most likely expect). The suggestion showed both that this is a highly controversial subject and that infix operators for matrix multiplication are crucial.
The next mention I could find is from 2014 which turned out to be a very fruitful thread. The ensuing discussion raises the question of handling numpy subclasses in general, which general theme is still very much on the table. There is also strong criticism:
What sparked this discussion (on Github) is that it is not possible to write duck-typed code that works correctly for:
- scipy.sparse sparse matrixes
The semantics of all three are different; scipy.sparse is somewhere between matrices and ndarrays with some things working randomly like matrices and others not.
With some hyberbole added, one could say that from the developer point of view, np.matrix is doing and has already done evil just by existing, by messing up the unstated rules of ndarray semantics in Python.
followed by a lot of valuable discussion of the possible futures for matrices. Even with no
@ operator at the time there is a lot of thought given to the deprecation of the matrix class and how it might affect users downstream. As far as I can tell this discussion has directly led to the inception of PEP 465 introducing matmul.
In my opinion, a "fixed" version of np.matrix should (1) not be a np.ndarray subclass and (2) exist in a third party library not numpy itself.
I don't think it's really feasible to fix np.matrix in its current state as an ndarray subclass, but even a fixed matrix class doesn't really belong in numpy itself, which has too long release cycles and compatibility guarantees for experimentation -- not to mention that the mere existence of the matrix class in numpy leads new users astray.
Eventually, first action to deprecate
numpy.matrix was taken in late November 2017. Regarding dependents of the class:
How would the community handle the scipy.sparse matrix subclasses? These are still in common use.
They're not going anywhere for quite a while (until the sparse ndarrays materialize at least). Hence np.matrix needs to be moved, not deleted.
while I want to get rid of np.matrix as much as anyone, doing that anytime soon would be really disruptive.
There are tons of little scripts out there written by people who didn't know better; we do want them to learn not to use np.matrix but breaking all their scripts is a painful way to do that
There are major projects like scikit-learn that simply have no alternative to using np.matrix, because of scipy.sparse.
So I think the way forward is something like:
Now or whenever someone gets together a PR: issue a PendingDeprecationWarning in np.matrix.__init__ (unless it kills performance for scikit-learn and friends), and put a big warning box at the top of the docs. The idea here is to not actually break anyone's code, but start to get out the message that we definitely don't think anyone should use this if they have any alternative.
After there's an alternative to scipy.sparse: ramp up the warnings, possibly all the way to FutureWarning so that existing scripts don't break but they do get noisy warnings
Eventually, if we think it will reduce maintenance costs: split it into a subpackage
It is no longer recommended to use this class, even for linear algebra. Instead use regular arrays. The class may be removed in the future.
And at the same time a
PendingDeprecationWarning has been added to
matrix.__new__. Unfortunately, deprecation warnings are (almost always) silenced by default, so most end-users of numpy will not see this strong hint.
Finally, the numpy roadmap as of November 2018 mentions multiple related topics as one of the "tasks and features [the numpy community] will be investing resources in":
Some things inside NumPy do not actually match the Scope of NumPy.
- A backend system for numpy.fft (so that e.g. fft-mkl doesn’t need to monkeypatch numpy)
- Rewrite masked arrays to not be a ndarray subclass – maybe in a separate project?
- MaskedArray as a duck-array type, and/or
- dtypes that support missing values
- Write a strategy on how to deal with overlap between numpy and scipy for linalg and fft (and implement it).
- Deprecate np.matrix
It's likely that this state will stay as long as larger libraries/many users (and in particular
scipy.sparse) rely on the matrix class. However, there's ongoing discussion to move
scipy.sparse to depend on something else, such as
pydata/sparse. Irrespective of the developments of the deprecation process users should use the
ndarray class in new code and preferably port older code if possible. Eventually the matrix class will probably end up in a separate package to remove some of the burdens caused by its existence in its current form.