Commonly, variance-based global sensitivity analysis methods are popular and applicable to quantify the impact of a set of input variables on output response. However, for many engineering practical problems, the output response is not single but multiple, which makes some traditional sensitivity analysis methods difficult or unsuitable. Therefore, a novel global sensitivity analysis method is presented to evaluate the importance of multi-input variables to multi-output responses. First, all functions of characterizing the relationship of input variables and output responses are decomposed based on the high dimensional model representation (HDMR). For each decomposed function, all its function sub-terms are orthogonal each other. Then, the covariance matrices are obtained by calculating the variances of each HDMR and the covariances between HDMRs. According to the integral condition of HDMR, the total covariance matrix of output responses is equal to the sum of all the covariance matrices of their corresponding function sub-terms for each HDMR. Further, a coefficient K is defined to makes that the determinant of the total covariance matrix is equal to the determinant of the covariance matrices for each function sub-terms. Based on the ratio of sub-determinant and K times total determinant, the first- or high-order sensitivity index can be calculated to evaluate the effect of first- or high-order function sub-term, respectively. The sensitivity of variable X_i can be obtained by calculating the sum of the first- and high-order sensitivity indices whose function sub-terms include the variable X_i at the same time. Finally, several numerical and engineering examples are employed to demonstrated the accuracy and practicality of the presented global sensitivity analysis method on multi-input variables and multi-output responses.