The synthesis of the required target spectrum with monochromatic light-emitting diodes (LED) is of great significance in reality. When multiple components of the LEDs are needed, and the accuracy required for the target spectrum synthesis is high, solving the problem of the proportion of the components becomes a combination optimization problem with a non-negative solution of the over-determined set of linear equations. Generally speaking, the approximately global solutions can be found for the heuristic-based methods. However, the convergence speed to the worldwide optimizer is low partly because the objective functions analysis properties, such as gradient information, are not used. Gradient-based algorithms converge to local solutions fast and with high accuracy, but the requirement for non-negative solutions in the problem limits their global convergence. Meanwhile, the least squares information of the objective function was not fully utilized in previous research. In this paper, a two-stage optimization algorithm, named LLR_LBFGS, is proposed based on the mathematical analytic properties, the quadratic nonlinear format of the objective function, and the non-negative requirements for the final solution. In the first stage, unconstrained linear fitting is carried out to obtain the unique solution of the least squares theory of the expressions. In the second stage, the non-negative global optimal solution of the problem is further obtained with the help of the constrained quasi-Newtonian method LBFGS. Taking the fitting of standard target spectra CIE-A, CIE-D65, CIE-D50, CIE-D55 and CIE-D75 as the research object, the new method is compared with the Lasso Regression Algorithm (LASSO), Ridge Regression Algorithm (RIDGE), Differential Evolution Algorithm (DE), Particle Swarm Optimization (PSO) and Genetic Algorithm (GA) in solving the same problem in terms of accuracy, running speed, and the decision coefficient R2. The numerical results based on actual industrial cases show that the LLR_LBFGS converge faster, and the solution accuracy is higher because the information on the objective function is utilized more efficiently. Its universality indicates excellent potential to be applied more generally to solve the problem of LED spectral fitting. A more flexible solution to solve the spectral matching problem can be set up according to the design ideas of this paper. This also has significant implications for improving the effectiveness of intelligent optimization methods for finding the optimal solution for LED spectral matching.