Laser shock hole-clinching is one of mechanical joining processes in which metal foils are joined together based on plastic deformation caused by laser-induced shock wave. In the process, fracture often occurs due to the unreasonable arrangement of laser process parameters, suggesting that it is urgent to seek a feasible way to improve the joint forming quality. However, it consumes a great deal of time and effort to obtain the optimal solution of process parameters from so many different combinations through experiments and numerical simulations. In this study, the mathematical models between laser process parameters and joint forming quality were established through response surface methodology (RSM). A finite element analysis model of laser shock hole-clinching process for pure copper and pre-pierced stainless steel foils was developed, and then it was used to perform the calculation scheme arranged by design of experiments. Analysis of variance (ANOVA) was implemented to evaluate the statistical significance of RSM models and the influence of process parameters on objectives. Multi-objective optimization was carried out to achieve the optimal combination of laser process parameters by using genetic algorithm (GA). It is revealed that the RSM-GA-integrated approach is an effective way to realize the modeling and optimization of laser process parameters for laser shock hole-clinching. The pulsed laser energy (E), number of laser pulses (N), and laser spot diameter (D) are statistically significant, and both the interlock value and the maximum thinning rate are sensitive to these parameters based on ANOVA. The fitting precision analysis indicates that the established RSM models can be used to navigate the design space of variables and predict the actual data with high accuracy. Moreover, it is found that the influence order of laser process parameters on the interlock value from strong to weak is N, E, and D, while this order changes to D, N, and E for the maximum thinning rate. According to the Pareto noninferior solutions and the corresponding values of satisfaction function, the optimal combination is E = 164 mJ, N = 25, and D = 2.1 mm in the given design space. The GA optimization result has been experimentally confirmed.