Similarly to what was recently done for rotate(). (skew is rather rarely used, but given that @timhoffm suggested this in #22108 (review), I may as well do it too :-))

To address @timhoffm's other comment ("Does it still make sense to keep Affine2D._mtx as a numpy array if we only do element-wise operations?"): I'm quite convinced the whole transform stack would be faster if the transformation matrix was not a numpy array (because 3x3 is a size where the numpy's overhead is generally bigger than the gains from vectorization), but changing everything at once (even better would be to move things to C, but using plain C structs (or equivalently C++ objects) to store the coefficients) would be quite a big PR. So I'm doing the easy parts first :-)

PR Summary

PR Checklist

Tests and Styling

• Has pytest style unit tests (and pytest passes).
• Is Flake 8 compliant (install flake8-docstrings and run flake8 --docstring-convention=all).

Documentation

• New features are documented, with examples if plot related.
• New features have an entry in doc/users/next_whats_new/ (follow instructions in README.rst there).
• API changes documented in doc/api/next_api_changes/ (follow instructions in README.rst there).
• Documentation is sphinx and numpydoc compliant (the docs should build without error).