We present a numerical method for polynomial approximation of multivariate functions. The method utilizes Gauss quadrature in tensor product form, which is known to be inefficient in high dimensions. Here we demonstrate that by using a new randomized algorithm and taking advantage of the tensor structure of the grids, a highly efficient algorithm can be constructed. The new method does not require prior knowledge/storage of the entire data set at all the tensor grid points, whose total number of points is excessively large in high dimensions. Instead the method utilizes one data point at a time and iteratively conducts the approximation. This feature allows the use of the method irrespective of the size of the data set. We establish the rate of convergence of this iterative algorithm and show that its operational counts can be lower than the standard methods, when applicable, such as least squares. Numerical examples in up to hundreds of dimensions are presented to verify the theoretical analysis and demonstrate the effectiveness of the method.