We complete the computation of the two-loop helicity amplitudes for the production of three photons at hadron colliders, including all contributions beyond the leading-color approximation. We reconstruct the analytic form of the amplitudes from numerical finite-field samples obtained with the numerical unitarity method. This method requires as input surface terms for all relevant five-point non-planar integral topologies, which we obtain by solving the associated syzygy problem in embedding space. The numerical samples are used to constrain compact spinor-helicity ansätze, which are optimized by taking advantage of the known one-loop analytic structure. We make our analytic results available in a public C++ library, which is suitable for immediate phenomenological applications. We estimate that the inclusion of the subleading-color contributions will decrease the size of the two-loop corrections by about 30% to 50%, and the NNLO cross sections by a few percent, compared to the results in the leading-color approximation.