IDA* Optimizations

We employ a number of tricks to improve the running time of the Cycle Combination Solver’s IDA* tree search.

We enhance the speed of puzzle operations through the use of puzzle-specific SIMD on AVX2 and Neon instruction set architectures. Namely, the VPSHUFBVPSHUFB instruction on AVX2 and the tbl.8tbl.8/tbl.16tbl.16 instructions on Neon perform permutation composition in one clock cycle, enabling for specialized SIMD algorithms to compose two Rubik’s Cube states and test for a Cycle Combination Solver solution. They have both been disassembled and highly optimized at the instruction level. Additionally, the puzzle-specific SIMD uses compacted representations optimized for the permutation composition instructions. For example, it uses a representation of a Rubik’s Cube state that can fit in a single YMMYMM CPU register on AVX2 and in the DD and QQ CPU registers on Neon.

Pruning table generation also uses puzzle-specific SIMD. To generate a pruning table on the corners orbit, we need to use a different Rubik’s Cube representation because we don’t want to waste CPU caring about what happens to edges. So, every orbit has its own specialized SIMD representation and SIMD algorithm modifications.

We leave the precise details at the prescribed references; we defer our discussion of how the SIMD algorithms work for a later revision.

At every increasing depth level of the IDA* search tree we explore times as many nodes. We formally call this number the branching factor—the average number of child nodes visited by a parent node. A few clever observations can reduce the branching factor.

We observe that we never want to rotate the same face twice. For example, if we perform followed by , we’ve just reversed the move done at the previous level of the tree. Similarly if we perform followed by another , we could have simply done straight away. In general, any move should not be followed by another move in the same move class, the set of all move powers. This reduces the branching factor of the child nodes from for all moves to . Additionally, we don’t want to search both and because they commute, and result in the same net action. So, we assume that (or ) never follows (or ), and in general, we only permit searching distinct commutative move classes strictly in a single order only. Move sequences that satisfy these two conditions are called canonical sequences. Canonical sequences are special because these two conditions make it easy to check if a move sequence in the search tree is redundant.

What does the second condition reduce our branching factor from to? We start by counting the number of canonical sequences at length , denoted , using a recurrence relation. We consider the last move of the sequence , the second to last move , and the third to last move . The recurrence relation can be constructed by analyzing two cases:

• Case 1: and do not commute.

  In this case, is simply multiplied by the number of possibilities of . Since and do not commute, cannot be () nor its opposite face (). Therefore, must be one of move classes, or one of the possible moves. We can establish that the first component in the recurrence relation for is .

• Case 2: and commute.

  We need to be careful to only count and , one time so we count them in pairs. In this case, is simply multiplied by the number of strictly ordered pairs. There are pairs of commutative move classes: . We have to discard one of these pairs because necessarily commutes with the move classes in one of these pairs since the union of all of these pairs is every move. Such a canonical sequence where the subsequence all commute cannot exist because one of those moves will always violate the strict move class ordering. For example, if is and is , then there is no possible option for that makes the full sequence a canonical sequence.

  Each move class in each pair can perform three moves, which implies that each pair contributes possible moves. Overall we find this number to be possible moves. We can establish that the second component in the recurrence relation for is .

can be thought of as the superposition of these two cases with the base cases (exercise to the reader: figure out where these come from). Hence, . The standard recurrence relation can be solved as follows:

The term dominates as approaches infinity; our new branching factor is approximately !

It turns out that is not an exact bound on the number of distinct positions at sequence length but merely an upper bound. This is because the formula overcounts, and the actual number is always lower: it considers canonical sequences that produce equivalent states such as and as two distinct positions. It turns out it is extremely nontrivial to describe and account for these equivalences, to the point where it’s not worth doing so: at shallow and medium depths, roughly stays within of the actual distinct position count. The Cycle Combination Solver considers the extra work negligible and searches equivalent canonical sequences anyways. The Big O time complexity of IDA* can be realized as , an improvement over from.

The Cycle Combination Solver uses an optimized finite state machine to perform the canonical sequence optimization.

We use a special form of symmetry reduction during the search we call sequence symmetry, first observed by Rokicki and improved by our implementation. Some solution to the Cycle Combination Solver conjugated by yields , which we observe to be a rotation of the original sequence as well as a solution to the Cycle Combination Solver by the properties of conjugation discussed earlier. Repeatedly applying this conjugation:

forms an equivalence class based on all the rotations of sequences that are all solutions to the Cycle Combination Solver. The key is to search a single representative sequence in this equivalence class to avoid duplicate work.

Similarly to symmetry conjugation, we choose the representative as the lexicographically minimal sequence on a move-by-move basis (with a move class ordering relation defined). Unlike symmetry conjugation, we don’t manually apply all sequence rotations to find the representative; rather, we embed sequence symmetry as a modification to the recursive IDA* algorithm such that it only ever searches the representative sequence. We do this by observing that if a representative sequence starts with move , then every other move cannot be lexicographically lesser than it. If this observation were to be false, we could keep on rotating the sequence until the offending move is at the beginning of the sequence, and since that move is lexicographically lesser than that sequence rotation would be the true representative. This contradicts the initial representative sequence assumption. We permit moves that are lexicographically equal to (i.e. in the same move class) but change the next recursive step to repeat the logic on the move after . The overall effect is that the IDA* algorithm only visits move sequences such that no later subsequence is lexicographically lesser than the beginning of the move sequence. This suffices for the complete sequence symmetry optimization.

The modification described is not yet foolproof. The sequence would technically be valid as there is no later subsequence lesser than the beginning, but the actual lexicographically minimal representative is the sequence rotation. The “later subsequence” of the true representative wraps around from the end to the beginning. So, extra care must be taken at the last depth to manually account for the wrapping behavior. We only apply this to the last depth, so sequences like are still searched by the next depth limit of IDA*.

We can extend our prior definition of canonical sequences to include sequence symmetry as a third condition. How does sequence symmetry affect the number of canonical sequences at depth ? Because a sequence of length has sequence rotations, sequence symmetry logically divides the total number of nodes visited by , but only in the best case. The canonical sequence only has members in its sequence rotational equivalence class, not , so the average value to divide by is actually a bit less than . It follows that the average number of canonical sequences at depth (and the IDA* asymptotic time complexity) is bound by and . Testing has shown this number to typically be right in the middle of these two bounds.

Furthermore, we take advantage of the fact that the optimal solution sequence almost never starts and ends with commutative moves. We claim that the IDA* algorithm almost never needs to test such that and commute for a solution. The proof is as follows.

We first observe that if is a solution, then is also a solution by a sequence rotation. This tells us that and cannot be in the same move class or else they could be combined to produce the shorter solution . Such a shorter solution would have been found at the previous depth limit, implying that never would have been explored, making this situation an impossibility. This also tells us that also cannot be in a greater move class than because would be a lexicographically lesser than , contradicting our earlier proof that IDA* only searches the lexicographically minimal sequence rotation (the representative). Therefore, must be in a lesser move class than .

If is a solution, then is also a solution because and commute. By the transitive property, if is a solution, then so is . Both of these sequences are independently searched and tested as a solution because there is no direct “commutative move ordering” or sequence symmetry relation between them. This is redundant work; we choose to discard the case. This completes our proof.

This optimization only applies to the last depth in IDA*, so it only prevents running the test to check if a node is a solution and does not affect the time complexity. It turns out to be surprisingly effective at reducing the average time per node because most of the time is spent at the last depth.

We alluded to an edge case when we said “almost never.” If doesn’t exist, or if every move from commutes with and , then this optimization will skip canonical sequences where every move commutes with each other; for example on the Rubik’s Cube. The number of skipped sequences is so small that we have the bandwidth to manually search and test these sequences for solutions before running IDA*.

We use a simple optimization described by Mérõ called pathmax to prune nodes with large child pruning heuristics. When a child node has a large pruning heuristic, we can set the current node cost to that value minus one and re-prune to avoid expanding the remaining child nodes. This larger heuristic is still admissible because it is one less than a known lower bound, and the current node is one move away from all of its child nodes. This is only effective when the heuristics are inconsistent, or, in this case, when the pruning table entries are the minimum of two or more other values. With exact pruning tables only, this optimization will never run because the entries are perfect heuristics that cannot exhibit this type of discrepency.

Our last trick is to enhance IDA* through the use of parallel multithreaded IDA* (PMIDA*). PMIDA* runs in two phases. In the first phase, we use BFS to explore the state space to a shallow depth, maintaining a queue of all of states at the last search depth. In the second phase, we use a thread pool to run IDA* in parallel for every state in that queue, utilizing of all of the CPU cores on the host machine. To uphold the optimality guarantee, PMIDA* synchronizes the threads using a barrier that triggers when they have all completed exploring the current level. It can be thought of as a simple extension to the familiar IDA* algorithm.

There have been many parallel IDA* algorithms discussed in literature; how do we know PMIDA* is the best one? We take advantage of the special fact that the Cycle Combination Solver starts searching from the solved state. In order to understand this, we compare the total Rubik’s Cube position counts with the Rubik’s Cube position counts that are unique by symmetry.

Recall that our theoretical branching factor is . In the table of Rubik’s Cube position counts, the branching factor roughly matches this number. However, at the shallow depths of the table of Rubik’s Cube position counts unique by symmetry antisymmetry, our branching factor is much less because there are duplicate positions when performing moves from the solved state. Intuitively, this should make sense: the Rubik’s Cube is not scrambled enough to start producing unique positions. It is easy to pick out two sequences of length two that are not unique by symmetry; for example and . The branching factor converges to its theoretical value as the Rubik’s Cube becomes more scrambled because symmetric positions become more rare. In fact, it was shown by Qu that scrambling the Rubik’s Cube can literally be modelled as a Markov chain (it’s almost indistinguishable from a random walk of a graph). Hence, it is unlikely for two random move sequences of the same length to produce positions equivalent by symmetry. We know that such collisions do happen because the branching factor doesn’t actually reach the value, but we consider them negligible.

The effectiveness of the PMIDA* algorithm stems from combining all of these observations. When our initial shallow BFS search is done, we filter out the many symmetrically equivalent positions from the queue to avoid redundant work before we start parallelizing IDA*. The savings are incredibly dramatic: at depth , for example, we symmetry reduce the number of nodes from to . This is a reduction by , a factor that is close to the familiar (the number of symmetries antisymmetries). Once we do that, and the cube starts to become sufficiently scrambled, we are confident to claim that each IDA* thread worker explores their own independent regions of the search space and duplicates a negligible amount of work.

We make note that there are almost always going to be more positions in the queue to parallelize than available OS threads. We use an optimized thread pool work stealing algorithm for our multithreaded implementation.

We squeeze out our last bit of juice by overlapping pruning table memory latency with the computation. It has been empirically observed that random access into the pruning table memory is the dominating factor for Rubik’s Cube solvers. Modern processors include prefetching instructions that tell the memory system to speculatively load a particular memory location into cache without stalling the execution pipeline to do so. Our PMIDA* implementation uses a technique described by Rokicki called microthreading to spend CPU time on different subsearches while waiting for the memory to come to a query. It splits up each thread into eight “slivers” of control. Each sliver calculates a pruning table query memory address, does a prefetch, and moves on to the next sliver. When that sliver gets control again, only then does it reference the actual memory. By handling many subsearches simultaneously, microthreading minimizes the CPU idle time.

How does PMIDA* affect the asymptotic time complexity? We established in an upper bound of . The time required by PMIDA* can be computed by adding the time of the first and second phases. In the first phase the time required for the BFS is where is the aforementioned shallow depth. In the second phase we symmetry reduce at the shallow depth, split the work across independent threads, and ignore nodes before depth . The time required is where is the number of symmetries antisymmetries. The PMIDA* time complexity is thus , but we consider to be very small and to be a negligible constant. As such the final time complexity becomes . We can apply the exact same logic to our lower bound, and we get .