The most important factor in a multivariate kernel density estimation is a~choice of a bandwidth matrix. Because of its role in controlling both the amount and the direction of multivariate smoothing, this choice is a particularly important. Considerable attention has been paid to constrained parameterization of the bandwidth matrix such as a diagonal matrix or pre-transformation of the data. General multivariate kernel density derivative estimators has been investigated in paper Chac\'on, Test, p. 375--398, Vol. 19, 2011. The present paper is focused on data-driven selectors of full bandwidth matrices for a density and its gradient. This method is based on an optimally balanced relation between integrated variance and integrated squared bias. The analysis of statistical properties shows the rationale of the proposed method. It is also given the relative rate of convergence to compare the method with cross-validation and plug-in methods. The utility of this method is illustrated through a~simulation study and application to real data.