Let $g_1 \cdots g_M$ be the basic gates that you are allowed to use. For the purposes of this $\operatorname{CNOT}_{12}$ and $\operatorname{CNOT}_{13}$ etc are treated as separate. So $M$ is polynomially dependent on $n$, the number of qubits. The precise dependence involves details of the sorts of gates you use and how $k$-local they are. For example, if there are $x$ single qubit gates and $y$ 2-qubit gates that don't depend on order like $CZ$ then $M = xn+\binom{n}{2}y$.
A circuit is then a product of those generators in some order. But there are multiple circuits that do nothing. Like $\operatorname{CNOT}_{12} \operatorname{CNOT}_{12} = \mathrm{Id}$. So those give relations on the group. That is it is a group presentation $\langle g_1 \cdots g_M \mid R_1 \cdots \rangle$ where there are many relations that we do not know.
The problem we wish to solve is given a word in this group, what is the shortest word that represents the same element. For general group presentations, this is hopeless. The sort of group presentation where this problem is accessible are called automatic.
But we can consider a simpler problem. If we throw out some of the $g_i$, then the words from before become of the form $w_1 g_{i_1} w_2 g_{i_2} \cdots w_k$ where each of the $w_i$ are words only in the remaining letters. If we managed to make them shorter using the relations that don't involve the $g_i$, then we will have made the entire circuit shorter. This is akin to the optimizing the CNOTs on their own made in the other answer.
For example, if there are three generators and the word is $aababbacbbaba$, but we don't want to deal with $c$, we will instead shorten $w_1=aababba$ and $w_2=bbaba$ to $\hat{w}_1$ and $\hat{w}_2$. We then put them back together as $\hat{w}_1 c \hat{w}_2$ and that is a shortening of the original word.
So WLOG (without loss of generality), let's suppose we are in that problem already $\langle g_1 \cdots g_M \mid R_1 \cdots \rangle$ where we now use all the gates specified. Again this is probably not an automatic group. But what if we throw out some of the relations. Then we will have another group that has a quotient map down to the one we really want.
The group $\langle g_1 g_2 \mid - \rangle$ no relations is a free group, but then if you put $g_1^2=\mathrm{id}$ as a relation, you get the free product $\mathbb{Z}_2 \star \mathbb{Z}$ and there is a quotient map from the former to the later reducing the number of $g_1$'s in each segment modulo $2$.
The relations we throw out will be such that the one upstairs (the source of the quotient map) will be automatic by design. If we only use the relations that remain and shorten the word, then it will still be a shorter word for the quotient group. It just won't be optimal for the quotient group (the target of the quotient map), but it will have the length $\leq$ to the length it started with.
That was the general idea, how can we turn this into a specific algorithm?
How do we choose the $g_i$ and relations to throw out in order to get an automatic group? This is where knowledge of the kinds of elementary gates we typically use comes in. There are a lot of involutions, so keep only those. Keep careful attention to the fact that these are only the elementary involutions, so if your hardware has a hard time swapping qubits that are vastly separated on your chip, this is writing them in only the ones that you can do easily and reducing that word to be as short as possible.
For example, suppose you have the IBM configuration. Then $s_{01},s_{02},s_{12},s_{23},s_{24},s_{34}$ are the allowed gates. If you wish to do a general permutation, decompose it into $s_{i,i+1}$ factors. That is a word in the group $\langle s_{01},s_{02},s_{12},s_{23},s_{24},s_{34} \mid R_1 \cdots \rangle$ that we wish to shorten.
Note that these don't have to be the standard involutions. You can throw in $R(\theta) X R(\theta)^{-1}$ in addition to $X$ for example. Think of the Gottesman-Knill theorem, but in an abstract manner that means it will be easier to generalize. Such as using the property that under short exact sequences, if you have finite complete rewriting systems for the two sides, then you get one for the middle group. That comment is unnecessary for the rest of the answer, but shows how you can build up bigger more general examples from the ones in this answer.
The relations that are kept are only those of the form $(g_i g_j)^{m_{ij}} = 1$. This gives a Coxeter group and it is automatic. In fact, we don't even have to start from scratch to code up the algorithm for this automatic structure. It is already implemented in Sage (Python based) in general purpose. All you have to do is specify the $m_{ij}$ and it has the remaining implementation already done. You might do some speedups on top of that.
$m_{ij}$ is really easy to compute because of the locality properties of the gates. If the gates are at most $k$-local, then the computation of $m_{ij}$ can be done on a $2^{2k-1}$ dimensional Hilbert space. This is because if the indices don't overlap, then you know that $m_{ij}=2$. $m_{ij}=2$ is for when $g_i$ and $g_j$ commute. You also only have to compute less than half of the entries. This is because the matrix $m_{ij}$ is symmetric, has $1$'s on the diagonal ($(g_i g_i)^1 = 1$). Also most of the entries are just renaming the involved qubits so if you know the order of $(\operatorname{CNOT}_{12} H_1)$, you know the order of $\operatorname{CNOT}_{37} H_3$ without doing the computation over again.
That took care of all relations that only involved at most two distinct gates (proof: exercise). The relations that involved $3$ or more were all thrown out. We now put them back in. Let's say we have that, then one can perform Dehn's greedy algorithm using new relations. If there was a change, we knock it back up to run through the Coxeter group again. This repeats until there are no changes.
Every time the word is either getting shorter or staying the same length and we are only using algorithms that have linear or quadratic behaviour. This is a rather cheap procedure so might as well do it and make sure you didn't do anything stupid.
If you want to test it out yourself, give the number of generators as $N$, the length $K$ of the random word you're trying out and the Coxeter matrix as $m$.
edge_list=[]
for i1 in range(N):
for j1 in range(i):
edge_list.append((j1+1,i1+1,m[i1,j1]))
G3 = Graph(edge_list)
W3 = CoxeterGroup(G3)
s3 = W3.simple_reflections()
word=[choice(list([1,..,N])) for k in range(K)]
print(word)
wTesting=s3[word[0]]
for o in word[1:]:
wTesting=wTesting*s3[o]
word=wTesting.coset_representative([]).reduced_word()
print(word)
An example with N=28 and K=20, the first two lines are the input unreduced word, the next two is the reduced word. I hope I didn't make a typo when entering the $m$ matrix.
[26, 10, 13, 16, 15, 16, 20, 22, 21, 25, 11, 22, 25, 13, 8, 20, 19, 19, 14, 28]
['CNOT_23', 'Y_1', 'Y_4', 'Z_2', 'Z_1', 'Z_2', 'H_1', 'H_3', 'H_2', 'CNOT_12', 'Y_2', 'H_3', 'CNOT_12', 'Y_4', 'X_4', 'H_1', 'Z_5', 'Z_5', 'Y_5', 'CNOT_45']
[14, 8, 28, 26, 21, 10, 15, 20, 25, 11, 25, 20]
['Y_5', 'X_4', 'CNOT_45', 'CNOT_23', 'H_2', 'Y_1', 'Z_1', 'H_1', 'CNOT_12', 'Y_2', 'CNOT_12', 'H_1']
Putting back those generators like $T_i$ we only put back the relations like $T_i^n = 1$ and that $T_i$ commutes with gates that do not involve qubit $i$. This allows us to make the decomposition $w_1 g_{i_1} w_2 g_{i_2} \cdots w_k$ from before have the $w_i$ as long as possible. We want to avoid situations like $X_1 T_2 X_1 T_2 X_1 T_2 X_1$. (In Cliff+T one often seeks to minimize T-count). For this part, the directed acyclic graph showing the dependency is crucial. This is a problem of finding a good topological sort of the DAG. That is done by changing precedence when one has a choice of what vertex to use next. (I wouldn't waste time optimizing this part too hard.)
If the word is already close to optimal length, there is not much to do and this procedure won't help. But as the most basic example of what it finds is if you have multiple units and you forgot there was an $H_i$ at the end of one and an $H_i$ at the beginning of the next, it will get rid of that pair. This means you can black box common routines with greater confidence that when you put them together, those obvious cancellations will all be taken care of. It does others that aren't as obvious; those use when $m_{ij} \neq 1,2$.
