At $e^+ \, e^-$ colliders the QED--initial state radiation forms a large part of the radiative corrections. Their precise and fast evaluation is an essential asset for the experiments at LEP, the ILC and the FCC-ee, operating at high luminosity. A long standing problem in the analytic calculation of the $O(\alpha^2)$ initial state corrections concerns a discrepancy which has been observed between the result of Berends et al. (1988) in the limit $m_e^2 \ll s$ and the result by Bl{\"u}mlein et al. (2011) using massive operator matrix elements deriving this limit directly. In order to resolve this important issue we recalculated this process by integrating directly over the phase space without any approximation. For parts of the corrections we find exact solutions of the cross section in terms of iterated integrals over square root valued letters representing incomplete elliptic integrals and iterations over them. The expansion in the limit $m_e^2 \ll s$ reveals errors in the constant $O(\alpha^2)$ term of the former calculation and yields agreement with the calculation based on massive operator matrix elements, which has impact on the experimental analysis programs. This finding also explicitly proofs the factorization of massive initial state particles in the high energy limit including the terms of $O(\alpha^2)$ for this process.