© 2020 ACM. For point sets A, B ⊂ Rd, |A| = |B| = n, and for a parameter ϵ > 0, we present a Monte Carlo algorithm that computes, in O(npoly(log n, 1/ϵ)) time, an ϵ-approximate perfect matching of A and B under any Lp-norm with high probability; the previously best-known algorithm takes Ω(n3/2) time. We approximate the Lp-norm using a distance function, d(⋅, ⋅) based on a randomly shifted quad-tree. The algorithm iteratively generates an approximate minimum-cost augmenting path under d(⋅, ⋅) in time proportional, within a polylogarithmic factor, to the length of the path. We show that the total length of the augmenting paths generated by the algorithm is O(n/ϵ)log n), implying that the running time of our algorithm is O(npoly(log n, 1/ϵ)).