  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.
