The substitution principle is interesting (and, being trained as an astronomer, not a computer scientist, I had not heard of it before). I think matrix is indeed obviously wrong here (with indexing being more annoying, but multiplication being a good example as well).

Perhaps more interesting as an example to consider is MaskedArray, which is much closer to a sensible subclass, though different from Quantity in that what is masked can itself be an ndarray subclass. In a sense, it is more of a container class, in which the operations are done on what is inside it, with some care taken about which elements are fixed. This becomes quite clear when one thinks of implementing __array_ufunc__ or __array_function__: for Quantity, calling super after dealing with the units is very logical, for MaskedArray, it makes more sense to call the (universal) function again on the contents [1].

For this particular class, if reimplemented, it might make most sense as a "mixin" since its attributes depend both on the masked class (.mask, etc.) and on what is being masked (say, .unit for a quantity). Thus, the final class might be an auto-generated new class (e.g., MaskedQuantity(MaskedArray, Quantity)). We have just added a new Distribution class to astropy which is based on this idea [2] (since this uses casting from structured dtypes which hold the samples to real arrays on which functions are evaluated, this probably could be done just as well or better with more flexible dtypes, but we have to deal with what's available in the real world, not the ideal one...).