Discovering underlying partial differential equations (PDEs) from observational data has important implications across fields. It bridges the gap between theory and observation, enhancing our understanding of complex systems in applications. In this paper, we propose a novel approach, termed physics-informed sparse optimization (PIS), for learning surface PDEs. Our approach incorporates both $L_2$ physics-informed model loss and $L_1$ regularization penalty terms in the loss function, enabling the identification of specific physical terms within the surface PDEs. The unknown function and the differential operators on surfaces are approximated by some extrinsic meshless methods. We provide practical demonstrations of the algorithms including linear and nonlinear systems. The numerical experiments on spheres and various other surfaces demonstrate the effectiveness of the proposed approach in simultaneously achieving precise solution prediction and identification of unknown PDEs.