Recent research found that multivariate shrinkage on multiwavelet transform coefficients further improves the traditional wavelet methods. It is because the multiwavelet transform, with appropriate initialization, provides better representation of signals so that their difference from noise can be clearly identified. In this paper, we consider the optimal threshold selection for multiwavelet denoising by using a multivariate shrinkage function. Firstly, we study the threshold selection using the Stein's unbiased risk estimator (SURE) for each resolution level when the noise structure is given. Then, we consider the method of generalized cross validation (GCV) when the noise structure is not known a priori. Simulation results show that the higher multiplicity (>2) wavelets usually give better denoising results. Besides, the proposed threshold estimators often suggest better thresholds as compared with the traditional estimators.