Efficient algorithm for computing a compositional square root

Question

Suppose that we have a finite alphabet, for example X = [1,2,3,..., n] and a map f with keys in X and values also in X (a mathematical function from X to X).

Definition: A functional square root, or compositional square root, or half-iterate, is another map g with keys in X and values in X such that g[g[x]]==f[x] for all x in X.

Example 1: If f = {'A':'A', 'B':'A', 'C':'A'} then g = {'A':'A', 'B':'A', 'C':'B'} can be a functional square root of f. The functional square root is not unique. The map f is also its own functional square root.

Example 2: The map f = ['A':'B', 'B':'A'] doesn't have a functional square root. In fact, if we define g['A'] = 'A', then 'A' would be a fixed point of g and therefore it would have to be a fixed point of f. If we define g['A'] = 'B', then we would need to define g['B'] = 'B' and 'B' would have to be a fixed point of f. So, sometimes the functional square root doesn't exist.

I have computed this by hand sometimes, but now I was thinking that I don't really have an efficient strategy to compute it.

Question: [Design an efficient algorithm to] Compute a functional square root given an input map.

I think it is fair to work out first the case that the input is guaranteed to have a functional square root.

Existence of an algorithm: There is always the brute force solution of generating all maps g from X to X and composing each with itself to check if they are equal to f.

---

Some observations I have seen so far that allow to reduce the problem for some inputs:

1. Each fixed point, a key x such that f[x] == x, allows one to define g[x] = x.
2. Each key y that is mapped to a fixed point, f[y] == x and f[x] == x, allows one to define g[y] = x.
3. No other keys, besides the ones satisfying the condition in point (2), should be mapped by g to a fixed point of f.
4. A case that there is an efficient algorithm: If f is injective (no two keys are sent to the same value) then f is a permutation of X in that case a square root can be computed from its cycle decomposition of the permutation. Finding the cycle decomposition can be done in one pass and a square root from that also.

The fact that for bijections the square root can be computed (it looks) linearly makes me hopeful that there may be some efficient way to compute the general case. If I got my definitions right, the problem is NP since a given g one can verify that g is a compositional square root of f in linear time. To me it doesn't smell like NP-complete, but I have very little experience with proving those reductions from one problem to another.

---

The problem can also be seen in the following way. The input map f can be seen as a directed graph of a finite automaton in which each state has out-degree 1. Computing the compositional square root would be to find an automaton with the same set of states that would behave like the original one when stepping two steps at a time.

Answer 1

There is a square matrix creatable for any differentiable function known as Carleman matrix, defined as:

carleman[j,k] = (dk/dx^k (f(x)^j)) / k!

Where n! is factorial of n and x=0. For example, the 3x3 Carleman matrix of e^x is:

[
  [1, 0, 0],
  [1, 1, 0.5],
  [1, 1, 2]
]

Because of the usage of the derivative, you can see the Taylor expansion of e^x on the second row, the Taylor expansion of (e^x)^2 or e^2x on the third row, and so on. This will be important later.

Note that matrix multiplying two Carleman matrices of functions f(x) and g(x), results in the Carleman matrix of f(g(x)), so if you SQUARE the matrix of f(x), you get matrix for f(f(x)). This means that if we can find the "square root" of a function's Carleman matrix, we get it's half iterate's Carleman matrix for free.

Luckily, people have already found numerical methods to solve this problem. One iterative method is called the Denman-Beavers Method, defined as:

Y[k+1] = (Y[k]*(Z[k]^-1))/2
Z[k+1] = (Z[k]*(Y[k]^-1))/2

Where Y[0] is initialized to the input matrix, and Z[0] is initialized to the n-by-n identity matrix. I‘ve implemented Carleman matrix and the Denman-Beavers square root (up to 500 iterations) in numerical.js Github library. Upon taking the square root of the previously mentioned 3x3 Carleman matrix of e^x, I got:

[
  [1, 0, 0],
  [0.5, 0.8944271891673152, 0.22360679791657895],
  [0.23606797789307893, 0.8944271883343149, 1.3416407866664737],
]

(Note that the irrational numbers here seem suspiciously close to sqrt(5)/10 and 4*sqrt(5)/2, probably a result of the Denman iteration.)

The Taylor expansion of the new half-iterate Carleman matrix will show up on the second row, so we can approximate the half-iterate of e^x as

f(x)=0.5+0.8944271883343149x+0.22360679791657895x^{2}

Here is a plot of f(f(x)) (blue) vs e^x (green):

You can see there is a tangent line at x=0.189, although e^x slowly catches up after it outraces the polynomial. This strategy can be applied to basically any function in Javascript with Require.js and my numerical library installed using:

require(["numerical"], function(numerical) {
console.log((numerical.Matrix.carleman(<function>,
<number of expansion terms>)).sqrt().array[1]);
}

But note that my Taylor polynomial is currently numerically unstable and it will explode after 4 or 5 terms (an update is coming soon), so don't try using more than that many terms.

Answer 2

The problem can be broken down into two problems.

If you take the intersection of the ranges of f, f^2, f^3, ... you'll come up with a finite set Y that f operates as a permutation on. As you noted in observation 4, we can find all of the square roots of a permutation. There may be many such, but we know what choices we could have made on the way to finding them. We either pair up cycles of the same length with some rotation, or we reorder a cycle of odd length. If we can't do that, there is no square root, and all of the ways of making a square root can be enumerated by the ways to make those two choices.

The second problem is the transient values, that form chains ending in Y. The chains are easy to find, and we know where they end. If g is to be a square root of f, then a chain of length 2k or 2k-1 for g has to be a pair of chains of f of length k and either k or k-1 such that the end of the first is g(x) where x is the end of the second. That is a very strong condition, because either those chains both end in the same cycle of odd length with the possibly shorter one ending almost halfway back in p, or else those chains both end in different cycles of the same length, and came from a cycle of even length.

I would therefore suggest tackling the problem as follows.

First turn X into a graph where you have a directed arrow from x to f(x). Find all cycles and chains. The union of all cycles is Y, and first verify that f has a square root in Y. After that, sanity check that the chains could be matched up. If not, there is no square root.

If you pass a couple of sanity checks, try to actually match up the chains recursively into a working square root for the part of the graph that contains them. If this succeeds, then you can easily match up the rest of the cycles and you have a square root. If it fails, then there is no square root.

In most cases, you'll find there is no square root fairly fast. The possibly slow approach is the matching of the chains recursively. Hopefully there aren't too many of those decisions, and they are fairly easy to make, but an exponential backtracking is not, a priori, impossible to find.