We consider the smoothness of solutions of a system of refinement equations written in the form phi = Sigma(alpha is an element of z) a(alpha)phi(2.-alpha), where the vector of functions phi = (phi(1),...,phi(r))(T) is in (L-p( R))(r) and a is a finitely supported sequence of r x r matrices called the refinement mask. We use the generalized Lipschitz space Lip*(nu, L-p(R)), nu >0, to measure smoothness of a given function. Our method is to relate the optimal smoothness, nu(p)(phi), to the p-norm joint spectral radius of the block matrices A(epsilon), epsilon = 0,1, given by A(epsilon) = (a(epsilon + 2 alpha - beta))(alpha,beta), when restricted to a certain finite dimensional common invariant subspace V. Denoting the p-norm joint spectral radius by rho(p)(A(0)\(V), A(1)\(V)), we show that nu(p)(phi) greater than or equal to 1/p - log(2)rho(p)(A(0)\(V), A(1)\(V)) with equality when the shifts of phi(1),...,phi(r) are stable and the invariant subspace is generated by certain vectors induced by difference operators of sufficiently high order. This allows an effective use of matrix theory. Also the computational implementation of our method is simple. When p = 2, the optimal smoothness is also given in terms of the spectral radius of the transition matrix associated with the refinement mask. To illustrate the theory, we give a detailed analysis of two examples where the optimal smoothness can be given explicitly. We also apply our methods to the smoothness analysis of multiple wavelets. These examples clearly demonstrate the applicability and practical power of our approach.