Several works in the past decade have used the ratio between total (rest 8-1000$mu$m) infrared and radio (rest 1.4~GHz) luminosity in star-forming galaxies (q$_{IR}$), often referred to as the infrared-radio correlation (IRRC), to calibrate radio emission as a star formation rate (SFR) indicator. Previous studies constrained the evolution of q$_{IR}$ with redshift, finding a mild but significant decline, that is yet to be understood. For the first time, we calibrate q$_{IR}$ as a function of textit{both} stellar mass (M$_{star}$) and redshift, starting from an M$_{star}$-selected sample of $>$400,000 star-forming galaxies in the COSMOS field, identified via (NUV-r)/(r-J) colours, at redshifts 0.1$<$z$<$4.5. Within each (M$_{star}$,z) bin, we stack the deepest available infrared/sub-mm and radio images. We fit the stacked IR spectral energy distributions with typical star-forming galaxy and IR-AGN templates, and carefully remove radio AGN candidates via a recursive approach. We find that the IRRC evolves primarily with M$_{star}$, with more massive galaxies displaying systematically lower q$_{IR}$. A secondary, weaker dependence on redshift is also observed. The best-fit analytical expression is the following: q$_{IR}$(M$_{star}$,z)=(2.646$pm$0.024)$times$(1+z)$^{(-0.023pm0.008)}$-(0.148$pm$0.013)$times$($log~M_{star}$/M$_{odot}$-10). The lower IR/radio ratios seen in more massive galaxies are well described by their higher observed SFR surface densities. Our findings highlight that using radio-synchrotron emission as a proxy for SFR requires novel M$_{star}$-dependent recipes, that will enable us to convert detections from future ultra deep radio surveys into accurate SFR measurements down to low-SFR, low-M$_{star}$ galaxies.