We propose the successive inverse polynomial interpolation method to optimize model parameters in subgrid parameterization for large-eddy simulation. This approach is illustrated for the Smagorinsky eddy-viscosity model used in homogeneous decaying turbulence. The optimal Smagorinsky parameter is resolution dependent and provides minimal total error in the resolved kinetic energy. It is approximated by starting with a "bracketing interval" that is obtained from separate "no-model" and "dynamic eddy-viscosity" large-eddy simulations. The total error level is reduced 3-6 times compared to the maximal initial errors. The computational overhead of the full optimization at resolution N³ is comparable to a single simulation at (3N/2)³ grid cells. The increased accuracy is higher than obtained with dynamic modeling at a resolution of (4N)³.