Feynman integrals obey linear relations governed by intersection numbers, which act as scalar products between vector spaces. We present a general algorithm for the construction of 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 decomposition of a few Feynman integrals at one and two loops, as first steps toward potential applications to generic multiloop integrals. The proposed method can be more generally employed for the derivation of contiguity relations for special functions admitting multifold integral representations.