First, let's define \$N\$ to be the number of cells and \$T\$ to be the number of tiles. For example, Pac-Man has about 36 unique tiles (\$T = 36\$) and has about \$28 \times 31\$ cells (\$N = 868\$).
In general, the upper bound is \$O(T^2 \cdot N^2)\$. In practice this probably can be reduced down to \$O(T^2 \cdot N)\$.
We can't really do better than \$O(T \cdot N)\$ because just to traverse all cells and remove all non-chosen tiles, we're already at \$T \cdot N\$. To see why it might be \$O(T^2 \cdot N^2)\$, we might be required to traverse the whole grid after a single tile removal to see if any tiles need to be removed from the given constraints.
In practice, we only have pairwise neighbor tile/cell constraints and so the "impact" of a single tile removal can be localized to only cells around impacted site.
The algorithm for this approach flags neighboring cells as needing to be inspected if a tile was removed near them. After inspection, remove them from consideration, marking any neighbors as needing inspection if any further updates (tile removals) have occured.
Here's some rough pseudo-code:
function propagateConstraints(grid, accessed) {
while (still culling) {
removeFlaggedCells()
foreach cell in accessed {
foreach neighbor_cell in cell {
if tiles can be removed from neighbor_cell {
add neighbor_cell to accessed
flag neighbor_cell
}
}
remove cell from accessed and unflag cell
}
}
}
In general, this only updates a kind of "wave-front" or boundary of where a tile has been removed. In the worst case this can cycle through the whole grid but this only happens if it needs to, whereas a naive version touches many cells whether it needs to or not.
I'm kind of waiving my hands here and only considering the naive way of testing whether tiles need to be updated by doing a complete pairwise comparison. In theory, if one were sufficiently motivated, there might be speedups there to try and reduce the \$O(T^2 \cdot N)\$ down to something closer to \$O(T \cdot N)\$.
I've further ignored the dimension of the grid and taken for granted that the number of neighbors is constant. There is another factor of \$D\$, the dimension that comes into the runtime. For fixed dimension, this obviously doesn't impact the "big-oh" but is still a relevant factor when practically considering run time.
I wrote a small article about my implementation of Wave-Function-Collapse which can be found here. The source code is also available here (though it might be a little hard to follow) and is libre/free, so feel free to inspect or use it as a reference implementation.