Regional seismic fragility assessment of bridge portfolios must address the embedded uncertainties and variations stemming from both the earthquake hazard and bridge attributes (e.g., geometry, material, design detail). To achieve bridge‐specific fragility assessment, multivariate probabilistic seismic demand models (PSDM) have recently been developed that use both the ground motion intensity measure and bridge parameters as inputs. However, explicitly utilizing bridge parameters as inputs requires numerous nonlinear response history analyses (NRHAs). In this situation, the associated computational cost increases exponentially for high‐fidelity bridge models with complex component connectivity and sophisticated material constitutive laws. Moreover, it remains unclear how many analyses are sufficient for the response data and the resulting demand model to cover the entire solution space without overfitting. To deal with these issues, this study integrates Gaussian process regression (GPR) and active learning (AL) into a multistep workflow to achieve efficient regional seismic fragility assessment of bridge portfolios. The GPR relaxes the probability distribution assumptions made in typical cloud analysis‐based PSDMs to enable heteroskedastic nonparametric seismic demand modeling. The AL leverages the varying standard deviation to select the least but most representative bridge‐model‐ground‐motion sample pairs to conduct NRHA with much‐improved efficiency. Both independent and correlated multi‐output GPRs are proposed to deal with bridge portfolios with seismic demand correlations among multiple components (column, bearing, shear key, abutment, unseating, and joint seal). Considering a single benchmark highway bridge class in California as the case study, the AL‐GPR framework and the associated component‐level fragility results are investigated in terms of their efficiency, accuracy, and robustness. The fragility results show that 70 AL‐selected samples would enable the GPR to derive bridge‐specific fragility models comparable to the ones using the multiple stripes analysis approach with 1950 ground motions considered for each individual bridge. The AL‐GPR model also successfully captures the physics of how bridge span length, deck area, column slenderness, and steel reinforcement ratio would change the damage state exceedance probabilities of different bridge components. The efficiency of AL stems from the fact that, with the multi‐output independent GPR, a stable and reliable fragility model can be achieved using 50 AL‐selected samples compared to at least 270 randomly chosen samples. The proposed methodology advances the state of the art in enabling more efficient and reliable regional seismic fragility assessment of multi‐component bridge portfolios.