We modify Feynman integrals by introducing an auxiliary parameter \eta into each Feynman propagator. The modified Feynman integrals are analytical functions of \eta, and physical Feynman integrals can be obtained by taking \eta -> 0^+. Asymptotic expansions of the modified Feynman integrals at \eta -> \infinity can be very easily calculated, which contain only equal-mass vacuum integrals. Due to uniqueness theorem of analytical functions, the asymptotic series determine both values of physical Feynman integrals and their linear relations. Based on the asymptotic expansion, we construct an algorithm to generate linear relations that can reduce arbitrary multi-loop Feynman integral to master integrals. Computation complexity of numerically solving these linear relations is O(N) when the number of linear equations N is large. Our method may overcome the difficulty of IBP method, where linear equations are coupled and thus computation complexity is O(N^3).