We present a simple formula for the overlap integrals of two sets of multi-dimensional harmonic oscillators. The oscillators have in general different equilibrium points, force constants, and natural vibration modes. The formula expresses the overlap matrix in the one-dimensional case, [m'\n "], as a so-called LU decomposition, [m'\n "] = [0'\0 "]Sigma LmtUtn, where the summation index has a range 0 less than or equal to t less than or equal to min(m,n), i.e., it is the matrix product of a lower-triangular matrix L with an upper-triangular Li. These matrices are obtained from simple recursion formulae. This form is essentially retained in the multi-dimensional case. General matrix elements are obtained by exact and finite expressions, relating them to matrix elements over a single set of harmonic oscillator wave functions. We present test calculations with error estimates, also comparing with literature examples. (C) 1998 Elsevier Science B.V.