The peak sidelobe level (PSL) of an orthogonal waveform set directly governs the range and azimuth ambiguities of synthetic aperture radar (SAR) images and is therefore a key figure of merit for the imaging quality of multiple-input multiple-output SAR (MIMO-SAR). Minimizing the PSL is an NP-hard problem with a non-differentiable objective, and existing approaches often suffer from high computational cost and limited scalability while achieving only suboptimal PSL suppression. We therefore propose a loss function built around the 1.5-entmax sparse transform, which is simultaneously differentiable, sparse, and adaptive. Analysis of its differentiability and derivation of its analytical gradient allow us to recast the NP-hard problem as a differentiable optimization problem that can be solved by well-established algorithms. Owing to the sparsity of the transform, the gradient is concentrated on the high-energy sidelobes, while the gradient contribution from low-energy sidelobes becomes exactly zero, removing the gradient noise contributed by low-energy sidelobes and yielding a lower PSL. To overcome the scale sensitivity of the 1.5-entmax function, the upper bound on the sidelobe magnitude is used to normalize the input, which makes the threshold adaptive and removes the need for additional loss-function hyperparameter tuning across waveform-design tasks of different sizes. To further improve computational efficiency, we combine gradient descent with a deep learning framework: the waveform phases are treated as the learnable parameters of a neural network-like model, thereby yielding a back-propagation-based optimization framework with graphics processing unit (GPU) parallelism. The algorithm is implemented for parallel execution on GPU, and the gradient is computed efficiently through the network’s automatic differentiation in conjunction with a custom analytical-gradient operator, leading to a substantial increase in computational speed. Without the need for the manual tuning of loss-function hyperparameters, the proposed algorithm achieves the lowest PSL across waveform sets of various sizes and reduces computation time by approximately two orders of magnitude in large-scale settings.