We present rigorous modeling and optimal design for 3D graded index photonic lens at the telecommunication wavelnegth. Based on our numerical results, the efficiency of the optimized designs can reach 87%. In this work, we use a global optimization method based on statistical learning in order to enhance the efficiency of a graded index photonic metalens (see Fig. 1(a)). This method belongs to the class of Bayesian optimization methods and is known as Efficient Global Optimization (EGO) [1, 2]. As a 3D fullwave solver, we use our rigorous Discontinuous Galerkin Time Domain (DGTD) solver from the DIOGENES software suite [3] in order to rigorously model such configuration. For the metalens Fig. 1. (a): schematic view of the 3D photonic metalens. The structure consists of Si region (green part) on top of subtrate made of SiO 2 (red part). The Si region is divided into three parts; the input port with width W = 3000 nm, the output port with width w = 300 nm, and in between we have several Si strips with length L and height h = 310 nm. The input mode is injected from the surface S in (solution of the 2D modal problem) and the objective function is computed at the output surface S out as the overlap between the obtained solution and the solution to the 2D modal problem at this surface. The widths of the strips are denoted by e i , i ∈ {0, 6}. (b) and (c) represent the e(H y) and the e(E x), of the optimized design, respectively. presented in Fig. 1(a), we optimize 8 parameters, which are the widths of the strips and their common length L. In addition we focus on the most challenging case, i.e., TE case (unlike other works in the literature, where the classical TM case is considered), where the input field is polarized perpendicular to the strips. In this case, the near field coupling has to be taken into account and a higher order fullwave solver is needed. The optimization results reveal that at least two global designs (different parameters) have been obtained in which the efficiency reaches 80%. Furthermore, we show also that the efficiency may reach 87% in case the lengths of the strips are assumed to be different from each others in our optimization process.