Josephson junctions in a two-dimensional electron gas with spin-orbit coupling are a promising candidate to realize topological superconductivity. While it is known that the geometry of the junction strongly influences the size of the topological gap, the question of how to construct optimal geometries remains unexplored. We introduce a greedy numerical algorithm to optimize the shape of Majorana junctions. The core of the algorithm relies on perturbation theory and is embarrassingly parallel, which allows it to explore the design space efficiently. By introducing stochastic variations in the junction Hamiltonian, we avoid overfitting geometries to specific system parameters. Furthermore, we constrain the optimizer to produce smooth geometries by applying image filtering and fabrication resolution constraints. We run the algorithm in various setups and find that it reliably produces geometries with increased topological gaps over large parameter ranges. The results are robust to variations in the optimization starting point and the presence of disorder, which suggests the optimizer is capable of finding global maxima.