Purpose: To improve the quality of quantitative MR images recovered from undersampled measurements, we incorporate the signal model of the variable-flip-angle (VFA) multi-echo 3D gradient-echo (GRE) method into the reconstruction of $T_1$, $T_2^*$ and proton density (PD) maps. Additionally, we investigate the use of complementary undersampling patterns to determine optimal undersampling schemes for quantitative MRI. Theory: We propose a probabilistic Bayesian formulation of the recovery problem. Our proposed approach, approximate message passing with built-in parameter estimation (AMP-PE), enables the joint recovery of distribution parameters, VFA multi-echo images, and $T_1$, $T_2^*$, and PD maps without the need for hyperparameter tuning. Methods: We conducted both retrospective and prospective undersampling to obtain Fourier measurements using variable-density and Poisson-disk patterns. We investigated a variety of undersampling schemes, adopting complementary patterns across different flip angles and/or echo times. Results: AMP-PE adopts a joint recovery strategy, it outperforms the state-of-the-art $l1$-norm minimization approach that follows a decoupled recovery strategy. For $T_1$ mapping, employing fixed sampling patterns across different echo times produced the best performance. Whereas for $T_2^*$ and proton density mappings, using complementary sampling patterns across different flip angles yielded the best performance. Conclusion: AMP-PE achieves better performance by combining information from both the MR signal model and the sparse prior on VFA multi-echo images. It is equipped with automatic and adaptive parameter estimation, and works naturally with the clinical prospective undersampling scheme.