Given a function $$f(x)$$ and $$\frac{\partial f(x)}{\partial x_i}=\frac{f^2(x1,...,x_i+\pi/2,...,x_n)-f^2(x1,...,x_i-\pi/2,...,x_n)}{f(x)}$$. When $$f(x)\to0$$, $$\frac{\partial f(x)}{\partial x_i}$$ could be infinitely large. ($$f^2(x1,...,x_i+\pi/2,...,x_n)-f^2(x1,...,x_i-\pi/2,...,x_n)$$ is always non-zero)

I have very little experience in deal with this situation in gradient descent process...In my code, $$f(x)$$ is in continuous domain but for purpose to simulate some real world process, $$f(x)$$ is sampled to be discrete and would return values uniformly distributed over $$[0,1]$$. Assume discrete $$f(x)$$ has $$N$$ identity values, at the beginning there is a training set of size $$M$$ ($$M$$ is very large), $$\{x_i, f(x_i)=\frac{k_i}{N}\}_{i=1..M} (k_i \in 1, 2, ..., N)$$.

I found that setting $$1/f(x)$$ to some value like $$0.01$$ when $$f(x)=0$$ would reach the optimizim easily but slightly slower than ideal process, while set to much smaller value like $$0.00001$$ would let $$f(x)=0$$ have a great impact on the process and failed to form a descent curve.

Is the method replacing infinitely large values to some large but finite values correct? Or there are any better ways to deal with the infinite gradient problem?

Yes. For example, the same problem happens for the logarithm in cross-entropy loss function, i.e. $$p_i \text{log}(p'_i)$$ when $$p'_i \rightarrow 0$$. This is avoided by replacing $$\text{log}(x)$$ with $$\hat{\text{log}}(x) = \text{log}(x+\epsilon)$$ for some small $$\epsilon$$.
Similarly, you are changing $$f(x)$$ in the denominator to $$\hat{f}(x) = max(\epsilon, f(x))$$.
However, I would suggest $$\hat{f}(x) = f(x) + \epsilon$$ instead of a cut-off threshold. This way, the difference in $$f(x_1) < f(x_2) < \epsilon$$ would not be ignored unlike the max cut-off.