In this paper, we propose a novel three-dimensional (3D) near-field beamforming (BF) design for Large Intelligent Surface (LIS). We firstly investigate the definitions of near-field and far-field of LIS, and derive the Fresnel near-field region where amplitudes variations are negligible and only phase variations worsen the harvested array-gains. We show that the Fresnel region which covers the majority part of near-field, can be enlarged by a factor of four when considering possible imperfectness from a conventional two-dimensional (2D) far-field BF. Therefore, it is of interest to design analog 3D-BF than can recover array-gain losses in this region. Secondly, with the proved decomposition theorem we show that the optimal 3D-BF can be decomposed into a 2D far-field BF and a one-dimensional (1D) near-field BF. The 2D far-field BF compensates phase variations from mismatches in the azimuth and elevation angles, while the 1D near-field BF compensates remaining phases variations caused by distance differences from the UE to different antenna-elements on LIS. Such a proposed “2D+1D” BF design reduces codebook size significantly and is fully compatible with the existing far-field BF in the fifth-generation new-radio (5G-NR) system. Thirdly, we analyze an optimal codebook design for the 1D near-field BF, and show that with a small codebook it can perform close to optimal. Numerical results further verify that our proposal is effective and robust to recover array-gains in the near-field of LIS.