How to batch-mask (multiply) multiband raster with mask of different cell size?

I have a multiband GeoTIFF raster of certain cell size and a binary GeoTIFF mask (1 or no-data, but can be recreated differently) with different cell size. I need to mask the raster such that only pixels "under" mask pixels are data (others thus become no-data). Both the raster and mask are same EPSG:4326 and overlay, the mask’s extent is only a fraction of the raster.

GDAL (Python) attempt

gdal_calc.py, which does not work, since the raster happen to have different cell sizes.

Rasterio attempt

I tried masking with rasterio, which requires raster to have same extent. I thus first converted the mask extent to Shapely.geometry.box, clipped the raster but bumped to the same problem as in GDAL.

QGIS attempt

I tried multiplying the raster and mask in QGIS native Raster Calculator, which works (despite the variable cell size), but is difficult to batch, since the raster is multiband. I do not properly know PyQGIS yet, but could adopt and tinker with a solution.