Feynman integrals obey linear relations governed by intersection numbers, which act as scalar products between vector spaces. We present a general algorithm for constructing multivariate intersection numbers relevant to Feynman integrals, and show for the first time how they can be used to solve the problem of integral reduction to a basis of master integrals by projections, and to directly derive functional equations fulfilled by the latter. We apply it to the derivation of contiguity relations for special functions admitting multi-fold integral representations, and to the decomposition of a few Feynman integrals at one- and two-loops, as first steps towards potential applications to generic multi-loop integrals.