The Active Subspace (AS) method is a widely used technique for identifying the most influential directions in high-dimensional input spaces that affect the output of a computational model. The standard AS algorithm requires a sufficient number of gradient evaluations (samples) of the input output map to achieve quasi-optimal reconstruction of the active subspace, which can lead to a significant computational cost if the samples include numerical discretization errors which have to be kept sufficiently small. To address this issue, we propose a multilevel version of the AS method (MLAS) that utilizes samples computed with different accuracies and yields different active subspaces across accuracy levels, which can match the accuracy of single-level AS with reduced computational cost, making it suitable for downstream tasks such as function approximation. In particular, we propose to perform the latter via optimally-weighted least-squares polynomial approximation in the different active subspaces, and we present an adaptive algorithm to choose dynamically the dimensions of the active subspaces and polynomial spaces. We demonstrate the practical viability of the MLAS method with polynomial approximation through numerical experiments based on random partial differential equations (PDEs).