  Feynman loop integrals are a key ingredient for the calculation of higher
order radiation effects, and are responsible for reliable and accurate
theoretical prediction. We improve the efficiency of numerical integration in
sector decomposition by implementing a quasi-Monte Carlo method associated with
the CUDA/GPU technique. For demonstration we present the results of several
Feynman integrals up to two loops in both Euclidean and physical kinematic
regions in comparison with those obtained from FIESTA3. It is shown that both
planar and non-planar two-loop master integrals in the physical kinematic
region can be evaluated in less than half a minute with $\mathcal{O}(10^{-3})$
accuracy, which makes the direct numerical approach viable for precise
investigation of higher order effects in multi-loop processes, e.g. the
next-to-leading order QCD effect in Higgs pair production via gluon fusion with
a finite top quark mass.
