During the recent COVID-19 pandemic, the instantaneous reproduction number, R ( t ), has surged as a widely used measure to target public health interventions aiming at curbing the infection rate. In analogy with the basic reproduction number that arises from the linear stability analysis, R ( t ) is typically interpreted as a threshold parameter that separates exponential growth ( R ( t ) > 1) from exponential decay ( R ( t ) < 1). In real epidemics, however, the finite number of susceptibles, the stratification of the population (e.g. by age or vaccination state), and heterogeneous mixing lead to more complex epidemic courses. In the context of the multidimensional renewal equation, we generalize the scalar R ( t ) to a reproduction matrix, R ( t ) , which details the epidemic state of the stratified population, and offers a concise epidemic forecasting scheme. First, the reproduction matrix is computed from the available incidence data (subject to some a priori assumptions), then it is projected into the future by a transfer functional to predict the epidemic course. We demonstrate that this simple scheme allows realistic and accurate epidemic trajectories both in synthetic test cases and with reported incidence data from the COVID-19 pandemic. Accounting for the full heterogeneity and nonlinearity of the infection process, the reproduction matrix improves the prediction of the infection peak. In contrast, the scalar reproduction number overestimates the possibility of sustaining the initial infection rate and leads to an overshoot in the incidence peak. Besides its simplicity, the devised forecasting scheme offers rich flexibility to be generalized to time-dependent mitigation measures, contact rate, infectivity and vaccine protection.