Motivated by the development of non-intrusive methods for high dimensional parametric PDE's, we study the stability of a sparse high dimensional polynomial interpolation procedure introduced in [6]. A key aspect of this procedure is its hierarchical structure: the sampling set is progressively enriched together with the polynomial space. The evaluation points are selected from a grid obtained by tensorization of a univariate sequence. The Lebesgue constant that quantifies the stability of the resulting interpolation operator depends on the choice of this sequence. Here we study the -Leja sequence, obtained by the projection of Leja sequences on the complex unit circle, with initial value 1, onto [−1, 1]. For this sequence, we prove cubic growth in the number of points for the Lebesgue constant of the multivariate interpolation operator, independently of the number of variable and of the shape of the polynomial space.