In order to reconstruct the 3D surface from gradient fields quickly and accurately, a new fast and accurate least squares integration algorithm was proposed. Compact finite difference scheme was introduced into optimization equation for better accuracy. Then the objective function was represented as a Sylvester function. With Hessenberg-Schur algorithm, the space and time complexity were reduced from O(N2) and O(N3) to O(N) and O(N3/2), respectively. The experiment result showed that when the 4th-order compact scheme is used, the accuracy of the new method is improved by one order higher than Higher-order Finite-difference-based Least-squares Integration(HFLI) and Global Least-Squares(GLS). While with 6th-order compact scheme, the accuracy is improved by one order higher than Spline-based Least-squares Integration(SLI). The robustness of the proposed method is weaker than that of HFLI and SLI, but better than GLS. The reconstruction speed was obviously faster than that of HFLI and SLI.