We present an algorithm computing recurrence relation coefficients for bivariate polynomials, orthonormal with respect to a discrete inner product. These polynomials make it possible to give the solution of a discrete least squares approximation problem. To compute these polynomials, we pose the inverse eigenvalue problem and solve it efficiently and in a stable way, using a sequence of Givens rotations. We also show how to generalize the algorithm for the case of polynomials in more variables. Several numerical experiments show the validity of the approach.