Four-dimensional (4D) constellations with up to 131 072 points (17 bit/4D-sym) are designed for the first time using geometric shaping. The constellations are optimized in terms of mutual information (MI) and generalized MI (GMI) for the additive white Gaussian noise (AWGN) channel, targeting a forward error correction (FEC) rate of 0.8 at finite signal-to-noise ratios. The presented 15-17 bit constellations are currently the highest-performing constellations in the literature, having a gap to the AWGN capacity as low as 0.17 dB (MI) and 0.45 dB (GMI) at 17 bit/4D-sym. For lower cardinalities, our constellations match or closely approach the performance of previously published optimized constellations. We also show that (GMI-)optimized constellations with a symmetry constraint, optimized for a FEC rate of 0.8, perform nearly identical to their unconstrained counterparts for cardinalities above 8 bit/4D-sym. A symmetry constraint for MI-optimized constellations is shown to have a negative impact in general. The proposed procedure relies on a Monte-Carlo-based approach for evaluating performance and is extendable to other (nonlinear) channels. Stochastic gradient descent is used for the optimization algorithm for which the gradients are computed using automatic differentiation. This article is part of the theme issue 'Celebrating the 15th anniversary of the Royal Society Newton International Fellowship'.