Meta-Minsky Circles and Guillochés

I just ran across this lovely page about guillochés, and something in it bothered me slightly. I reflected on that feeling, and I realized what it was: I felt like the complexity of using transcendental functions was somehow too large for the spare and parsimonious loveliness that is a guilloché. That is, no reflection on Max Irwin, but we ought to be able to generate these designs on a computer without invoking floating-point sums of Tschebysheff polynomials or table lookups, because of course the natural world does not work in that way.

After some thought, this led to the JS implementation of guillochés as linear dynamical systems that you see in the gray canvas here.

However, this system is much more general than guillochés. You can explore its 12-dimensional space by randomizing the parameters repeatedly and mousing over the canvas; if it stops being interesting, you can reset it to return to the original guilloché.

The HAKMEM-149 algorithm

My initial uneasiness led me to think about Minsky’s circle algorithm in HAKMEM item 149, (cf. AIM-239 item 149 on “page 73”, the 74th page of 107) which approximates sines and cosines discretely using a simple linear dynamical system:

NEW X = OLD X - ε * OLD Y
NEW Y = OLD Y + ε * NEW(!) X

Which is to say, p⃗_i =

	1	-ε
	ε	1-ε²

p⃗_i-1. (And, as Minsky points out, if ε is a power of 2, you can evaluate that without even doing any multiplication, which mattered more to efficiency when he wrote that in 1972 than it does today.)

The First HAKMEM-149 Generalization: With Matrices

There are a couple of interesting things about this matrix. First, it’s nearly diagonal: a third of its diagonal elements are 1, and another third are nearly 1, while of its 30 off-diagonal elements, 22 are 0 and four more are nearly 0. (This results in its spatial behavior being fairly coherent, since the state-variable changes from one step to the next are small, and its causality graph being fairly simple.) Second, as a result, it is a fairly sparse matrix. Third, it is a function of only two variables, ε₀ and ε₁, (and nearly a linear one); these two represent the angular velocities of its two oscillators. The space of basic guilloché rosettes up to uniform scaling and translation is two-dimensional, but these two variables do not represent its two dimensions; rather, their ratio represents one of the two dimensions, a common factor between them is a step size parameter controlling the granularity of the approximation (which may result in chaotic behavior!) while the relative radii or amplitudes of the oscillators is determined by their initial state. And it is singular, because its rank is only 4, which is to say that it collapses a six-dimensional input space into four dimensions (because its first two columns are entirely zero). But within those four dimensions (that is, if we amputate the first two rows and columns) it is unitary.

These curious properties imply that we are not likely to stumble across a matrix like this if we simply randomly explore the 36-dimensional space of 6×6 matrices, because almost all matrices are invertible and almost no invertible matrices are unitary. (I hope the reader will forgive me if this is too basic; my grasp of linear algebra is fairly weak. Conversely, I hope the reader will forgive me if this is too complicated, and perhaps take a detour through the basics of linear algebra, which is super fun!) But, if we want to explore similar systems, we could perhaps decompose it into a nearly-diagonal N×N hidden-state-transition submatrix and an arbitrary 2×N projection into pixel space. In this way, the original Minsky circle algorithm becomes capable of drawing arbitrarily scaled and rotated ellipses while its state space remains two-dimensional.

The Second HAKMEM-149 Generalization: With Reversible Instruction Sequences

A slightly different generalization of the circle algorithm which assigns a lower Kolmogorov complexity to the “basic rosette guilloché” variant above is a sequence of instructions of the form xᵢ ← xᵢ ± xⱼ >> k, where i ≠ j, followed by a selection of disjoint subsets of the xᵢ for the pixel coordinates. Among other advantages, this constrains the system to be invertible (as long as (a+b)-b = a is valid, which is true for integer math but not floating point), since this basic instruction preserves xⱼ, and converts a uniform distribution on k into a roughly exponential distribution on ε. It has the disadvantage that To represent this as a series of numeric (i, j, k) tuples, we can arbitrarily choose to represent the ± as the sign of k, so the original Minsky circle system is simply the sequence [(0, 1, -lg ε), (1, 0, lg ε)].

The Third HAKMEM-149 Generalization: With Matrices as Reversible Instruction Sequences

It happens that this particular sequence of instructions has an exhaustive set of (i, j) pairs lexicographically sorted, as does the “basic guilloché” matrix above, and if we restrict our attention to sequences that share this property, we can elide the (i, j) pairs entirely, leaving us with only the ks, which additionally need only 4 to 7 bits each. (If k is 5 bits, its minimum value is -32, and we can take -xⱼ >> 32 to be 0 if our state variables are 32-bit integers; but 16-bit integers may be adequately interesting.)

This reduces the representation of the underlying dynamical system of Minsky’s circles to some 8 to 14 bits, of which some 2⁴ to 2⁷ instances are instances of Minsky’s circles, which is small enough that you could imagine finding the Minsky behavior either by exhaustive enumeration or simply by search. The “basic guilloché” 4×4 matrix, then, becomes 4×3×[4 to 7] = [48 to 84] bits.

We should expect the constraint to invertibility to result in a more interesting state space because it excludes systems with point attractors — the system can have fixed points, but no other states can feed into them. This should become less important as the dimensionality increases, but I expect many of the more interesting systems to appear at low dimensionalities because our minds and especially visual cortices are not equipped to infer large numbers of hidden variables from projections into two dimensions.

If we were to allow the ks to take arbitrary real values instead of simple integers, then we could explore continuous variation over this state space, but because the orbits are generally going to be chaotic, this will not result in continuous variation of the image.

Systems in this family are eminently suited to simple digital hardware implementation, and audio output from them should also be interesting, since it should be periodic and produce a variety of somewhat distorted sinewaves. Continuous-time analog implementations may also be possible, for example, with each state variable represented by a capacitor and each non-diagonal matrix element implemented by the gain on an op-amp, or in current-mode with each state variable represented by an inductor. The correspondence between such analog implementations and the discrete-time algorithm is necessarily imperfect for a variety of reasons, but both may be interesting systems.

Now, therefore, be it resolved that here I will write code to systematically explore and visualize this family of low-Kolmogorov-complexity dynamical systems.

Periodicity and Reversibility

I posted this comment on Mark VandeWettering’s post about the algorithm:

With fixed-point (or generally non-arbitrary-precision) arithmetic, this algorithm is a deterministic finite-state automaton with no input and, say, 16 bits of state in X and 16 bits of state in Y. Any DFA with no input has eventually-periodic behavior, perhaps with a degenerate period-1 single-state loop, but it doesn't necessarily follow that its initial state is part of whatever periodic attractor it ends up in; it could be a sequence like A, B, C, B, (repeat). But for this to be true, some state must have more than one predecessor — in this example, B is preceded by both A and C. This means that the DFA's transition function does not have a deterministic inverse: if you construct a time-reversed version of it, it's nondeterministic. If the transition function is reversible (information-preserving), then you can run the DFA backward to the initial state (and beyond!) from any eventual successor state.

But the state transition function implemented by the HAKMEM-149 algorithm is reversible, as noted in HAKMEM-150, because each of the two steps in the loop, “x -= ε·y” and “y += ε·x”, are individually invertible — in fixed-point arithmetic, anyway. (In floating point, maybe you lose low-precision bits if the change results in floating the point one binary place to the left, so this doesn't apply. And of course in “modern C” signed overflow is undefined, because to compiler vendors, the correctness of your code is less important than fractional percentage point performance increases.)

Therefore, with fixed-point arithmetic with the traditional overflow semantics, every possible initial state is a point on its own attractor — we eventually return to the original point, QED.

(However, as so vividly illustrated in “Minskys and Trinskys”, the periodic behavior from roundoff error and overflow need not look anything like a circle.)

In floating-point, the algorithm is not necessarily reversible, but the above argument demonstrates that it is still necessarily periodic — although conceivably with periods of up to 2⁶⁴ or 2¹²⁸. In practice, it converges to periodicity very rapidly.

This reversibility argument applies equally well to the second and third generalizations of Minsky’s algorithm described above; as long as the instruction that computes the amount by which to change xᵢ does not depend on the previous value of xᵢ, it is guaranteed to be reversible and periodic and therefore eventually return to its starting point. However, you don’t need to increase the dimensionality much before the potential periodicity time becomes large compared to the lifetime of the universe.

The Method of Finite Differences

Babbage designed his Difference Engine to tabulate polynomials using the Method of Finite Differences, a sort of generalization of Pascal’s Triangle sometimes called the Method of Vanishing Triangles or the Method of Difference Tables. This is a very simple way to extrapolate a polynomial sequence. First, for a degree-N sequence, you write down N+1 terms of the sequence. Suppose we want to tabulate y = x² - 2x + 1; first we write down its values for x = 0, 1, 2, and just for fun, 3 and 4:

1 0 1 4 9

Now we calculate the first forward differences Δy of this sequence y, writing below each number its difference from the number before:

1	0	1	4	9
1	-1	1	3	5

These are the first differences Δyᵢ = yᵢ - yᵢ₋₁, which are in a precise sense a discrete analogue of the derivative dy/dx. So now we continue and calculate Δ²y:

1	0	1	4	9
1	-1	1	3	5
1	-2	2	2	2

And Δ³y:

1	0	1	4	9
1	-1	1	3	5
1	-2	2	2	2
1	-3	4	0	0

At this point our finite differences have petered out in an infinite sequence of zeroes; Δ³y will never start having nonzero values if we follow it further, because the underlying function here is a second-order polynomial.

This difference table is also known as a “vanishing triangle” because (if you ignore the gray elements on the left of each row that depend on things before the beginning of the sequence, which are normally not considered part of it) it has a triangular shape, and if the input sequence is polynomial, sooner or later you get to a row where the triangle is all zeroes. Then the leftmost diagonal column of this triangle down to where it vanishes — here, (1, -1, 2) — provides enough information to reconstruct the whole original sequence.