This study tailors an efficient optimization procedure for multiple-leak identification in the frequency domain. Firstly, the maximum likelihood estimation is adopted to build a multi-dimensional optimization problem whose objective function is only in terms of the leak locations. Then, closed-form formulations to compute the gradient of the proposed objective function are presented. Finally, to demonstrate the effectiveness of the developed analytical formulations, the steepest-descent method with multiple initializations is implemented to solve the optimization problem. An initialization algorithm based on a set of points with maximum leak reflections is proposed. The use of the most classical gradient-based technique to find the optimum allowed for a deeper understanding of the multi-dimensional optimization required to pinpoint multiple leaks in transient-based methods in terms of complexity and solvability of the problem. The cases of three leaks and various numbers of resonant frequencies with a clear manifestation of the performance of gradients and their applications in the steepest descent method are investigated and discussed. The procedure is also validated against a laboratory case study having two leaks. The main advantages of the proposed technique, which distinguishes it from previous studies, are (i) the versatility of the process to find leaks with even a single spatial measurement vector and (ii) the efficiency and accuracy of the underlying optimization, as closed-form formulations are employed in the scheme. This work may provide a benchmark for investigating and assessing other methods dealing with multiple leaks due to its comprehensive presentation of the local extrema and the entire identification process.