Paraxial diffraction modeling based on the Fourier transform has seen widespread implementation for simulating the response of a diffraction-limited optical system. For systems where the paraxial assumption is not sufficient, a class of algorithms has been developed that employs hybrid propagation physics to compute the propagation of an elementary beamlet along geometric ray paths. These “beamlet decomposition” algorithms include the well-known Gaussian beamlet decomposition (GBD) algorithm, of which several variants have been created. To increase the computational efficiency of the GBD algorithm, we derive an alternative expression of the technique that utilizes the analytical propagation of beamlets to tilted planes. We then use this accelerated algorithm to conduct a parameter-space search to find the optimal combination of free parameters in GBD to construct the analytical Airy function. The experiment is conducted on a consumer-grade CPU, and a high-performance GPU, where the new algorithm is 34 times faster than the previously published algorithm on CPUs, and 67,513 times faster on GPUs.