Despite many years of development in integration-by-parts reduction, reconstructing all reduction coefficients, which are rational functions of kinematic variables and the space-time dimension, remains a non-trivial problem. The main difficulty comes from the large number of unknowns in a general ansatz, which can lead to a linear system that is too large to solve. In this paper, we present an algorithm for reconstructing reduction coefficients through a sequence of intermediate bases. The resulting analytic reduction coefficients are products of a few analytic matrices, whose non-zero entries are simple rational functions. We demonstrate the efficiency of this algorithm with two cutting edge examples: a three-point massive box-triangle and a four-point massive pentagon-triangle. Reconstructing all reduction coefficients for the box-triangle (pentagon-triangle) requires 3289 (13013) numerical samplings, significantly fewer than the number of unknowns in the general ansatz, 1407406 (21638331).