  We propose a nonperturbative framework to study general correlation functions
of single-trace operators in $\mathcal{N}=4$ supersymmetric Yang-Mills theory
at large $N$. The basic strategy is to decompose them into fundamental building
blocks called the hexagon form factors, which were introduced earlier to study
structure constants using integrability. The decomposition is akin to a
triangulation of a Riemann surface, and we thus call it hexagonalization. We
propose a set of rules to glue the hexagons together based on symmetry, which
naturally incorporate the dependence on the conformal and the R-symmetry cross
ratios. Our method is conceptually different from the conventional operator
product expansion and automatically takes into account multi-trace operators
exchanged in OPE channels. To illustrate the idea in simple set-ups, we compute
four-point functions of BPS operators of arbitrary lengths and correlation
functions of one Konishi operator and three short BPS operators, all at one
loop. In all cases, the results are in perfect agreement with the perturbative
data. We also suggest that our method can be a useful tool to study conformal
integrals, and show it explicitly for the case of ladder integrals.
