The frequency domain is allegedly preferred because it replaces convolution of complexity $n^2$ with a diagonal matrix multiplication. Yet, I see that in z-domain we have multiplication of polynomials, $(a_0z+a_1z^{-1}+a_2z^{-2}+\ldots)(b_0z+b_1z^{-1}+b_2z^{-2}+\ldots)$, which is, despite its name, is nothing more than convolution with $n^2$ operations: you need multiply one polynomial with every term of the other and combine the terms. This looks like square convolution matrix rather than diagonal one. What is the benefit?
In this regards, I recalled my old question: why nobody speaks about eigenfunctions of the s- and z-domains? Everybody speaks that transfer function is linear and diagonal in the frequency (s- and z-domains). But, when you speak about linear operators, you must bring up their natural basis, especially when speak about diagonalization. How do eigenvectors look like in the Laplace and z-domain? What is the basis of this domain? Might be you do not need to prove the convolution theorem then?