Fast computation of large scale marginal extremes with multi-dimensional covariates