We study the problem of maximum-likelihood (ML) estimation of an approximate common eigenstructure, i.e. an approximate common eigenvectors set (CES), for an ensemble of covariance matrices given a collection of their associated i.i.d vector realizations. This problem has a direct application in multi-user MIMO communications, where the base station (BS) has access to instantaneous user channel vectors through pilot transmission and attempts to perform joint multi-user Downlink (DL) precoding. It is widely accepted that an efficient implementation of this task hinges upon an appropriate design of a set of common "virtual beams", that captures the common eigenstructure among the user channel covariances. In this paper, we propose a novel method for obtaining this common eigenstructure by casting it as an ML estimation problem. We prove that in the special case where the covariances are jointly diagonalizable, the global optimal solution of the proposed ML problem coincides with the common eigenstructure. Then we propose a projected gradient descent (PGD) method to solve the ML optimization problem over the manifold of unitary matrices and prove its convergence to a stationary point. Through exhaustive simulations, we illustrate that in the case of jointly diagonalizable covariances, our proposed method converges to the exact CES. Also, in the general case where the covariances are not jointly diagonalizable, it yields a solution that approximately diagonalizes all covariances. Besides, the empirical results show that our proposed method outperforms the well-known joint approximate diagonalization of eigenmatrices (JADE) method in the literature.