Abstract: Direction-of-arrival (DOA) estimation is a fundamental technique in array signal processing. DOA estimation methods based on programmable metasurfaces enable single-channel reception architectures, significantly reducing system cost and hardware complexity. However, most existing approaches rely on spatial discretization, and their performance is severely degraded by grid mismatch errors. To address this issue, this paper proposes an off-grid DOA estimation method for programmable metasurfaces based on joint sparse recovery. By employing a first-order Taylor expansion, a joint overcomplete dictionary incorporating both grid terms and derivative terms is constructed, enabling explicit modeling and correction of off-grid errors. Furthermore, to reduce the computational and storage complexity of the joint model in two-dimensional scenarios, a Kronecker-decomposition-based model simplification strategy is introduced, transforming high-dimensional operations into low-dimensional matrix computations. Based on the formulated models, a joint orthogonal matching pursuit (JOMP) algorithm and its Kronecker-based efficient implementation (KJOMP) are developed. Simulation and experimental results demonstrate that the proposed JOMP algorithm achieves significantly improved DOA estimation accuracy in both one-dimensional and two-dimensional scenarios under various metasurface sizes and signal-to-noise ratios. While maintaining estimation accuracy comparable to that of the JOMP algorithm, the KJOMP algorithm achieves significant reductions in both computational time and spatial complexity. The proposed method is well suited for applications in radar, Internet of Things (IoT), sixth-generation (6G) communications, positioning, navigation, and unmanned aerial vehicles.