PURPOSE: Diffusion MRI (dMRI) data typically suffer of significant cross-site variability, which prevents naively performing pooled analyses. To attenuate cross-site variability, harmonization methods such as the rotational invariant spherical harmonics (RISH) have been introduced to harmonize the dMRI data at the signal level. A common requirement of the RISH method is the availability of healthy individuals who are matched at the group level, which may not always be readily available, particularly retrospectively. In this work, we propose a framework to harmonize dMRI without matched training groups. METHODS: Our framework learns harmonization features while controlling for potential covariates using a voxel-based generalized linear model (GLM). RISH-GLM allows us to simultaneously harmonize data from any number of sites while also accounting for covariates of interest, thus not requiring matched training subjects. Additionally, RISH-GLM can harmonize data from multiple sites in a single step, whereas RISH is performed for each site independently. RESULTS: We considered data of training subjects from retrospective cohorts acquired with three different scanners and performed three harmonization experiments of increasing complexity. First, we demonstrate that RISH-GLM is equivalent to conventional RISH when trained with data of matched training subjects. Second, we demonstrate that RISH-GLM can effectively learn harmonization with two groups of highly unmatched subjects. Third, we evaluate the ability of RISH-GLM to simultaneously harmonize data from three different sites. CONCLUSION: RISH-GLM can learn cross-site harmonization both from matched and unmatched groups of training subjects and can effectively be used to harmonize data of multiple sites in one single step.