I want to start with misconceptions:
Modern GPUs (NVIDIA for quite a while, and AMD since Southern Islands) do not meaningfully support vector/matrix operations natively in hardware. They are vector architectures in a different direction: each component of a vector (x, y, z) are generally 32- or 64-valued, containing values for each element in a lane. So a 3D dot product is not usually an instruction, it is a multiply and two multiply-adds.
Additionally, counting primitive operations like multiply-add, transforming a vector by a quaternion is more expensive than transforming a vector by a matrix. Transforming a vector by a 3x3 matrix is 3 multiplies and 6 multiply-adds, and transforming a vector by a quaternion is two quaternion multiplies, each of which consist of 4 multiplies and 12 multiply-adds. (You can get less naïve than this—here's a writeup on a faster way—but it's still not as cheap as multiplying a vector by a matrix.)
However, performance is not always determined simply by counting the number of ALU operations it performs. Quaternions require less space than the equivalent matrix (assuming you are only doing pure rotation/scale), and that means less storage space and less memory traffic. This is often important in animation (which is conveniently also often where the nice interpolation properties of quaternions show up).
Other than that:
- Matrices use more space because they support more operations. A 3x3 matrix can contain nonuniform scale, skew, reflection, and orthogonal projection.
- Matrices can be naturally thought of as basis vectors, and easily constructed from those vectors.
- Multiplying one quaternion by another (composing two rotations) is less operations than multiplying one matrix by another.