We present a method for constructing multivariate refinable Hermite interpolants and their associated subdivision algorithms based on a combination of analytical and numerical approaches. Being the limit of a linear iterative procedure, the critical $L^2$ Sobolev smoothness of a refinable Hermite interpolant is given by the spectral radius of a matrix dependent upon the {\it refinement mask}. The design question: given certain constraints (support size, symmetry type, refinement pattern etc.), how can one choose the refinement mask so that the resulting refinable function has optimal smoothness? This question naturally gives rise to a spectral radius optimization problem. In general, the objective function is not convex, and may not be differentiable, or even Lipschitz, at a local minimizer. Nonetheless, a recently developed robust solver for nonsmooth optimization problems may be applied to find local minimizers of the spectral radius objective function. In fact, we find that in specific cases that are of particular interest in the present context, the objective function is smooth at local minimizers and may be accurately minimized by standard techniques. We present two necessary mathematical tricks that make the method practical: (i) compression of matrix size based on symmetry; (ii) efficient computation of gradients of the objective function. We conclude by reporting some computational results.