Remove usage of np.matrix

@eteq

Numpy matrices are quite notorious for, in hindsight, being a relatively bad design. One thing that is particularly annoying is that one cannot have stacks of matrices, which leads to work-around like those in #5085 -- https://github.com/astropy/astropy/pull/5085/files#diff-11cd3ae0a0afb774b06628d3912308a0R36

Currently, we only use np.matrix in a few places. I suggest we convert all of these to regular arrays so that they can be properly vectorized. I think none are exposed, so we should just be able to make the changes, but I thought I would get opinions before proceeding cc @eteq, @astrofrog, @StuartLittlefair.

p.s. A bit off-topic, but in many places in coordinates, we really need cartesian coordinates with x,y,z as last dimension rather than first. We could make an alternative property on CartesianCoordinate (say, get_xyz(self, xyz_axis=0) or to_vector) or perhaps even implement a Vector class (which CartesianCoordinate could be a subclass of). The slight relevance of the latter is that if we go that route, we obviously could also have a good Matrix class...

mhvk@swan:~$ find packages/astropy -name \*.py -exec grep np.matrix {} \; -print        X = np.matrix(np.vstack([x, np.ones(len(x))]).T)
        beta = np.linalg.inv(X.T * X) * X.T * np.matrix(y).T
packages/astropy/astropy/modeling/tests/test_fitters.py
            cdelt = np.matrix(np.diag(self.wcs.get_cdelt()))
            pc = np.matrix(self.wcs.get_pc())
                cdelt = np.matrix(self.wcs.cd) * np.matrix(np.diag(self.wcs.cdelt))
                cdelt = np.matrix(np.diag(self.wcs.cdelt))
                pc = np.matrix(self.wcs.pc)
packages/astropy/astropy/wcs/wcs.py
        return np.matrix(((c, s, 0),
        return np.matrix(((c, 0, -s),
        return np.matrix(((1, 0, 0),
        return R.view(np.matrix)
packages/astropy/astropy/coordinates/angles.py

For reference, multi-D versions of angles.rotation_matrix and angle_axis work if one uses an ndarray subclass which overrides __mul__ with __matmul__; see https://github.com/astropy/astropy/compare/master...mhvk:coord-no-np-matrix?expand=1

In turn, that means these can work over the whole coordinate code if the places where matrices are multiplied are adjusted to make those use np.matmul (or a suitably backported version for older numpy).

👍 to dropping np.matrix in favor of matmul. I'm less sure about some of the more subtle consequences of the off-topic bit... (E.g., #4926 (comment))

And if I understand right, this will prepare us to have an all-@ codebase whenever we decide only py 3.5 will be supported ;)

@eteq

@eteq, @StuartLittlefair - as proof is often in the pudding, I made two simple PRs that try either replacing np.matrix with a new AstropyMatrix that can be multidimensional, or by using matmul instead of multiplication where needed. I'm slightly leaning towards the latter, in particular because the erfa matrix routines always return regular arrays, so one might as well be similar to that.

p.s. From a quick timing tests, it seems that even np.einsum('...ij,...jk->...ik', m1, m2) is faster than m1 * m2 (where m1 and m2 are matrices). So even just for performance reasons these are wins...

@mhvk

Very nice @mhvk.

I just admit I'm torn by the two PRs. The former leads to much cleaner and readable code. However, it shares a fundamental flaw with the matrix class, which PEP465 puts better than I can:

Code that expects an ndarray and gets a matrix, or vice-versa, may crash or return incorrect results. Keeping track of which functions expect which types as inputs, and return which types as outputs, and then converting back and forth all the time, is incredibly cumbersome and impossible to get right at any scale. Functions that defensively try to handle both types as input and DTRT, find themselves floundering into a swamp of isinstance and if statements.

For this reason I think I prefer the second PR.

@eteq

I think the most basic mistake of the np.matrix class (as well as my quick AstropyMatrix class) was that the shape includes the shape of the matrix itself. I think much more logically it would just be the shape of the matrix stack. This would be similar for a vector (think CartesianRepresentation, where the shape reflects the shape of the coordinates, not also that there are three parts).

But even so, I'm coming around to seeing the use of regular arrays and matmul as the fastest and clearest path for now. And of course, as @eteq pointed out, once we only support python >=3.5, we can swap out all those matmul with @ and then the notation looks OK again... (In the meantime, maybe I should create a new function that just calculates the products of a whole slew of matrices, as
matmul(matmul(matmul(a, b), c), d) really does look ugly; matrix_product(a, b, c, d)?

fixed by #5104.

@mhvk

@mhvk - Use of np.matrix in Astropy now shows up as a warning in the logs of affiliated package tests, e.g. for reproject:

software/anaconda3/envs/astropy37/lib/python3.6/site-packages/reproject/tests/test_mosaicking.py::TestOptimalWCS::()::test_resolution
  /Users/deil/software/anaconda3/envs/astropy37/lib/python3.6/site-packages/astropy/wcs/wcs.py:3043: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    cdelt = np.matrix(np.diag(self.wcs.get_cdelt()))
  /Users/deil/software/anaconda3/envs/astropy37/lib/python3.6/site-packages/astropy/wcs/wcs.py:3044: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    pc = np.matrix(self.wcs.get_pc())
  /Users/deil/software/anaconda3/envs/astropy37/lib/python3.6/site-packages/numpy/matrixlib/defmatrix.py:68: PendingDeprecationWarning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.
    return matrix(data, dtype=dtype, copy=False)

Looks like this is the only remaining use of np.matrix in Astropy:

$ ack 'np\.matrix'
astropy/wcs/wcs.py
3048:            cdelt = np.matrix(np.diag(self.wcs.get_cdelt()))
3049:            pc = np.matrix(self.wcs.get_pc())
3053:                cdelt = np.matrix(self.wcs.cd) * np.matrix(np.diag(self.wcs.cdelt))
3055:                cdelt = np.matrix(np.diag(self.wcs.cdelt))
3058:                pc = np.matrix(self.wcs.pc)

astropy/modeling/tests/test_fitters.py
500:        X = np.matrix(np.vstack([x, np.ones(len(x))]).T)
501:        beta = np.linalg.inv(X.T * X) * X.T * np.matrix(y).T

@mhvk - Do you the use of np.matrix in astropy/wcs/wcs.py should also be replaced?

Actually, it seems to be only one method pixel_scale_matrix where it's used:
https://github.com/astropy/astropy/pull/2832/files#diff-d4c5b9e74c6c26000684075497033504R2741

Should be easy to rewrite using np.array and np.dot, no?

@cdeil

@cdeil I think we should discuss this in a new issue (or PR). The topic is related but the NumPy warning is relatively new, so I'm a bit hesitant to reopen this issue. But we shouldn't let that important comment linger in a closed issue.

Odd, I thought I had gotten them all! That's why I felt safe in introducing a warning in numpy... Anyway, it indeed shouldn't be too hard to write a PR that fixes this.

@mhvk

@mhvk - Do you have time to send a PR? I could also give it a try. Should at the end np.dot or np.matmul be called?

@cdeil

@cdeil - if you could tackle it, that would be great - somehow in my test environment pytest is not catching the pending deprecation warnings.

And, yes, it should just entail removing matrix and replacing the relevant * with @ (allowed in python 3.5! though perhaps using np.matmul is better for backporting...)

Navigation Menu

Repository navigation

Remove usage of `np.matrix` #5088

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Remove usage of np.matrix #5088

Description

Activity

mhvk commented on Jun 13, 2016

eteq commented on Jun 15, 2016

mhvk commented on Jun 15, 2016

StuartLittlefair commented on Jun 16, 2016

mhvk commented on Jun 16, 2016

mhvk commented on Jul 29, 2016

cdeil commented on Sep 2, 2018

MSeifert04 commented on Sep 2, 2018

mhvk commented on Sep 3, 2018

cdeil commented on Sep 3, 2018

mhvk commented on Sep 3, 2018

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

Remove usage of `np.matrix` #5088