In order to effectively solve the low precision problem of the single gravity density inversion and the magnetic susceptibility inversion, and the limitation of the gravity-magnetic joint inversion method based on the petrophysical parameter constraint, this paper studies the three-dimensional gravity-magnetic cross-gradient joint inversion based on the structural coupling and the fast optimization algorithm. Based on the forward and inversion modeling of three-dimensional gravity density and three-dimensional magnetic susceptibility using the same underground grid, along with cross-gradient coupling as the structural constraint, we propose a new gravity-magnetic joint inversion objective function including the data fitting term, the total variation regularization constraint term and the cross-gradient term induced by the structural coupling. The depth weighted constraint and the data weighting constraint are included into the objective function, which requires different physical property models to minimize their respective data residuals. At the same time, the cross-gradient term tends to zero, so that the structure of the gravity and magnetic models tends to be consistent. In realization, we address a fast and efficient gradient algorithm to iteratively solve the objective function. We apply this new joint inversion algorithm to the 3D gravity-magnetic model inversion test and compare it with the results of a single inversion algorithm. The experimental tests of synthetic data indicate that the gravity-magnetic cross-gradient joint inversion method can effectively improve the accuracy of the anomaly position and numerical accuracy of the inverted anomaly physical parameters compared with the single physical inversion method.