I want to make a clustered heat map for the differentially expressed genes I identified with EdgeR where the batch effect was modeled as an experimental factor. To do this, I used the limma package removeBatchEffect() after TMM-normalization and log2 CPM transformation on the DE subset of the raw expression matrix and it worked.

Is there any concrete practical benefit in doing removeBatchEffect() on the entire dataset first (but with some independent filtering first), and then take out the rows corresponding to the DE genes for visualization? In terms of robustness and accuracy?

In other words, is it better to take out differentially expressed genes from the raw matrix first, and then run the batch correction described above or first run the batch correction procedure on the full raw matrix and only then take out the genes you want to visualize? My DE set is on the order of a few thousand de novo contigs. How much does the size of this DE set matter for the answer?

Doesn't matter. The batch correction is done separately for each gene, so whether it's subsetted or not is irrelevant to the function as it'll give the same result for each gene:

y <- matrix(rnorm(10*6),10,6)
colnames(y) <- c("A1","A2","A3","B1","B2","B3")
y[,1:3] <- y[,1:3] + 10

# All together:
batch <- c("A","A","A","B","B","B")
out <- removeBatchEffect(y, batch=batch)

# One at a time:
out2 <- removeBatchEffect(y[1,,drop=FALSE], batch=batch)
all.equal(out[1,,drop=FALSE], out2) # should be TRUE

If you subset first, it'll be faster, but that shouldn't be a major issue unless you've got lots of features.