We present a comprehensive study of the velocity operator, \hat{\boldsymbol{v}}=\frac{i}{\hbar} [\hat{H},\hat{\boldsymbol{r}}]\,𝐯̂=iℏ[Ĥ,𝐫̂], when used in crystalline solids calculations. The velocity operator is key to the evaluation of a number of physical properties and its computation, both from a practical and fundamental perspective, has been a long-standing debate for decades. Our work summarizes the different approaches found in the literature, but never connected before in a comprehensive manner. In particular we show how one can compute the velocity matrix elements following two different routes. One where the commutator is explicitly used and another one where the commutator is avoided by relying on the Berry connection. We work out an expression in the latter scheme to compute velocity matrix elements, generalizing previous results. In addition, we show how this procedure avoids ambiguous mathematical steps and how to properly deal with the two popular gauge choices that coexist in the literature. As an illustration of all this, we present several examples using tight-binding models and local density functional theory calculations, in particular using Gaussian-type localized orbitals as basis sets. We show how the the velocity operator cannot be approximated, in general, by the k-gradient of the Bloch Hamiltonian matrix when a non-orthonormal basis set is used. Finally, we also compare with its real-space evaluation through the identification with the canonical momentum operator when possible. This comparison offers us, in addition, a glimpse of the importance of non-local corrections, which may invalidate the naive momentum-velocity correspondence.