In additional to the phonon (massless Goldstone mode) in Galilean invariant superfluid, there is another type of mode known as the Higgs amplitude mode in superfluid with emergent Lorentz invariance. In two dimensions, due to the strong decay into phonons, whether this Higgs mode is a detectable excitation with sharp linear response has been controversial for decades. Recent progress gives a positive answer to this question. Here, we review a series of numerical studies of the linear response of a two-dimensional Lorentz invariant superfluid near the superfluid-Mott insulator quantum critical point (SF-MI QCP). Particularly, we introduce a numerical procedure to unbiasedly calculate the linear response properties of strongly correlated systems. The numerical procedure contains two crucial steps, i.e., one is to use a highly efficient quantum Monte Carlo method, the worm algorithm in the imaginary-time path-integral representation, to calculate the imaginary time correlation functions for the system in equilibrium; and then, the other is, based on the imaginary time correlation functions, to use the numerical analytical continuation method for obtaining the real-time (real-frequency) linear response function. Applying this numerical procedure to the two-dimensional Bose Hubbard model near SF-MI QCP, it is found that despite strong damping, the Higgs boson survives as a prominent resonance peak in the kinetic energy response function. Further investigations also suggest a similar but less prominent resonance peak near SF-MI QCP on the MI side, and even on the normal liquid side. Since SF-MI quantum criticality can be realized by ultracold aotms in optical lattice, the Higgs resonance peak can be directly observed in experiment. In addition, we point out that the same Higgs resonance peak exists in all quantum critical systems with the same universality, namely (2 + 1)-dimensional relativistic U(1) criticality, as SF-MI QCP.