As touched upon in in #13797

matrix was a collection of behaviors, some good and some bad:

• .I and .H attributes. These were handy.
• * means matmul - not inherently bad, but perhaps surprising, and certainly less necessary in the presence of @. I think there might still be some value in allowing "multiply (*)" on matrix objects to mean what @ does on array objects - it allows them to be dropped into some equations in place of scalars
• Everything is always exactly 2d. This was the cause of almost all the pain within numpy when handling subclasses - it broke liskov substitution

I think there might be a good use case for a matrix with some subset of the following rules:

• Construction fails if the input array is not ndim >= 2
• Helpers for (stacks of) row and column vectors:
  • matrix.rowvec(rows) == matrix(rows[..., None, :])
  • matrix.colvec(cols) == matrix(cols[..., :, None])
• Simpler indexing rules:
  • Shape of the returned array always matches ndarray
  • m[..., i, j], where i and j index and preserve the last two dimensions, returns a matrix object. Any type of indexing that discards the last two dimensions downcasts to array
• .T means linalg.transpose(self), .H means .T.conj()
• .I either means linalg.inverse(self) or perhaps an object with __mul__ = linalg.solve

Such as class would not be compatible with np.matrix, but it might be a good candidate for a new numpy/matrix package.