Algorithms for Skein Manipulation in a Genus-2 Handlebody

Rachel Kinard
Air Force Research Laboratory, Sensors Directorate (AFRL/RYAT), 2241 Avionics Cir, WP AFB, OH 45433

Răzvan Gelca
Texas Tech University

Paul T. Schrader
Air Force Research Laboratory, Information Directorate (AFRL/RIGC), 525 Brooks Road, Rome, NY 13441

[2022]Primary  18D10, Secondary  16T05, 17A99

We present a series of algorithms for skein manipulation in a genus-2 handlebody, implementing a novel strand sorting method to reduce any skein to a skein in a 2-punctured disk. This reduction guarantees resolution as a linear combination of basis elements of the Kauffman Bracket Skein Module. Manually, these skein manipulations prove to be computationally intensive due to the inherent exponential nature of skein relations (i.e., a skein diagram with $n$ crossings yields $2^n$ new skein diagrams, each in $\mathbb{C}[t,t^{-1}]$, the Laurent polynomials with complex coefficients). Thus, as the number of crossings in a skein diagram increases, manual computations become intractable and automation desirable. We enable the automation of all skein computations in the genus-2 handlebody by first converting the skein diagram into an equivalent array, reducing the task of performing skein computations to that of implementing array operators, and then proving that we can always recover the resulting complex Laurent polynomial.

§ INTRODUCTION

§.§ Motivation

Traditional skein manipulation (the art of producing algebraic representations from a diagrammatic array) is a delicate combination of creativity, care, and complexity employing very carefully executed manual calculations. Knots with only a few crossings (e.g., the trefoil) still require significant computational effort. Even an expert in these computations must be vigilant; the introduction of so much as a minor error or a faulty assumption risks propagation throughout the entire computation, producing incorrect resulting polynomial representations. Skein computations are becoming increasingly important to scientific application in knot theoretical studies, quantum field theory, and quantum information science, arising from the study of knot and link invariants. In these applications, the exponential nature of skein computations and the scale at which these computations must be performed makes automation desirable, producing consistent, reliable results and refocusing the mathematician's effort on more productive theoretical pursuits. Within skein theory exist surprisingly simple patterns and remarkable relationships awaiting discovery; it often happens that a skein computation involving numerous terms simplifies to a relation with very few terms.

§.§ History

Historically, skein theory began with John Conway's discovery of a skein relation for the Alexander polynomial, and was formalized by Vaughn Jones' knot polynomial in 1984 [1]. Unlike the Alexander polynomial, the Jones polynomial could only be defined by skein relations. Following Jones, Louis Kauffman defined a relation and link invariant, the Kauffman Bracket, also computed via skein relations [2]. Edward Witten [3] gave an intrinsic definition of the Jones polynomial using quantum field theory based on the Chern-Simons functional (1989), and this theory was rigorously developed by Nikolai Reshetikhin and Vladmir Turaev [4], [5]. (We note that there is a parallel theory constructed by Christian Blanchet, Nathan Habegger, Gregor Masbaum, and Pierre Vogel which makes use of the Kauffman bracket in [6]; see also [7]). In particluar, we replicate the skein computations in the genus 2 handlebody of J. Przytycki [8]. Skein computations in the genus-2 handlebody are well understood, however, the ideal way to demonstrate the computational dilemmas arising in even these simple skein manipulations is by providing some rudimentary examples.

§.§ Our Setting

To begin, let $H$ be a genus-2 handlebody. A knot in $H$ is an embedding of a circle $\mathbb{S}^{1} \hookrightarrow H$. A link is an embedding of several disjoint copies of $\mathbb{S}^{1}$, ${\bigsqcup}_{i\in I}\mathbb{S}^{1} \hookrightarrow H$, where the index set $I$ numbers the components of the link. We refer to these embeddings as “skeins”. Similar to skeins of yarn, they are equivalent up to ambient isotopy (strands can be manipulated isotopically by combinations of Reidemeister moves without affecting the result). Figure <ref> shows these embeddings in the genus-2 handlebody.

Recall the three Reidemeister moves shown in Figure <ref>. Application of Reidemeister moves can affect the number of crossings we see in a given projection of the skein system. However, the projection itself also affects the seen number of crossings and any computation or series of computations we perform has to be consistent regardless of the projection we choose.

A framed knot in $H$ is the embedding of an annulus, $S^{1} \times I \hookrightarrow H$, where $I$ is the unit interval. It is customary to represent a framed knot by the same diagram as a unframed knot with the convention that the framing is parallel to the plane of the diagram (the "blackboard" framing). We manipulate framed knots in the same manner (by Reidemeister moves) except that we must take the framing into account whenever we twist/untwist any portion of the skein (R1).

Skeins embedded in a genus-2 handlebody, H

Skein computations are traditionally performed diagrammatically. A skein relation, shown in Figure <ref>, resolves an apparent crossing in a skein as a diagrammatic combination of two other skeins. Resolving all crossings in a skein produces a Laurent polynomial with skein terms. This requires $2^{n}$ computations, where $n$ is the number of crossings in the projection of the original skein. Thus, the resulting Laurent polynomial with coefficients in $\mathbb{C}$, or $\mathbb{C}[t,t^{-1}]$, has $2^{n}$ skein terms. For some of the resulting skeins, we obtain curves which are basis elements of the Kauffman Bracket Skein Module of the genus-2 handlebody (Przytycki); these fundamental curves are shown in Figure <ref>: $x,y,z$. However, in many cases, resolving all crossings in a skein produces a curve which has no self-intersections, but is not a fundamental curve, such as the skein terms at the botton of Figure <ref>. For such curves, how do we proceed? In section <ref>, we propose a method for resolving all skeins in the genus-2 handlebody.

Reidemeister Moves: $R1$, $R2$, and $R3$. Skeins are considered up to isotopy.

Skein Relations (framed) and Skein Computations; resolving a crossing produces two skein terms, each with a coefficient in $\mathbb{C}[t,t^{-1}]$.

Fundamental Curves

As a quick sidenote, although many of the skeins resulting from a skein computation are not fundamental curves, some can be easily realized as combinations of fundamental curves through simple manipulation, application of isotopy, and Reidemeister moves. Curves like those shown in Figure <ref> occur so often we give them a variable name as well. This is strictly unnecessary, but very convenient!

Finally, we note the use of power notation on variables to represent parallel copies of fundamental curves. For example, $x^{2}$ indicates 2 parallel copies of $x$, while $z^{3}$ indicates 3 parallel copies of $z$, etc. Thus, our task of resolving a skein in the genus-2 handlebody amounts to reducing skeins to linear combinations of powers of $x, y, z, y', y''$ and the unknot.

Auxilary curves: $y'=-t^{-2}xz-t^{-4}y$ and $y''=-t^{4}y-t^{2}xz$. These curves occur frequently, and are simple combinations of fundamental curves.

Next, we provide a motivational example of how such skein computations are performed for framed skeins in the genus-2 handlebody.

§.§ Example: A Sample Skein Computation

The skein shown in Figure <ref> can be resolved as a polynomial combination of fundamental curves in $\mathbb{C}[t,t^{-1}]$. In the first row of Figure <ref> the crossings are resolved one at a time using the skein relations from Figure <ref>. Resolving these 2 crossings produces $2^2=4$ terms. We note that the first three terms are collections of the fundamental curves given in Figure <ref> and Figure <ref>. In the second row of Figure <ref>, we consider the last term; this skein can be resolved by performing a Reidemeister $R1$ move to induce a crossing by twisting the lower portion of the skein. We present the result at the bottom of Figure <ref>; we have fully resolved the skein as a linear combination of fundamental curves with coefficients in $\mathbb{C}[t,t^{-1}]$. We can, of course, further simplify the resulting polynomial by distributing and combining like terms; we have produced the skein polynomial associated to the given skein.

In the first row, we resolve both crossings. The second row induces a crossing by a $R1$ move. The third row expresses the skein as its resulting $\mathbb{C}[t,t^{-1}]$ representation.

We note here that inducing the crossing in the lower part of the skein was a choice. We could have instead induced the crossing in the upper portion of the skein. In this case, the choice of induction location does not affect our result, but in many cases it does, sometimes negatively.

§.§ Example: A Problematic Skein Computation

In this case, the choice of induction location affects the result in a negative manner. Consider the given skein in the upper left of Figure <ref>. As in Example <ref> we resolve all existing crossings (1 crossing produces $2^1=2$ terms) and determine how to manipulate the skein to produce fundamental curves. The first term is resolved as in Example <ref> (by inducing a crossing on either the upper portion or lower portion by a Reidemeister (R1) move). However, In the second row of Figure <ref>, we consider the second term; there are multiple places we may induce a crossing. Observe that when we induce a crossing and resolve in the second row of Figure <ref> we recover the original skein. This was a bad location to induce a crossing! Not only did we end up back at the curve we started with, but we have also introduced another term to the computation.

The skein computation described in Example <ref>. A bad choice of crossing location produces the original skein - bad choices compound exponentially!

§.§ A Method for Skein Computations

A more complicated skein: Where do we induce crossings?

Thus far, our method for performing skein computations is to first resolve all the given crossings, cleverly induce new crossings via a series of Reidemeister moves, and then resolve the new crossings we induce. In Example <ref> and <ref>, we see that there is some choice as to where to induce new crossings using Reidemeister moves; some choices are good choices, and result in the resolution of skeins as fundamental (basis) curves. Bad choices result in the recovery of the skein we started with (or an equivalently unknown curve) and the computation produces additional powers and terms which complicate the system even more. Furthermore, given each crossing resolves the skein into 2 terms, bad choices compound exponentially!

More complicated skeins, like the one shown in Figure <ref>, introduce even more potential locations to induce a new crossing (after all existing crossings have been resolved). In practice, and in more complicated settings, only expertise and years of experience dictate how a skein should be cleverly isotoped to be realized as a combination of fundamental curves. The process can be frustrating at best for a novice!

In this paper, we propose a series of algorithms to resolve any skein in a genus-2 handlebody as a combination of fundamental curves given in Figure <ref>, removing the need for clever realization and crossing induction by implementing a sorting method across each strand in the skein projection diagram. This effectively reduces a skein in the genus-2 handlebody $H$ to a skein in a 2-punctured disk $D$. Given that all curves in the $D$ are the desired fundamental curves, this procedure gives us the resolution of any skein in $H$ as basis elements of the Kauffman Bracket Skein Module of $H$. In addition to providing an upper bound on the number of crossings induced (while minimal, we do not guarantee minimum) these algorithms provide a method for automating skein computations (i.e., skeins resolve themselves!) This generalized method removes the need for cleverly realizing more complicated skeins as combinations of known skeins. Additionally, performing all computations with thoughtful reductions at each stage manages the overabundance of terms that increase in computational complexity very efficiently.

In this case, the main difficulty we face in automation is proper input; the current input of our computations are diagrammatic. In order to convert from a diagram to a more algorithmic input (e.g. an array), we must record all the necessary information about our skein system in the given projection. We will establish a notation similar to the format of a braid representation describing the self-intersections of the skein and the relative position of each component to the other components. We then show this representation recovers a unique skein (Lemma <ref>) and apply Algorithms <ref>-<ref> to fully resolve it as a combination of fundamental curves. The result is a complex Laurent polynomial in $x$, $y$, $z$, and the unknot ($t^{-2}+t^{2}$). Note that, although unnecessary, the auxiliary curves $y'$ and $y''$ can be replaced (rather than resolved) by their values in Figure <ref> for a more efficient computation.

§ SKEIN REPRESENTATION: DIAGRAM TO ARRAY AND SKEIN INPUT PARAMETERS

§.§ Diagram to Array

We start by converting our skein system to a braid diagram as shown in Figure <ref>. Note that nothing about the skein itself has changed, and we view the holes in the genus-2 handlebody as fixed strands of an underlying base braid. When we perform skein computations, the strands of the base braid stay fixed; crossings between the skein and these strands need not be resolved - only crossings between the skein and itself (self-intersections) are resolved using the given skein relations.

Braid Representation

Next, we define a reference line at the bottom of the braid diagram (Figure <ref>) and index the crossings between the skein and the base braid from bottom to top along each strand of the braid. We will also index crossings between the skein and itself. The indexing sets between strands 1, 2, and the set of skein self-intersections do not have to be related, and, using a different indexing scheme along each strand will not affect the final result; it is only important that elements along each strand are ordered from bottom to top. Skeins with multiple components, Figure <ref>, are indexed along each strand from bottom to top regardless of to which component that portion of the skein belongs. We will identify each of the crossings in the skein by labeling them on the diagram. The order of the crossings determines the order in which they will be resolved. Note that this ordering is inherently different than the index ordering we assigned. To accentuate this difference, we will label the crossings with rational numbers $(0.01, 0.02, ..., 0.0n)$, not to be confused with the counting numbers we use for indices based on our reference line.

Indexing a single component

Indexing multiple components

In general, we can label our crossings by $1 \times 10^{-k_{r}}, 2 \times 10^{-k_{r}}, ..., r \times 10^{-k_{r}}$ where $k_{r}$ is the magnitude of the number of crossings. For example, a skein with 1-99 crossings are labeled as 0.01, 0.02, ..., 0.99, a skein with 100-999 crossings are labeled 0.001, 0.002,...,0.999, etc. Finally, at the locations where the skein intersects the strands of the base braid we will note whether the skein goes over ($+n$) or under ($-n$) the base braid $n$. With this information and the labeling on the set of skein self-intersections we can form an array, $E^{l}$, for each component $l$ of the skein. Among all components of the skein if there are $s$ self-intersections and $r$ intersections between the skein and the braid then there are $r+2s$ points to record. i.e., the label on each self-intersection is recorded twice, once for the overcrossing and once for the undercrossing, although these over/undercrossings may be in the same or different components of the skein.

We outline how to record skein information at every point along the skein. The points we consider are locations where the skein crosses the strands of the base braid in the given projection. In each component $l$ of the skein $K$, and at each point $p$, we note the skein intersection information ($e^l_p$), skein index information ($i^l_p$), and skein orientation information ($q^l_p$). Recording this information over every component, we prove (Lemma <ref>) that this representation $(E,I,Q)$ uniquely recovers the skein $K$

§.§ Crossing Input Information (E)

In each component $l$ of the skein $K$, let $p_{l}$ be the number of points of interest (places where the skein intersects the strands of the base braid or skein self-intersections). If there are $k$ total components in the skein, for each component $0 < l \leq k$, we form an array:

\begin{equation*}
E^{l}=[e_{0}^{l},e_{1}^{l}, ..., e_{p_{l}-1}^{l}]
\end{equation*}

where $e_{j}^{l}$, $0 \leq j \leq (p_{l}-1)$ is determined by the following, letting $S$ be the set of skein self-intersections:

Skein Intersection Input Information     Value:

If the skein crosses over strand 1 of the base braid:     $e_{j}^{l}=1$

If the skein crosses under strand 1 of the base braid:     $e_{j}^{l}=-1$

If the skein crosses over strand 2 of the base braid:     $e_{j}^{l}=2$

If the skein crosses under strand 2 of the base braid:     $e_{j}^{l}=-2$

If the skein crosses over a self-intersection ($s \in S$)     $e_{j}^{l}=s$

If the skein crosses under a self-intersection ($s \in S$)     $e_{j}^{l}=-s$

We record $E^{l}$ for each component $l$, $0 < l \leq k$, and from this we form the skein intersection array:

\begin{equation*}
\end{equation*}

where $k$ is the number of components in the skein $K$.

We note that we must choose a reference point and direction to list elements in each component of the skein. That is, we choose an element to begin our array arbitrarily. This choice does not affect our final answer. Figure <ref> shows how to label skein intersection information every time the skein crosses over or under the base braid (that is, passes in front or behind the hole in the genus-2 surface from our given projection).

The value of entries in the skein intersection array

§.§ Index Input Information (I)

Next, we form the index array for each component $l$, starting with the point of interest corresponding to $e_{0}^{l}$, (the first entry in $E^{l}$) after our assigned reference point. We proceed in the same order as the elements of $E^{l}$ so that the entries of $I^{l}$ correspond to the entries of $E^{l}$. We record $i^{l}_{j}$ point by point as assigned relative to our line of reference. For each component $0 < l \leq k$:

\begin{equation*}
I^{l}=[i_{0}^{l},i_{1}^{l}, ..., i_{p_{l}-1}^{l}]
\end{equation*}

\begin{equation*}
\end{equation*}

where $k$ is the number of components in the skein $K$.

§.§ Orientation Input Information (Q)

We form an array for the orientation of each component:

\begin{equation*}
Q^{l}=[q_{0}^{l},q_{1}^{l}, ..., q_{p_{l}-1}^{l}]
\end{equation*}

While it is only necessary to define whether the component is traversed clockwise or counterclockwise from the assigned reference point (i.e., the assignment of a single value for each component), we have to make a local decision about orientation when comparing any two elements. Thus, we proceed to assign each element $e_{j}^{l},i_{j}^{l}$ a value $q_{j}^{l}$, $0 \leq j \leq (p_{l}-1)$, corresponding to the location of that point in the given projection, and compare these values between two points to determine a motion from left to right or right to left.

Skein Orientation Input Information     Value:

If $e_{j}^{l}=1$, and the skein goes from left to right:     $q_{j}^{l}=3$

If $e_{j}^{l}=1$, and the skein goes from right to left:     $q_{j}^{l}=4$

If $e_{j}^{l}=2$, and the skein goes from left to right:     $q_{j}^{l}=4$

If $e_{j}^{l}=2$, and the skein goes from left to right:     $q_{j}^{l}=5$

If $e_{j}^{l}=\pm0.0k$:     $q_{j}^{l}=0$

We record $Q^{l}$ for each component $l$, $0 < l \leq k$, and from this we form the skein intersection array:

\begin{equation*}
\end{equation*}

where $k$ is the number of components in the skein $K$.

§.§ Self-intersections (U)

Finally, when performing skein computations for oriented skeins, we record the sign on each crossing. To determine the sign of the crossing, we look at the orientation of the overcrossing and then undercrossing. We keep track of the signs on the crossings to be resolved because the resolution depends on whether the crossing is positive or negative (as defined by the right hand rule for skeins). We record the signs of the crossings at each stage in an array $U$. When a crossing is resolved, the corresponding entry of $U$ is deleted.

When we resolve a crossing, we smooth it in both directions (as shown in Figure <ref>); this smoothing reverses the orientation of one of the sections of the skein. The sign of any other self-intersection (yet to be resolved) that lies along that section may be affected. If both components (the overcrossing and undercrossing) of another crossing lie on this segment, the orientation of both changes, and no sign change is necessary. However, if only one of the components (either the overcrossing or the undercrossing) of another crossing lies on this segment, the sign of that crossing will change in the next stage. Thus, for each new skein term, we record the signs of the remaining crossings for that term with the other skein information.

Thus, we store the signs of crossings to be resolved in an array, called "U", as follows:

\begin{equation*}
\end{equation*}

where $r$ is the number of self-intersections in the original skein. As we resolve these crossings, we delete their corresponding entries from $U$; when the skein $K$ is fully resolved, $U=\emptyset$.

We pause the narrative for two examples demonstrating how to record this information for a skein with one component <ref> and a skein with two components <ref>. We then summerize the skein input array in <ref>.

§.§ Example (Single Component): Crossing, Index, and Orientation Inputs

The skein in Figure <ref> has only one component ($l=1$). We set a reference point, and record the skein intersection information $E^{1}=[1,-0.1,2,-2,-0.2,2,-2,0.2,0.1,-1]$ and $E=[E^1]$. From the reference point shown, the index information is recorded $I^1=[1,8,6,5,7,3,4,7,8,2]$, and $I=[I^1]$. For orientation, we record $Q^1=[3,0,4,5,0,4,5,0,0,4]$, and $Q=[Q^1]$. (Note that $dim(E^l)=dim(I^l)=dim(Q^l)$, since $e_{j}^{l}$ corresponds one-to-one with $i_{j}^{l}$ and $q_{j}^{l}$, and $dim(E)=dim(I)=dim(Q)$.) There are 2 crossings to resolve, labeled 0.1 and 0.2. The array $U$ stores the signs of these crossings: $U=[-1,1]$.

Single Component Crossing, Index, and Orientation Inputs

§.§ Example (Multi-Component): Crossing, Index, and Orientation Inputs

Multi-Component Crossing, Index, and Orientation Inputs

The skein in Figure <ref> has two components ($l=2$). We set a reference point in the first component, assign an orientation direction (counterclockwise), and record the skein intersection information $E^{1}=[-1,1,-1,-0.2,0.3,-1,-1,1]$. We do the same for the second component, $E^{2}=[1,0.1,-2,2,-2,2,0.2,-0.3,-0.1,-1]$. Then, $E=[E^1,E^2]$. The index is given by $I^1=[1,2,3,9,15,4,7,8]$ for component 1, $I^2=[5,14,13,12,11,10,9,15,14,6]$ for component 2, and $I=[I^1,I^2]$. For orientation, $Q^{1}=[3,4,3,0,0,4,3,4]$, $Q^{2}=[3,0,4,5,4,5,0,0,0,4]$ and $Q=[Q^1]$.

§.§ Skein Input Parameters

In describing the initial skein, we produce three arrays: the crossing array $E_0$, the index array $I_0$, and the orientation array $Q_0$ as described above (NOTE: We are using the subscript 0 to indicate that this information is associated to the original skein, before the resolution of any crossings). We can include all this information in an initial array of the form:

where $s_0=0$ is the initial power of the variable $t$ (i.e., there is no variable $t$ prior to the application of the skein relations to resolve a crossing) and $c_0=1$ (the sign and value of the initial coefficient). The array $U_0$ is an array that stores the crossings to be resolved; as crossings are resolved elements are removed from $U$ until $U=\emptyset$. We thus recover the skein as $c_0t^{s_0}K_0 = 1t^0K_0$ where $K_0=[E_0,I_0,Q_0]$ is the complete description of the unique skein recovered by its crossing, index, and orientation arrays (as shown, up to a choice of orientation, in Proposition REF). At each stage in the skein computation, we resolve one of the $r$ crossings, producing $2^{r}$ terms. The final result, after the resolution of all crossings in the initial diagram, is a polynomial with $2^r$ diagrammatic terms, which can be described by array of the form:

\begin{equation*}
\end{equation*}

where each component $Z_z$ represents one terms of the final Laurent polynomial, $Z_z= [s_z,c_z,U_z,E_z,I_z,Q_z]$. We recover a linear combination of skeins in the variable $t$ from $Z$ as:

\begin{equation*}
\end{equation*}

where $K_z=[E_z,I_z,Q_Z]$, and the additional information ($c_z,s_Z$) stores the value of the coefficient in $\mathbb{C}[t,t^{-1}]$ as $c_zt^{s_z}$. Table 1 summarizes each of the entries in the skein array. The fact that the skein can be uniquely recovered from this the information given here is formalized in Proposition <ref>.

Summary of the skein array

2|c|$Z_z= [s_z,c_z,U_z,E_z,I_z,Q_z]$

Skein Array Information:     Description:

$s_z$     power of the variable $t$

$c_z$     sign and value of the constant, $s_zc_z \in \mathbb{C}[t,t^{-1}]$

$U_z=[\pm1,\pm1,...,\pm1]$     signs on crossings to be resolved ($dim(U_z)=r$)

$K_z=[E_z,I_z,Q_z]$     skein input information

$E_z=[E_{z}^{0},E_{z}^{1},...E_{z}^{l},...,E_{z}^{k}]$     skein crossing array

$I_z=[I_{z}^{0},I_{z}^{1},...I_{z}^{l},...,I_{z}^{k}]$     skein index array

$Q_z=[Q_{z}^{1},Q_{z}^{2},...Q_{z}^{l}...,Q_{z}^{k}]$     skein orientation array

The resolution of crossings $1 \times 10^{-k_{r}}, 2 \times 10^{-k_{r}}, ..., w \times 10^{-k_{r}}$ by manipulations of the array $Z=[Z_0,Z_1,...,Z_z,...Z_{2^r-1}]$ produces the linear combination of skeins in $t$ with coefficients in $\mathbb{C}[t,t^{-1}]$:

\begin{equation}
\sum_{z=0}^{2^w-1}c_zt^{s_z}K_z
\end{equation}

where $K_z=[E_z,I_z,Q_z]$ and $w$ is the integer corresponding to the label on the last crossing to be resolved. When $w=r$,

\begin{equation}
\sum_{z=0}^{2^r-1}c_zt^{s_z}K_z
\end{equation}

where $U=\emptyset$ and all crossings in the original crossings in the skein have been resolved.

Finally, we note that any skein embedded in a genus-2 handlebody can be uniquely recovered from the input information required in this paper ($E$, $I$, and $Q$). This is not the minimum representation of such a skein, but it is one which allows for automated computation and efficient manipulation.

The representation $K=[E,I,Q]$ recovers a unique skein in a genus-2 handlebody. Any skein in a genus-2 handlebody has multiple representations, $K$, equivalent up to permutation of the array elements.

We recover the skein $K$ from the given arrays $E$ (skein intersection array), $I$ (index array), and $Q$ (orientation array). In each component $l$, we recover $K^{l}=[E^l,I^l,Q^l]$ as follows: the dimension of the arrays $dim(E^l)=dim(I^l)=dim(Q^l)=h$ is the number of points of interest we consider. Beginning with the array $I^l$ we label points along each strand of the base braid by their sequential index values, from bottom to top, by partitioning $I^l$ into three sets: $I^{l}_{1}$, $I^{l}_{2}$, and $I^{l}_{3}$ according to the corresponding values in $E^{l}$. $I^{l}_{1}$ contains all elements $i^{l}_{j}$ for which $|e^{l}_{j}|=1$. $I^{l}_{2}$ contains all elements $i^{l}_{j}$ for which $|e^{l}_{j}|=2$. $I^{l}_{3}=I^{l} \setminus (I^{l}_{1} \cup I^{l}_{2})$; thus $I^{l}_{3}$ contains all indices corresponding to skein self-intersections. This induces a partition on $E^{l}$ (into $E^{l}_{1}$, $E^{l}_{2}$, and $E^{l}_{3}$). Let $n_{1}=dim(I^{l}_{1})$, $n_{2}=dim(I^{l}_{2})$, and $n_{3}=dim(I^{l}_{3})$. The values of the elements of $I^{l}_{1}$ determine an ordering along strand 1. Similarly for $I^{l}_{2}$. Begin with the minimal element of $I^{l}_{1}$ (the element with the lowest index) and consider the corresponding value of $E^{l}_{l}$; WLOG, call this minimal element $i^{l}_0 \in I^{l}_1 \subset I^{l}$. It's corresponding value is $e^{l}_0 \in E^{l}_1 \subset E^{l}$. (Note: we assume that the array $E^{l}$ begins at the over/undercrossing of strand 1 with the lowest index; starting anywhere else, we can rearrange our array such that this element is the first entry by permuting the elements. The starting point we choose to begin recording skein information is arbitrary; we can always choose to record skein information beginning with the lowest index along strand 1, thus, we assume $|e^{l}_0|=1$). If $e^{l}_0=1$, we draw an overcrossing and if $e^{l}_0=-1$ we draw an undercrossing at this location. Returning to the original array $I^l$, beginning with $i^{l}_{0}$ we connect $i^{l}_{0}$ to $i^{l}_{1}$ by inspecting the value of $e^{l}_{1}$. Continuing in this fashion for all  $(e^{l}_{j},i^{l}_{j})$, we recover the unoriented skein. Now, fix any value $j$ and consider $(e^{l}_{j},i^{l}_{j},q^{l}_{j})$ and $(e^{l}_{j+1},i^{l}_{j+1},q^{l}_{j+1})$. The comparison of $q^{l}_{j}$ and $q^{l}_{j+1}$ determines the orientation of component $l$. Thus, the skein in uniquely determined for each component from the representation $K^{l}=[E^{l},I^{l},Q^{l}]$.

We note that each skein $K$ may have multiple equivalent representations. Indexing the skein based on height we only require the index to be sequential along each strand. The index labeling on the base braid does not affect the labeling on any other strand (independent between strands). Thus, the array $I^l$ used in each component $l$ is not unique.

§ ALGORITHMS FOR SKEIN COMPUTATIONS

We first introduce an algorithm to resolve all crossings $1 \times 10^{-k_r}, 2 \times 10^{-k_r},..., r \times 10^{-k_r}$ in a given skein system. This is split between two sub-algorithms; the first (Algorithm 1) addresses the structure of the algorithm on array $Z$, resolving crossings one by one. The second (Algorithm 2) pulls out the details of that resolution on individual components of $Z$.

$Z_i$ decomposes as a sum of $Z'_i$ and $Z"_i$, coefficients in $\mathbb{C}[t,t^{-1}]$

In Algorithms <ref> and <ref>, beginning with $Z=[Z_0]$ (the initial skein diagram, $1t^0K_0$ and continuing for $Z_i=c_it^{s_i}K_i$) the burden of our computation at each stage is to apply the skein relations to $Z_i$ given in Figure <ref> to produce two new skeins: $Z'_i$ and $Z''_i$. We then replace $c_it^{s_i}K_i$ by the combination $(c_i+c'_{i})t^{s_i+s'_{i}}K'_i+(c_i+c''_{i})t^{s_i+s''_{i}}K''_i$. The way we edit the array $Z_i$ to produce $Z'_i$ and $Z''_i$ (shown in Figure <ref>) depends on several factors: if the crossing to be resolved is a self-intersection in the same component $Z_i$ of the skein (Case 1) only that component needs to be manipulated. If, however, the crossing is between two components of the skein, $Z_i$ and $Z_j$, (Case 2) both components are manipulated and replaced by $Z'_i$ and $Z''_i$.

We perform these array manipulations (representing skein crossing resolution) via a series of operators: $\mathcal{C}$, $\mathcal{R}$, $\mathcal{F}$, and $\mathcal{B}$. In Algorithms 1 and 2 we resolve crossing $r$ in $E_i$ and then apply the same manipulations to $I_i$ and $Q_i$.  For ease of reference, we call the crossing to be resolved “0.r” (assuming $k_r=1$). In Definition <ref>, we define these 4 operators, the input of which are arrays or pairs of arrays with given indices $j'$ and $j''$.

Given arrays $x$ ($dim(x)=m$), $y$ ($dim(y)=n$), $z$ ($dim(z)=k$), and indices $j'$, $j''$, we define the operators $\mathcal{C}$, $\mathcal{R}$, $\mathcal{F}$, $\mathcal{B}$ as follows:

\begin{equation*}
\mathcal{C}: x \xrightarrow[\text{split section, concatenate}]{\mathcal{C}} y, z
\end{equation*}

$\mathcal{C}(x,j',j'') = y,z$

$\text{where}\; y = [x_{j'+1},...,x_{j''-1}]$
$\text{and}\; z = [x_0,...,x_{j'-1},x_{j''+1}...,x_{m-1}]$,

\begin{equation*}
\mathcal{R}: x \xrightarrow[\text{reverse section}]{\mathcal{R}} z \\
\end{equation*}

$\mathcal{R}(x,j',j'') = z$

$\text{where}\; z =[x_0,...,x_{j'-1},x_{j''-1},...,x_{j'+1},x_{j''+1}...,x_{m-1}]$,

\begin{equation*}
\mathcal{F}: x,y \xrightarrow[\text{forward section insert}]{\mathcal{F}} z \\
\end{equation*}

$\mathcal{F}(x,y,j',j'') = z$

$\text{where}\; z =[x_0,...,x_{j'-1},y_{j''+1},..., y_{n-1},y_0,...,y_{j''-1},x_{j'+1},...,x_{m-1}]$,

\begin{equation*}
\mathcal{B}: x,y \xrightarrow[\text{backward section insert}]{\mathcal{B}} z \\
\end{equation*}

$\mathcal{B}(x,y,j',j'') = z$

$\text{where}\; z =[x_0,...,x_{j'-1},y_{j''-1},...,y_0, y_{n-1},...,y_{j''+1},x_{j'+1},...,x_{m-1}]$.

We adjust the array $U_i$ to record signs affected by the orientation changes in the recombination of skeins due to the resolution of a crossing. These changes affect future resolutions; the signs of crossings are only affected when the components in the array are reversed (i.e., when $\mathcal{R}$ and $\mathcal{B}$ are applied). As described in <ref>, we check for crossings in $E_i$ in the section of the skein that is reversed. In Case 1 ( i.e., crossings to be resolved are in the same component) we check for crossings between indices $j'$ and $j"$ (i.e., , any crossings in the array $[x_{j'+1},...,x_{j''-1}]$). If both the overcrossing and undercrossing of the same label lie in this section of the array, no change is needed. If only one part (either the overcrossing or the undercrossing) is in this section of the array we change the corresponding value in $U$ (specifically, the element of $U$ corresponding to sign of that crossing) to the opposite sign. In Case 2, the same condition holds except that we must check for crossings in the array $[y_0,...y_{n-1}]$.

Resolve Crossings in Skein

[1]
Input: $Z$; $dim(Z)=1$
Resolve crossing $r$ in $Z_i$
$Z_i\xrightarrow{\text{Algorithm 2}} Z'_i, Z''_i$
delete $Z_i$
insert $Z'_{i}$ at index $i$
insert $Z''_{i}$ at index $i+1$
increase $i$ by 2
Output: $Z$; $dim(Z)=2^r$

Resolve Crossings in Skein

[1]
Input: $Z_i=[s_i,c_i,U_i,E_i,I_i,Q_i], r$
Locate overcrossing $0.r$ and undercrossing $-0.r$ in $E_i$:
implement any array search for value
$l' \leftarrow$ component of $\pm0.r$ in $E_i$
$j' \leftarrow$ index of $\pm0.r$ in $E_i^{l'}$
$l'' \leftarrow$ component of $\mp0.r$ in $E_i$
$j'' \leftarrow$ index of $\mp0.r$ in $E_i^{l''}$
Case 1: over/undercrossing in same component
$EA_i: E_i^l \xrightarrow[\text{split section, concatenate}]{\mathcal{C}} \mathcal{C}_1(E_i^l),\mathcal{C}_2(E_i^l)$
delete $E_i^l$ in $E_i$
insert $\mathcal{C}_1(E_i^l)$ at index $l$ of $E_i$
insert $\mathcal{C}_2(E_i^l)$ at index $l+1$ of $E_i$
form corresponding $IA_i$, $QA_i$
Repeat Lines 3-6 for $IA_i$ and $QA_i$
adjust $U_i$ to form $UA_i$
Use Definition DEF
define $Z'_i=[s_i+1,c_i,UA_i,EA_i,IA_i,QA_i]$
$EB_i: E_i^l \xrightarrow[\text{reverse section}]{\mathcal{R}} \mathcal{R}(E_i^l)$
delete $E_i^l$ in $E_i$
insert $\mathcal{R}(E_i^l)$ at index $l$ of $E_i$
form corresponding $IB_i$, $QB_i$
Repeat Lines 10-12 for $IB_i$ and $QB_i$
adjust $U_i$ to form $UB_i$
Use Definition DEF
define $Z''_i=[s_i-1,c_i,UB_i,EB_i,IB_i,QB_i]$
Case 2: over/undercrossing in different components
WLOG, assume $l'<l''$
$EA_i: E_i^{l'},E_i^{l''} \xrightarrow[\text{forward section insert}]{\mathcal{F}} \mathcal{F}(E_i^{l'},E_i^{l''})$
delete $E_i^{l''}$ in $E_i$
delete $E_i^{l'}$ in $E_i$
insert $\mathcal{F}(E_i^{l'}),E_i^{l''}$ at index $l'$ of $E_i$
form corresponding $IA_i$, $QA_i$
Repeat Lines 18-21 for $IA_i$ and $QA_i$
adjust $U_i$ to form $UA_i$
Use Definition DEF
define $Z'_i=[s_i+1,c_i,UA_i,EA_i,IA_i,QA_i]$
$EB_i: E_i^{l'},E_i^{l''} \xrightarrow[\text{backward section insert}]{\mathcal{B}} \mathcal{B}(E_i^{l'},E_i^{l''})$
delete $E_i^{l''}$ in $E_i$
delete $E_i^{l'}$ in $E_i$
insert $\mathcal{B}(E_i^{l'}),E_i^{l''}$ at index $l'$ of $E_i$
form corresponding $IB_i$, $QB_i$
Repeat Lines 25-28 for $IB_i$ and $QB_i$
adjust $U_i$ to form $UB_i$
Use Definition DEF
define $Z''_i=[s_i+1,c_i,UB_i,EB_i,IB_i,QB_i]$
Output: $Z'_i$,$Z''_i$

Induce Crossings in Skein

[1]
Input: $Z$; $dim(Z)=2^r$, where $r$ is the number of crossings in the original skein
Induce crossings in $Z_k$
$Z_k$ has no skein self-intersections, $U_k=\emptyset$ (output of Algorithm 1,2)
Form $P_n$, $N_n$, $B_n$
Use Definition <ref>
$P_n \neq \emptyset$
$P_n, N_n, B_n \xrightarrow[\text{induction decision}]{\text{Algorithm 4}} (l',j'), (l'',j''), P'_n, N'_n, B'_n$

Use Algorithm 4 to locate where to induce crossing(s)

$P_n$, $N_n$, $B_n$ are modified in Algorithm 4
Update $P_n \leftarrow P'_n$, $N_n \leftarrow N'_n$, $B_n \leftarrow B'_n$
$Z_k, (l',j'), (l'',j'') \xrightarrow[\text{induce crossings}]{\text{Algorithm 5}} \hat{Z}_k$

Use Algorithm 5 to induce crossings at location $(l',j')$ and $(l'',j'')$

$Z_k$ now has either one (0.1) or two (0.1, 0.2) crossings
Assign $r \leftarrow 1$ Identify induced crossing $0.1$
$\hat{Z}_k, r \xrightarrow[\text{resolve crossings}]{\text{Algorithm 2}} Z'_k, Z''_k$
Use Algorithm 2 to resolve crossing 0.1
delete $Z_k$
insert $Z'_k$ at index $k$
insert $Z''_k$ at index $k+1$

Lines 10-13: the kth term is replaced by the two new skein terms (resolved)

This increases $dim(Z)$ by 1; continue with the new $Z_k=Z'_k$
If 2 crossings were induced in Algorithm 5
$r \leftarrow x$
Identify induced crossing 0.2 in $Z'_K$ of line 12 (now at index k)
$Z_k, r \xrightarrow[\text{resolve crossings}]{\text{Algorithm 2}} Z'_k, Z''_k$
Use Algorithm 2 to resolve crossing $0.2$
delete $Z_k$
insert $Z'_k$ at index $k$
insert $Z''_k$ at index $k+1$
This increases dim(Z) by 1

Identify induced crossing 0.2 in term $Z''_k$ of line 13 (now at index k+2)
$Z_{k+2}, r \xrightarrow[\text{resolve crossings}]{\text{Algorithm 2}} Z'_{k+2}, Z''_{k+2}$
Use Algorithm 2 to resolve crossing $0.2$
delete $Z_{k+2}$
insert $Z'_{k+2}$ at index $k+2$
insert $Z''_{k+2}$ at index $k+3$

Lines 14-23: Each of the terms from Lines 10-13 is replaced by two new terms (resolved)

Thus, the 4 new terms are located at indices k, k+1, k+2, and k+3
Continue with modified $Z_k$ (and associated $P_n, N_n, B_n$) until $P_n=\emptyset$
$Z_k$ is fully resolved; increase $k$ by 1 and consider next component ($Z_{k+1}$)
Output: $Z$; $dim(Z)>>2^r$, each $Z_k$ fully resolved as basis elements $x,y,z$ and the unknot

Induction Decision

[1]
Input: $n$, $B_n$, $P_n$, $N_n$
define $b=min(B_n)$, $c=min(P_n)$

Compare $min(P_n)$ to $min(B_n)$
delete $c$ from $P_n$ $\rightarrow P'_n$
$P'_n=P_n \setminus \{c\}$
$c \leftarrow min(P_n)$ recompute minimum of modified $P_n$
delete $b$ from $B_n$ $\rightarrow B'_n$
$B'_n=B_n \setminus \{b\}$
$b \leftarrow min(B_n)$
recompute minimum of modified $B_n$
Since, $b \neq c$, $min(B_n) \neq min(P_n) \rightarrow min(B_n) = min(N_n)$
Define $M_n=N_n$
$M_n$ is a copy of $N_n$; $M_n$ is modified and discarded
$a \geq c$
Compare $min(P_n)$ to $max(M_n)$
delete $a$ from $M_n$
$a \leftarrow max(M_n)$ (recompute maximum of $M_n$)
now, $a < c$, $max(M_n)<min(P_n)$

Induction Decision: switch indices $a=i_{j'}^{l'}$ and $c=i_{j''}^{l''}$

Output: $j', l', j'', l''$ where $a=i_{j'}^{l'}$ and $c=i_{j''}^{l''}$, $P'_n, N_n, B'_n$

Induce Crossings

[1]
Input: $Z_k$, $l', j', l'', j''$
$l'=l''=l \And |j''-j'|=1$
[same component (consecutive)]
insert $0.1$, $0$, $0$ at index $j''+1$ (after $j''$) of $E_k^l$, $I_k^l$, $Q_k^l$, respectively
insert $-0.1$, $0$, $0$ at index $j'$ (before $j'$) of $E_k^l$, $I_k^l$, $Q_k^l$, respectively

insert $-0.1$, $0$, $0$ at index $j'+1$ (after $j'$) of $E_k^l$, $I_k^l$, $Q_k^l$, respectively
insert $0.1$, $0$, $0$ at index $j''$ (before $j''$) of $E_k^l$, $I_k^l$, $Q_k^l$, respectively
[different components or same component (non-consecutive)]
insert $\pm0.1$, $0$, $0$ at index $j'+1$ of $E_k^{l'}$, $I_k^{l'}$, $Q_k^{l'}$, respectively
insert $\mp0.2$, $0$, $0$ at index $j'+1$ of $E_k^{l'}$, $I_k^{l'}$, $Q_k^{l'}$, respectively
insert $\mp0.1$, $0$, $0$ at index $j''+1$ of $E_k^{l'}$, $I_k^{l'}$, $Q_k^{l'}$, respectively
insert $\pm0.2$, $0$, $0$ at index $j''$ of $E_k^{l''}$, $I_k^{l'}$, $Q_k^{l'}$, respectively
Assign $x \leftarrow 2$ (Indicating 2 crossings have been introduced)
Update $u_k$ in $Z_k$ to reflect new signs:
Output: $j', l', j'', l''$ where $a=i_{j'}^{l'}$ and $c=i_{j''}^{l''}$

For Algorithm <ref>, we implement an index sorting procedure along each strand ($n=1,2$). First we determine where to induce this crossing (“Induction Decision” or Algorithm <ref>) by  forming three sub-arrays ($P$, $N$, and $B$) and comparing entries of $P$ and $N$ against $B$ to determine which indices to switch. We implement the following rule: overcrossings are moved to lower indices and undercrossings are moved to higher indices by either the application of an $R1$ or $R2$ Reidemeister move. This induces a crossing ($R1$) or pair of crossings ($R2$) at each index switch which is implemented in Algorithm <ref>.

The skein polynomial $Z$ is made up of individual skein terms $Z_k$, each containing a multi-curve $E_k$. Following the application of Algorithms 1-2, this multi-curve is made up of single curve components, $E_k^l$ with no-self intersections. Given a strand number $n$, we define unordered sets $P_n$, $N_n$, and $B_n$ from elements $(i_k^l)_j \in I_k^l$ based on the value of the corresponding elements $(e_k^l)_j \in E_k^l$ as follows:

\begin{equation*}
P_n=\{(i_k^l)_j\}, \text{such that } (e_k^l)_j=n,
\end{equation*}

\begin{equation*}
N_n=\{(i_k^l)_j\}, \text{such that } (e_k^l)_j=-n,
\end{equation*}

\begin{equation*}
B_n=\{(i_k^l)_j\}, \text{such that } |(e_k^l)_j|=n.
\end{equation*}

This amounts to partitioning the indices of the multi-curve $E_k$ into indices with a positive crossing on strand $n$ ($P_n$), indices with a negative crossing on strand $n$ ($N_n$), and a set that contains all indices along strand $n$ ($B_n$). Note that this partition is a single set across all single curve components $E_k^l$ (rather than a separate set for each component $l$).

The induced crossings are resolved by application of Algorithms <ref> and <ref>. The result of Algorithms <ref> and <ref> is the resolution of the skein into a linear combination of the fundamental curves given in REF. This is proved in Theorem <ref> below.

The strand sorting method of Algorithm <ref> produces only basis elements of the Kauffman Bracket Skein Module of the genus-2 handlebody: $x,y,z$ and the unknot.

The same method is applied to both strands of the base braid separately so consider elements along the first strand of the base braid ($|e^l_j|=1$). The result of the strand sorting method is that all overcrossings are moved to lower indices and undercrossings are moved to higher indices. This is equivalent to requiring the following condition along each strand: each skein term has the property that for all entries $e^l_j=1$ and all entries $e^{l'}_{j'}=-1$ the corresponding indices satisfy the relation $i^l_j<i^{l'}_{j'}$. We then realize that a skein with this requirement in the genus 2 handlebody is equivalent to a skein in the 2-punctured disk and all non-intersecting skeins in the 2-punctured disk are either $x$, $y$, $z$, or the unknot.

Consecutive index elements along a strand are considered pairwise from bottom to top. If the pair is a pair of overcrossings (or a pair of undercrossings), no change is made. If the pair is one overcrossing and one undercrossing, and the overcrossing has a higher index, a crossing is induced by a Reidemeister $R2$ move between the two strands, switching their positions and thus, their indices. If the overcrossing and undercrossing pair are in the same component, a single crossing is induced by a Reidemeister $R1$ move. The crossing(s) are resolved, producing either 2 or 4 skein terms in which the indices of the pair have been switched (the overcrossing now has a lower index than the undercrossing) or the same pair no longer exists (if the skein has been reduced). Continuing in this manner, all such pairs are switched until the condition above is satisfied.

§ CONCLUSION AND FUTURE WORK

We produced a method to resolve all crossings in a skein embedded in a genus-2 handlebody (Algorithm <ref> and <ref>) and then methodically induce crossings on the resulting linear combination of skeins (Algorithms <ref> and <ref>) to produce Laurent polynomial in the Pryzycki basis $x,y,z$ of the Kauffman Bracket Skein Module. This method, called "Resolve Original Diagram; Induce Crossings After (RODICA)" is named in memory of Dr. Rodica Gelca (Texas Tech University, Lubbock TX). This method provides a robust way to encode diagrammatic skein information and a means for the automation of all skein computations in a genus-2 handlebody. We further developed this method for the automation of skein computations in the complement of all 2-bridge knots, the topic of a follow on paper "Algorithms for Skein Manipulation and Automation of Skein Computations."

[1]
Vaughan F. R. Jones, “A polynomial invariant for knots via von Neumann
algebras,” Bulletin (New Series) of the American Mathematical
Society, vol. 12, no. 1, pp. 103 – 111, 1985.

[2]
Louis H. Kauffman, “State models and the Jones polynomial,” Topology, vol. 26, no. 3, pp. 395–407, 1987.

[3]
Edward Witten, “Quantum field theory and the jones polynomial,” Communications in Mathematical Physics, vol. 121, no. 3, pp. 351–399, 1989.

[4]
Nicolai Reshetikhin and Vladmir G Turaev, “Invariants of 3-manifolds via
link polynomials and quantum groups,” Inventiones mathematicae,
vol. 103, no. 3, pp. 547–598, 1991.

[5]
Vladimir G. Turaev, Quantum Invariants of Knots and 3-Manifolds.
Berlin, Boston: De Gruyter, 2016.

[6]
Christian Blanchet, Nathan Habegger, Gregor Masbaum, and Pierre Vogel,
“Topological quantum field theories derived from the kauffman bracket,”
Topology, vol. 31, pp. 685–699, 1992.

[7]
W. B. Raymond Lickorish, “The skein method for three-manifold
invariants,” Journal of Knot Theory and Its Ramifications,
vol. 2(02), pp. 171–194, 1993.

[8]
Jozef H Przytycki, “Skein modules of 3-manifolds,” Bulletin of the
Polish Academy of Sciences, pp. 91–100, 1991.

Skeins arise naturally in the study of knots and links as linear combinations of embedded (framed) circles in a 3-dimensional manifold identified modulo a set of locally defined "skein relations". Skein computations prove to be computationally intensive due to the exponential complexity of skein relations; try computing any knot polynomial with 5 or more crossings by hand will convince you why automation of these computations is so tempting!

In this paper, we have developed an automated method for performing skein computations in the complement of all 2-bridge knots. This family of knots is sufficiently restrictive as to make the automation process straightforward, but sufficiently rich in results. We use the Kauffman bracket skein relations as they relate to skeins constructed via the quantum group associated to SU(2) in Chern-Simons quantum field theory, but any skein relations can be used in the same why (with minor modifications).

It is important to stress that skein computations are not difficult to perform at any single stage. Given a fixed set of skein relations, one need only make the proper replacement, splitting one term into two, each with a different coefficient. Rather, it is the scale at which we wish to compute that is the primary motivation for this work. The manual organization required to perform such a computation by hand, carrying terms and signs correctly at each stage, is beyond what can be expected of most students, besides the fact that the diagrammatic nature of these computations is particularly demanding on the human brain, requiring a fairly involved mental manipulation of 3D geometric objects.

Finally, it is natural to ask whether it is worthwhile to perform such computations. The answer to this question is a resounding yes; amid skein theory's tedious computations exist surprisingly simple patterns and formulas. It often happens that a skein computation involving numerous terms simplifies to a relation with very few terms. Further, in the case of trefoil and 3-twist knots, the polynomial results always appear as Chebyshev polynomials, implying a deeper underlying structure at work, awaiting discovery. Our automation methods promote such discovery, revealing this structure through the generation of many examples in a reasonable timeframe.

In Section <ref>, we provide a brief introduction to skeins and skein computations, as well as a brief reference to the Kauffman Bracket Skein Module. In Section <ref>, the main body of this paper, we outline all the details regarding this novel method for the automation of skein computations in the genus 2 handlebody (Program Project) as well as in the complement of the 2-bridge knot (Program Resolve/Induce). In short, a multi-curve in the complement of a 2-bridge knot is projected to a specific handlebody embedded in the knot complement, and there manipulated based on indices assigned to each portion of the skein; strands which cross over in the given projection are moved to lower indices, and strands that cross under are moved to higher indices. The resolution of crossings introduced in this way produces basis elements of the Kauffman Bracket Skein Module with no more effort than brute force, although the steps are shown to be close to minimum thanks to thoughtful reductions after each stage of the resolution. In Section <ref>, we implement these methods in a Python script, applying our methods to several well-known applications. We conclude by providing this script to the reader for use in the automation of skein computations for the discovery of new relations!


11institutetext: Quantitative Methods Department, Universidad Pontificia
Comillas, Madrid, Spain 11email<EMAIL_ADDRESS>22institutetext:
Grupo Mente-Cerebro at ICS, Universidad de Navarra, Navarra, Spain 22email:
<EMAIL_ADDRESS>
# Optimizing Integrated Information with a Prior Guided Random Search
Algorithm

Eduardo C. Garrido-Merchán 1 Javier Sánchez-Cañizares2

(July 2022)

###### Abstract

Integrated information theory (IIT) is a theoretical framework that provides a
quantitative measure to estimate when a physical system is conscious, its
degree of consciousness, and the complexity of the qualia space that the
system is experiencing. Formally, IIT rests on the assumption that if a
surrogate physical system can fully embed the phenomenological properties of
consciousness, then the system properties must be constrained by the
properties of the qualia being experienced. Following this assumption, IIT
represents the physical system as a network of interconnected elements that
can be thought of as a probabilistic causal graph, $\mathcal{G}$, where each
node has an input-output function and all the graph is encoded in a transition
probability matrix. Consequently, IIT’s quantitative measure of consciousness,
$\Phi$, is computed with respect to the transition probability matrix and the
present state of the graph. In this paper, we provide a random search
algorithm that is able to optimize $\Phi$ in order to investigate, as the
number of nodes increases, the structure of the graphs that have higher
$\Phi$. We also provide arguments that show the difficulties of applying more
complex black-box search algorithms, such as Bayesian optimization or
metaheuristics, in this particular problem. Additionally, we suggest specific
research lines for these techniques to enhance the search algorithm that
guarantees maximal $\Phi$.

##  1 Introduction

Together with the Global Neuronal Workspace, Integrated Information Theory
(IIT) [5, 30, 23, 32, 31] is currently the most prominent approach to
obtaining a scientific theory of consciousness, i.e, a theory that can
identify the neural correlates of consciousness through some direct or
indirect measure. Difficulties in agreeing on experimental setups that may
discriminate between different theories have recently led to a sort of
adversarial collaboration between them. [18] Within said collaborative
framework, each theory provides the conditions for it to be experimentally
falsified. However, despite such partnership and the obtaining of promising
results in clinical tests, [9, 8, 15] reaching definitive answers may prove
more difficult than initially expected because of the multifactorial
complexity of the biological substrate of consciousness, namely the living
brain. [4]

The appeal of IIT for computation science and, in general, the AI community
relies on its straightforward definition of consciousness as a formal property
of any physical system, composed of nodes and interactions, that instantiates
it. Consciousness is maximally integrated information $\Phi$; [30] the latter
quantity being amenable to a strict mathematical definition in a network of
interconnected elements that can be thought of as a probabilistic causal
graph. Were said statement true, IIT would be providing not only a
quantifiable prediction for a system to be conscious but also a heuristic
approach to generate conscious systems —the holy grail of the AI program. Yet
it could also happen that IIT offers just a necessary, but not sufficient,
condition for the emergence of consciousness in a system. As a logical
consequence, a system implementing an artificial intelligence algorithm need
not be conscious whatsoever [12] and, hence, need not understand the meaning
of the problems that it solves. [10]

Interestingly, IIT derives its mathematical definition of consciousness from
five axioms that fully account for what, phenomenologically, being conscious
implies. To wit: The present experience of consciousness exists, is
informative, composed, integrated, and excludes other conscious experiences
[23]. Said axioms should naturally determine the emergence of a physical scale
in which the system, from its own perspective, presents a maximally integrated
structure of intrinsic information, in which the information of the whole is
maximally irreducible to the sum of information in subsets of any partition of
the system. Additionally, the physical state of the substrate instantiating
maximal $\Phi$, for the same token, presents a maximal cause-and-effect
composition (for the details of the link between measures of integrated
information and causality in an interconnected network, one may see [1, 2]).

Whereas IIT continues developing as a research program that strives for
finding the best definition of distance between informational contents [6] and
including in its framework many other features of consciousness, [14] it has
been criticized on diverse grounds [27, 24, 20]. More specifically for our
interests in this paper, some authors [21, 25] have raised the issue of the
so-called ’intrinsicality problem’ (IP), namely, the fact that consciousness
cannot be an intrinsic property of the system because maximal $\Phi$ crucially
depends on the possible existence of bigger values of $\Phi$ if the initial
system is appropriately linked to or embedded in larger systems. IP’s threat,
if true, could manifest that the causality involved in the constitution of the
system goes beyond the conceptualization of causality provided by IIT in its
current form. [25]

Even though the literature on the computation of $\Phi$ in different networks
is rapidly increasing and gaining insights on the logical architectures
favoring maximal $\Phi$ (basically those consisting of modular, homogeneous,
and specialized networks with feedback connections [23]), the prohibitively
large computation times for big networks, scaling with the number of nodes
($n$) as $\mathcal{O}(n53^{n})$, [17] make it necessary the search for
alternative heuristics and/or new methods for measuring $\Phi$ in larger
systems. [22, 29] Not much is known about the actual distribution of $\Phi$
over different network types and topologies. This paper thus explores the
benefit of random search algorithms in order to compute $\Phi$. If, according
to IIT’s creators ”the quantity and quality of consciousness are what they
are,” and consequently ”the cause–effect structure cannot vary arbitrarily
with the chosen measure of intrinsic information” [6], how maximal $\Phi$
varies with the number of systems nodes must be a relevant matter to get a
glimpse of how consciousness could eventually scale up in networks. Hence, our
paper studies the expected increase of maximal $\Phi$ with the number of
nodes, paving the way for gaining insights on the network’s architectures that
favor consciousness and assessing the seriousness of the IP, if IIT holds, in
future works. In such an endeavor, our paper will also show the limits of some
optimization techniques for IIT as the number of network nodes increases.

This paper is organized as follows. First, we provide a cursory definition of
integrated information, a problem definition of $\Phi$ optimization, and a
technical description of the methodology that we have implemented to optimize
$\Phi$. Most critically, we provide arguments to justify why our proposed
algorithm does not use smarter black-box global optimization techniques such
as Bayesian optimization. Then, we show in an experiments section empirical
evidence to support the claim that our method is able to optimize $\Phi$.
Finally, we present some concluding remarks and refer to possible research
lines to optimize $\Phi$ in further work.

##  2 Proposed Methodology

In this section, we begin formally defining integrated information and the
problem of optimizing integrated information in the space of transition
probability matrices. Then, we show how this problem cannot be solved with
common metaheuristics or techniques such as vanilla Bayesian optimization [11]
and introduce a random search algorithm to solve the problem that can be used
as a baseline for further research.

###  2.1 Integrated information

Throughout this paper, we use the definitions of integrated information
employed in [17]. Details regarding the differences between integrated
information ($\phi$) for a mechanism —a set of nodes and interactions in a
graph forming a concept— and conceptual integrated information ($\Phi$) for a
constellation of concepts forming a maximally irreducible conceptual
structure, can be consulted in [23]. Computations in our paper always refer to
$\Phi$.

We can define integrated information for a mechanism in state $s$ as

$\displaystyle\phi(s)=D[p(s_{0}\to s)\|p^{MIP}(s_{0}\to s)]\,,$ (1)

where $D$ is a measure of the distance between the probability distribution
$p$ of obtaining state $s$ as the cause (effect) of future (past) states
$s_{0}$ and the correspondent probability distribution $p^{IMP}$ when the
system has been optimally partitioned in its weakest links, i.e. decomposed
into those parts that are most independent (least integrated), providing a
minimum of effectively integrated information for the mechanism as a whole.

When one finds maximal $\phi$ for a mechanism in state $s$ in a system, one
has a concept. After determining the structure of concepts in a candidate
system, one is allowed to define the conceptual integrated information
($\Phi$) as $\Phi(C)=D(C\|C^{MIP})$, where $C$ refers to the structure of
concepts, $C^{MIP}$ refers to the optimal partition of the structure of
concepts that makes $\Phi(C)$ minimal, and $D$ is an appropriately adapted
measure of the distance between the corresponding probability distributions.
Different versions of IIT have employed different definitions of the distance
between probability distributions —the Kullback-Leibler divergence for $\phi$
[5], the Wasserstein Distance, also known as Earth Mover’s Distance, for both
$\phi$ and $\Phi$ [23], or the Intrinsic Difference [6]. As the PyPhi
algorithm favors the Earth Mover’s Distance as the measure of integrated
information [17], which quantifies how much two distributions differ by taking
into account the distance between system states, we will use it here too.

###  2.2 Problem definition

We now provide a formal definition of the global black-box optimization
problem that we need to solve. Concretely, we need to tackle the optimization
of a function $\Phi:\mathbb{R}^{2^{D},D}\to\mathbb{R}$, from a
D-$2^{D}$-dimensional transition probability matrix, with conditional
independencies, space to a real number.

$\displaystyle\omega^{\star}=\arg\max_{\omega\in\Omega}\Phi(\omega)\quad
s.t.\quad\omega\in\Upsilon\,,$ (2)

where $\omega^{\star}$ is the TPM that maximizes $\Phi(\cdot)$, $\omega$ is a
TPM of the D-dimensional TPM space, $\Omega\in\mathbb{R}^{2^{D},D}$ is the
D-dimensional TPM space, $\Upsilon$ is the space of TPM with conditional
independencies and $D\in\mathbb{N}^{+}$ is the number of dimensions of the
TPM, whose range in this problem is $\mathbb{N}$. Recall that the objective
function $\Phi$ is expensive to evaluate when the number of nodes of the TPM
is high. If $\omega\notin\Upsilon$, then $\Phi(\omega)=\emptyset$.
Consequently, $\mathbb{R}^{2^{D},D}$ is constrained to $\Upsilon$. This is
especially problematic, because of the curse of dimensionality. In particular,
as $D$ grows, the $\Omega_{D}$ space will contain higher $\Phi$ TPMs, which is
desirable. However, the number of valid TPMs will be lower with respect to all
the possible TPMs than in lower dimensional TPM spaces. Consequently, as the
evaluation of $\Phi$ is expensive, we have a tradeoff between risking to find
high $\Phi$ solutions in a higher dimensional TPM space where the volume of
valid TPMs is low or searching in a lower dimensional TPM space where we have
a high probability of obtaining a valid TPM but its $\Phi$ is potentially low.
More formally, if $\delta(\omega)$ is a delta function whose value is $1$ if
$\omega\in\Upsilon$, $D_{min}$ is the lowest dimensionality considered for the
number of nodes of the TPM, $D_{max}$ is the highest dimensionality considered
for the number of nodes of the TPM and $p(\omega|\Omega_{D})$ is the
probability of sampling each possible $\omega\in\Omega_{D}$, we have that:

$\displaystyle\int\delta(\omega)p(\omega|\Omega_{D_{max}})d\omega\leq...\leq\int\delta(\omega)p(\omega|\Omega_{D_{min}})d\omega\,.$
(3)

Moreover, as $D$ grows, the cost of evaluating $\Phi$ also grows. Hence, if
the budget is constrained by time, and not by iterations, we must be sure
about searching in a high dimensional $\Omega$ space. We tackle this issue and
define a new algorithm to solve this problem, among others that we further
explain, in the following subsection.

###  2.3 Proposed algorithm

As we have seen in the previous section, the global optimization problem of
the $\Phi$ function $\omega^{\star}=\arg\max_{\omega\in\Omega}\Phi(\omega)$ is
not smooth with respect to the smallest variations $\nabla\omega$ that could
be performed in the transition probability matrix space $\Omega$. For example,
in the case of binary nodes, switching just one node from zero to one can make
the $\Phi$ from $\mathbb{R}$ to $\emptyset$. In particular, this is
troublesome for some advanced black-box global optimization techniques, as we
describe in this section.

Bayesian optimization [11, 28] is the state-of-the-art class of methods to
solve black-box global optimization problems that iteratively recommend a
suggestion $\mathbf{x}$ making a tradeoff between exploration of promising
unknown regions of the search space $\mathcal{X}$ and potential good explored
solutions based on the predictions made by a probabilistic surrogate model as
a Gaussian process [26], which formally is defined as:

$\displaystyle\mathbf{x}^{\star}=\arg\min_{\mathbf{x}\in\mathcal{X}}f(\mathbf{x})\,,$
(4)

where $\mathbf{x}\in\mathbb{R}^{D}$ is a candidate point of the hypercube
$\mathcal{X}\subseteq\mathbb{R}^{D}$ and $f(\cdot)$ is the multi-modal
function to be optimized. In order to be considered a black-box global
optimization problem, $f$ must have no known analytical expression, and hence
we are not able to extract gradients; it can be noisy and it is expensive to
evaluate. Most critically, if the function $f(\cdot)$ is cheap to evaluate,
then it is recommended to use metaheuristics as genetic algorithms [16], that
also perform an exploration-exploitation tradeoff. The reason is that their
execution time to suggest recommendations $\mathbf{x}$ is faster than in the
case of Bayesian optimization, whose complexity is a cubic on the number of
previous recommendations $\mathcal{O}(N^{3})$, where $N$ is the number of
previous recommendations.

Data: $p(\omega),D_{min},D_{max},\mu,\epsilon,T$

Result: $\Phi^{\star},\omega^{\star}$

$i\leftarrow 0$;

$\Phi^{\star}\leftarrow 0$;

while  _$i\leq T$_ do

$j\leftarrow 0$;

$\Phi_{list}=[]$;

while  _$j\leq\epsilon$_ do

$D_{i}\leftarrow multinomial\\_sample(p(\omega),D_{min},D_{max})$;

$\omega_{i}\leftarrow uniform\\_TPM\\_sample(D_{i})$;

$\Phi^{i}=\Phi(\omega_{i})$;

if  _$\Phi^{\star}\geq\Phi_{i}$_ then

$\Phi^{\star}=\Phi_{i}$;

$\omega^{\star}=\omega_{i}$;

end if

$\Phi_{list}\leftarrow append\\_result(\Phi_{i},D_{i},\Phi_{list})$;

end while

$p(\omega)\leftarrow update\\_prior(p(\omega),\mu,\Phi_{list})$;

end while

Algorithm 1 Prior guided random search algorithm to optimize $\Phi$ on the
D-dimensional TPM space $\Omega$

The problem definition presented in the previous section involves the
computation of $\Phi$, which it is very expensive if the number of nodes of
the TPM is high. Consequently, this scenario does not present a good fit for
metaheuristic algorithms, as they require the evaluation of the objective
function to be cheap. Moreover, vanilla Bayesian optimization requires that
the function does not vary with respect to its parameters up to a certain
$K$-Lipschitz value. [13] Critically, in this problem the smallest variation
on any elements of the TPM matrix can make $\Phi$ jump from $\mathbb{R}$ to
$\emptyset$. Hence this assumption, required by Bayesian optimization, is
violated and would need to be specifically modified to tackle this particular
problem, whose complexity is even bigger as the space being optimized is
$\mathbb{R}^{2^{D},D}$, making Bayesian optimization also hard to use. Having
reviewed these alternatives, there also exists another reason why those
techniques may be inefficient. Concretely, if the $K$-Lipschitz value
assumption of Bayesian optimization is violated, the algorithm cannot perform
exploitation. Hence, any exploitative behavior on the D-dimensional TPM space
$\Omega$ would be inefficient. Consequently, the only alternative that remains
is to perform a prior guided random search algorithm, adapted to this
particular problem, and encode into the prior our knowledge about the $\Phi$
optimization problem. In particular, we enhance a random search procedure
instead of a grid search because random search significantly outperforms grid
search when the percentage of explained variance on the variable to optimize
is not uniform over the explanatory variables [7]. As assuming that the
variance of $\Phi$ is uniformly explained by all the nodes for all the
possible TPMs is a heavy assumption, we can not rely on uniform grid search,
as random search would be a smarter option in this case. Recall that random
search outperforms grid search when some variables of the input space explain
the objective function better than others. Consequently, we decided to enhance
the random search procedure. The prior guided random search illustrated on
Algorithm 1 will obtain $\Phi^{\star}$, which is an approximation to
$\Phi^{best}$, that is, $\Phi^{\star}\approx\Phi^{best}$, which could only be
found using an exhaustive search throughout all the TPM D-dimensional space
$\Omega$. Where, in Algorithm 1, $\mu$ is the adaptation rate of $p(\omega)$,
the multinomial a priori distribution of the probability of obtaining
$\Phi^{\star}$ in a certain $D$-dimensional space, $\epsilon$ is the number of
iterations that are needed to update $p(\omega)$, $T$ is the budget of total
iterations that the user can afford, $D_{min}$ is the lower dimension of the
TPM D-dimensional space and $D_{max}$ is the higher dimension of the TPM
D-dimensional space. The sampling functions define a probabilistic graphical
model where the TPM is modeled as a random variable that depends on the
multinomial random variable that encodes the probability distribution assigned
to the dimensions of the D-dimensional TPM space. More formally, we obtain
$\omega_{i}$ sampling from the following joint probability distribution:

$\displaystyle\omega=uniform\\_matrix(\omega\mid D)\quad
multinomial(D\mid\theta,D_{min},D_{max})\,.$ (5)

Sampling from the previous joint distribution of both random variables we can
sample TPMs conditioned to the probability of obtaining a TPM belonging to the
space of valid matrices $\Upsilon$. The initial state of the TPM necessary to
compute $\Phi$ is tested in an iterative fashion until a valid state is found.

At the end of a batch of $\epsilon$ iterations, we update our beliefs about
the $\theta$ parameter values of our multinomial distribution
$p(\omega|\theta)$ up to a learning rate value $\mu$. We do so in a Bayesian
learning fashion, basically defining a likelihood function
$f(\Phi_{list},\mu)$ that penalizes the dimensions whose results are worse and
rewards those whose values are better up to $\mu$. This likelihood function is
multiplied by the prior beliefs $p(\omega|\theta)$ and the posterior is
obtained by normalizing the previous operation with a normalization constant
$\sum_{D_{min}}^{D_{max}}p(\omega|\theta)f(\Phi_{list})$ to ensure that the
resulting distribution is a valid multinomial distribution and the prior is
conjugate. Finally, we smooth the maximum and the minimum of the vector of
probabilities by a quantity $\kappa=0.02$ that could be also modified by the
user. The specific parametric likelihood function $f(\Phi_{list},\mu)$ can be
set by the user. Concretely, we just provide as default to linearly penalize
or reward each dimension in $0.2\mu r$ points where $r$ is the ranking of each
dimension $D$ that consists of a $[-1,1]$ linear space vector. More formally,
if $p(\omega|\theta)_{new}$ is the posterior distribution obtained after the
batch of $\epsilon$ iterations, then, to obtain this distribution we perform
the following update:

$\displaystyle
p(\omega|\theta)_{new}=\frac{f(\Phi_{list},\mu)p(\omega|\theta)}{\sum_{D_{min}}^{D_{max}}f(\Phi_{list},\mu)p(\omega|\theta)}$
(6)

Lastly, we remark that the values of the set of prior beliefs $\theta$, the
learning rate $\mu$, and the batch size of iterations $\epsilon$ are hyper-
parameters that the user can set and will modify the behavior of the
algorithm. Additionally, they also allow the user to set the value of
$D_{min},D_{max}$ to define the TPM search space $\Omega$. In the next
section, we illustrate the algorithm performance, being compared with a grid
search method and a random search method as baselines.

##  3 Experiments

In order to estimate the integrated information of a transition probability
matrix, we use the PyPhi tool implemented on Python3 [17]. In the computation
of $\Phi$, to measure distances between probability distributions, we use the
Earth Mover’s Distance [3], which is equivalent to the Hamming Distance if the
nodes of the graph are binary. We split the experiment section into two parts.
First, we show the results obtained by the prior guided random search
algorithm with respect to other baselines, just to illustrate how we can
obtain graphs whose $\Phi$ score is bigger algorithmically. Driven by the
results of these experiments and because of the size of the TPM space that
makes it computationally unfeasible to explore, we append another section
where we try to statistically make inference about $\Phi$ as the number of
nodes grows. We exhaustively search for an initial state for the TPM and
retrieve the first feasible state for a randomly sampled TPM, making a point
estimation about the mean $\Phi$ of the TPM with the first feasible state
found, for computational reasons. All the code, written on Python3, and the
experiments are uploaded into the Github repository
(https://github.com/EduardoGarrido90/iit_opt).

###  3.1 Prior Guided Random Search experiments

The following experiments model the obtained $\Phi^{\star}$ by the various
methods as a Gaussian random variable with mean $\mu$ and standard deviation
$\sigma$. All the experiments are repeated $25$ times to infer their
respective values of $\mu$ and standard deviation $\sigma$. Namely, we compare
our proposed prior guided random search method with respect to a grid search
method and a random search method. The random search method samples uniformly
both the biggest dimension of the TPM and the values of the matrix. However,
the grid search method uniformly splits the dimensions of the TPM across all
the space range, and the binary TPMs are selected by obtaining the binary
representations of a linear space of integers from $0$ to $2^{n}$ where $n$ is
the number of nodes which are the columns of the matrices.

First, we begin with the following experiment, to verify whether all the
methods work and provide valid $\Phi$ results. We perform 25 repetitions of
the experiment with each method to ensure statistical significance. We set the
lower dimension, or the number of nodes, to 3 and the higher dimension to 4.
We execute a total of 50 iterations with $\epsilon=5$ for the prior guided
search algorithm. We set $\mu=0.1$ and the prior distribution values on the
dimensions to $p(\omega)=[0.2,0.8]$. With these values, we obtain the $\Phi$
results across the iterations shown in Figure 1. The best $\Phi$ found in the
search was $3.034082$.

Figure 1: Average results of $\Phi$ with the toy experiment configuration of
the Prior Guided Random Search method (blue) compared to the Grid (green) and
Random search (orange) baselines. Better seen in color.

As we can see in Figure 1, the prior guided random search method, plotted in
blue, outperform the baselines. The results are statistically significant as
the plotted results are the means of $25$ different experiments. As the prior
guided random search method performs almost twice as better as the grid search
method, a statistical hypothesis testing method will provide a p-value lower
than $.01$ to reject the null hypotheses that the differences are random,
accepting the alternative hypothesis that both methods have different
distributions. Concretely, random search also outperforms grid search. And
lastly, the more interesting conclusion is that the prior guided random search
method outperforms random search. Regarding grid search, we have found that
the restriction of the matrices of the linear grid seems to provide bad
results, resulting mostly on invalid TPMs and not being a method that is worth
using to optimize $\Phi$ in the $\Omega$ space. Hence, for the next
experiments, and also for computational reasons due to the complexity of
$\Phi$ with respect to the number of nodes, we will only compare the prior
guided random search method with the random search method.
We have performed a little experiment, which we will call Graphical
Experiment, of one repetition just to show graphically the obtained transition
probability matrix, the connectivity matrix, and the initial state. This
experiment involves 100 repetitions with the minimum dimension of nodes being
$3$ and the maximum number of nodes $4$. The initial prior is $[0.3,0.7]$. We
obtained, with the prior guided random search algorithm, a graph with
$\Phi=2.6823$, with a final prior with value $[0.12,0.88]$, the connectivity
matrix shown on Figure 2, the transition probability matrix shown on Figure 3
and the state shown on Figure 4.

Figure 2: Connectivity matrix associated to the best result obtained in the
Graphical Experiment. Figure 3: Transition probability matrix associated to
the best result obtained in the Graphical Experiment. 0 is painted in purple
and 1 in yellow. Figure 4: Initial state associated with the best result
obtained in the Graphical Experiment. 0 is painted in purple and 1 in yellow.

In the following experiment, we want to verify whether the behavior
differences in the method generalize to more than 2 dimensions in the space of
TPMs, this time, from 2 nodes to 5. With a prior distribution with weights
$p(\omega)=[0.1,0.1,0.3,0.5]$ and only $5$ repetitions to save computational
resources and $50$ iterations. We obtain the results plotted in Figure 5 where
we see how the Prior Guided Random Search algorithm is also outperforming the
random search method.

Figure 5: Average results of $\Phi$ with the small experiment configuration of
the Prior Guided Random Search method (blue) compared to the Random search
(orange) baseline. Better seen in color.

We can see, once again, that both searches are able to find better TPMs
according to the $\Phi$ measure once the iterations are executed and that
prior guided random search once again outperforms random search, not only on
average but in absolute value. In particular, the maximum $\Phi$ found by the
prior guided random search method is $3.1681$ and the best result obtained by
random search is $2.0832$.

As we have stated previously, the complexity of evaluating $\Phi$ with respect
to the nodes is high, in particular, $\mathcal{O}(n53^{n})$, which must be
multiplied by $2^{n}$ in these searches as we look for a feasible state for a
randomly sampled TPM. Consequently, and having shown empirical evidence of the
usefulness of the random searches method and how our proposed method even
outperforms random search for this specific task, our last experiment only
involves $1$ repetition of each method, for computational reasons. Being
unable to test $\Phi$ in TPMs with $6$ or more nodes due to the complexity of
the algorithm, in this experiment, we will check whether better results of
$\Phi$ are obtained if the number of nodes and iterations is bigger. The
minimum dimension is now $4$ and the maximum dimension is $5$. The prior
distribution is $[0.2,0.8]$ and the number of iterations is $500$, with
$\mu=1$ and $\epsilon=5$, obtained the results plotted in Figure 6.

Figure 6: Average results of $\Phi$ with the small experiment configuration of
the Prior Guided Random Search method (blue) compared to the Random search
(orange) baseline. Better seen in color.

We can see how better results are obtained if we have a bigger budget of
iterations or time. Moreover, this experiment also shows that the Prior Guided
Random Search algorithm outperforms Random search once the prior distribution
fits the optimal number of nodes. As the $\Omega$ space is very sparse, a much
higher number of iterations would be required to obtain a better approximation
to the optimal $\Phi$. More critically, for further work, with more
computation power that would enable searches for graphs of more than $6$
nodes, the space would even be more sparse, so much more computation would be
needed to find graphs with bigger $\Phi$. However, we have shown with these
experiments how our algorithm would adapt to the sparsity of every dimension
and keep searching for solutions in the number of nodes that represents the
better trade-off between the probability of finding better $\Phi$ and the
sparsity of the dimension that makes difficult to find a feasible $\Phi$.

###  3.2 Experiments to explore the graph space and make statistical inference
of $\Phi$

The complexity of the space of transition probability matrices makes it
computationally intractable to be explored exhaustively. Moreover, in order to
measure $\Phi$, we must not only take into account the matrix that encodes the
transitions but also the connectivity matrix, that represents the connections,
or edges, of the graph. Consequently, the complexity, or size measure, of the
space is even higher.

However, we can frame this problem into the classical statistical inference
framework. The population is the total of possible combinations for a given
number of graphs and we can extract a random sample from it. The goal is to
make inferences about the mean $\Phi$ parameter $\mu$ associated with each
number of nodes. Then, we can perform statistical hypothesis testing to
determine whether the difference of $\Phi$ is statistically significant when
$N$ grows. By doing this procedure, we can validate, with the previous set of
experiments, that our algorithm is going to work even better when the number
of nodes grows, assuming that the same inference will hold when the number of
nodes grows, which we can not empirically test since the complexity of
computing $\Phi$ when the number of nodes grows is exponential.

The first toy experiment involves 100 graphs of 3 nodes and 100 graphs of 4
nodes. We obtain an empirical $\Phi$ mean of $0.0517$ in the case of 3 nodes
and, in the case of 4 nodes, an empirical mean of $0.1711$. The confidence
intervals at 95% are $[0.0212,0.0822]$ in the case of 3 nodes and
$[0.1093,0.2328]$ in the case of 4 nodes. We also perform a 2 sample t-test
whose null hypothesis is that both population means are equal and the
alternative hypothesis is that they are different. We obtain a t-statistic of
$-3.38$ whose associated p-value is $0.0008$. At a confidence level of
$\alpha=0.01$ we can reject the null hypothesis that both means are equal,
hence accepting the alternative hypothesis that they are different, being the
$\Phi$ population mean of graphs of 4 nodes bigger than the $\Phi$ population
mean of graphs of 3 nodes. Also, the percentage of infeasible solutions found
in the case of $3$ nodes has been $72.6775\%$ and in the case of $4$ nodes has
been $84.7094\%$, providing empirical evidence of Eq. 2. To test consistency,
we performed the same experiment with a population of $200$ in both groups. In
this case, we obtained an even bigger, in absolute terms, a t-test statistic
of value $-4.6420$ with associated p-value $4.69e-6$, rejecting the null
hypothesis. The obtained empirical means are $0.0569$ and $0.1847$, pretty
similar to the ones obtained in the previous experiment with similar
confidence intervals.

The second experiment involves a population of 100 graphs of 4 nodes and 100
graphs of 5 nodes. We obtain an empirical $\Phi$ mean of $0.1432$ in the case
of 4 nodes and, in the case of 5 nodes, an empirical mean of $0.3267$. The
confidence intervals at 95% are $[0.0822,0.2043]$ in the case of 4 nodes and
$[0.2217,0.4317]$ in the case of 5 nodes. We also perform a 2 sample t-test
whose null hypothesis is that both population means are equal and the
alternative hypothesis is that they are different. We obtain a t-statistic of
$-2.94$ whose associated p-value is $0.0036$. At a confidence level of
$\alpha=0.01$ we can reject the null hypothesis that both means are equal,
hence accepting the alternative hypothesis that they are different, being the
$\Phi$ population mean of graphs of 5 nodes bigger than the $\Phi$ population
mean of graphs of 3 nodes. Also, the percentage of infeasible solutions found
in the case of $4$ nodes has been $84.98\%$ and in the case of $5$ nodes has
been a $91.84\%$

##  4 Philosophical discussion

This work presents an improvement on how to compute $\Phi$ for arbitrary
networks and TPMs. Whereas the present calculations still work with toy models
with a maximum of 6 nodes, the obtained results with a guided random method
strongly suggest that the latter outperforms other methods to explore optimal
TPMs as the number of nodes increases. In this sense, it reveals a promising
procedure for further progress in the time computation of $\Phi$ and,
eventually, for the ascertainment of the validity of IIT as the scientific
theory of consciousness in realistic networks. We aim at answering questions
like these: Can we quantitatively provide a method to generate consciousness?
May this methodology be adapted to deep neural networks if the assumptions of
IIT were true, providing deep learning machine consciousness architectures
[19]?

On the other hand, one should keep in mind the intrinsic limitations of the
PyPhi algorithm, which obviously affect our results. According to its
creators, “[T]he analysis can only be meaningfully applied to a system that is
Markovian and satisfies the conditional independence property. These are
reasonable assumptions for the intended use case of the software: analyzing a
causal TPM derived using the calculus of perturbations. However, there is no
guarantee that these assumptions will be valid in other circumstances, such as
TPMs derived from observed time series (e.g., EEG recordings).” [17]

The ideas of the last paragraph allow us to dive into additional philosophical
problems related to the overall IIT procedures. One of those problems is the
origin of networks. In the case of artificial networks (1), the unity of the
network is artificially created by network designers, avoiding external
influences that can disrupt the workings of the network. In the case of
natural networks (2), however, one may wonder who or what causes the unity of
the system. Moreover, whereas in case (1) possible consciousness as defined by
IIT remains extrinsic, as a bigger maximal $\Phi$ can be reached by adding new
nodes to the initial network, it is no longer clear —indeed very
controversial— that the same procedure works in case (2). Case (2) could be
the only scenario in which, in keeping with IIT tenets, consciousness is
intrinsic. But if that is the case, IP still holds for case (2).

Last but not least, our procedure to obtain bigger $\Phi$ as the number of
nodes is increased and the system’s architecture becomes more complex might
not reflect the ontology of case (2). The bottom-line problem is that an
optimization procedure might just be an epistemic tool to obtain the magnitude
which maximizes an extant, meaningful, quantity. However, the procedure itself
need not unravel the ontological processes underlying the reality of said
quantity.

##  5 Conclusions and further work

Throughout this paper, we have dealt with integrated information theory as a
relatively novel approach to quantifying truly integrated information in a
network of interacting nodes. Beyond its technical interest in the field,
research in IIT is exploding as it aims to become a scientific theory of
consciousness by validating its axioms in artificial and natural networks.
Nevertheless, at present, procedures to compute $\Phi$ remain too lengthy.
Even if the creation of the PyPhi software has provided IIT researchers with a
magnificent tool, computation times increase too quickly with the number of
nodes in a network. Yet, studying how maximal $\Phi$ varies when new nodes and
interactions are added to the system may prove crucial to see how, allegedly,
consciousness scales up and whether the intrinsicality problem presents a true
threat to IIT.

Our paper aims to provide a new method to evaluate maximal $\Phi$ with the
help of state-of-the-art optimization procedures of network quantities.
However, as stated in the technical sections, the requirements of conditional
independence for the TPMs charge our procedure with a burden that may not be
eased out while keeping all its computational strengths. For the time being,
ours remain a toy model that seems promising to pursue in future works.

This work thus opens up several lines of research. To enhance the behavior of
the method, we hypothesize that several actions could be taken. First, it may
be possible to optimize on a smooth, continuous space of conditional TPMs to
use Bayesian optimization or genetic algorithms. This approach would need to
formulate an analytical transformation from the space of TPMs into a smooth
and continuous space that can be modeled with a Gaussian process, that is the
probabilistic surrogate model of Bayesian optimization. By doing so, we can
exploit the exploitation-exploration behavior of Bayesian optimization, which
may outperform the current prior-guided random search algorithm. If this
transformation cannot be eventually found, we can extract the best-observed
result by the prior-guided random search algorithm and execute a combinatorial
optimization algorithm with the best-observed result as the initial point.
This may also provide good results if $\Phi$ varies smoothly as a result of
varying the TPM or the CM in a single element. Lastly, as we have seen, $\Phi$
is computationally costly. An interesting line of research would be to find a
cheaper approximation of $\Phi$ with respect to a particular graph with a
regression graph machine learning algorithm. In order to fit this algorithm it
would be necessary to create a huge dataset of graphs annotated by its $\Phi$
values and study whether a regression graph machine learning algorithm is able
to accurately approximate $\Phi$.

## References

  * [1] Albantakis, L., Marshall, W., Hoel, E., and Tononi, G. What Caused What? A Quantitative Account of Actual Causation Using Dynamical Causal Networks.  Entropy 21, 5 (may 2019), 459.
  * [2] Albantakis, L., and Tononi, G. Causal Composition: Structural Differences among Dynamically Equivalent Systems.  Entropy 21, 10 (oct 2019), 989.
  * [3] Andoni, A., Indyk, P., and Krauthgamer, R. Earth mover distance over high-dimensional spaces.  In SODA (2008), vol. 8, pp. 343–352.
  * [4] Arsiwalla, X. D., and Verschure, P. Measuring the Complexity of Consciousness.  Frontiers in Neuroscience 12, JUN (jun 2018), 1–6.
  * [5] Balduzzi, D., and Tononi, G. Integrated Information in Discrete Dynamical Systems: Motivation and Theoretical Framework.  PLoS Computational Biology 4, 6 (jun 2008), e1000091.
  * [6] Barbosa, L. S., Marshall, W., Albantakis, L., and Tononi, G. Mechanism Integrated Information.  Entropy 23, 3 (mar 2021), 362.
  * [7] Bergstra, J., and Bengio, Y. Random search for hyper-parameter optimization.  Journal of machine learning research 13, 2 (2012).
  * [8] Bodart, O., Gosseries, O., Wannez, S., Thibaut, A., Annen, J., Boly, M., Rosanova, M., Casali, A. G., Casarotto, S., Tononi, G., Massimini, M., and Laureys, S. Measures of metabolism and complexity in the brain of patients with disorders of consciousness.  NeuroImage: Clinical 14 (2017), 354–362.
  * [9] Casali, A. G., Gosseries, O., Rosanova, M., Boly, M., Sarasso, S., Casali, K. R., Casarotto, S., Bruno, M.-A., Laureys, S., Tononi, G., and Massimini, M. A Theoretically Based Index of Consciousness Independent of Sensory Processing and Behavior.  Science Translational Medicine 5, 198 (aug 2013), 198ra105–198ra105.
  * [10] Garrido-Merchán, E. C., and Blanco, C. Do artificial intelligence systems understand?  arXiv e-prints (2022), arXiv–2207.
  * [11] Garrido Merchán, E. C., et al. Advanced methods for bayesian optimization in complex scenarios.
  * [12] Garrido Merchán, E. C., and Lumbreras, S. On the independence between phenomenal consciousness and computational intelligence.  arXiv e-prints (2022), arXiv–2208.
  * [13] Hager, W. W. Lipschitz continuity for constrained processes.  SIAM Journal on Control and Optimization 17, 3 (1979), 321–338.
  * [14] Haun, A., and Tononi, G. Why does space feel theway it does? Towards a principled account of spatial experience.  Entropy 21, 12 (2019).
  * [15] Haun, A. M., Oizumi, M., Kovach, C. K., Kawasaki, H., Oya, H., Howard, M. A., Adolphs, R., and Tsuchiya, N. Conscious Perception as Integrated Information Patterns in Human Electrocorticography.  eneuro 4, 5 (sep 2017), ENEURO.0085–17.2017.
  * [16] Kumar, M., Husain, D., Upreti, N., Gupta, D., et al. Genetic algorithm: Review and application.  Available at SSRN 3529843 (2010).
  * [17] Mayner, W. G., Marshall, W., Albantakis, L., Findlay, G., Marchman, R., and Tononi, G. Pyphi: A toolbox for integrated information theory.  PLoS computational biology 14, 7 (2018), e1006343.
  * [18] Melloni, L., Mudrik, L., Pitts, M., and Koch, C. Making the hard problem of consciousness easier.  Science 372, 6545 (may 2021), 911–912.
  * [19] Merchán, E. C. G., and Molina, M. A machine consciousness architecture based on deep learning and gaussian processes.  In International Conference on Hybrid Artificial Intelligence Systems (2020), Springer, pp. 350–361.
  * [20] Merker, B., Williford, K., and Rudrauf, D. The integrated information theory of consciousness: A case of mistaken identity.  Behavioral and Brain Sciences 45, May (may 2022), e41.
  * [21] Mørch, H. H. Is Consciousness Intrinsic? A Problem for the Integrated Information Theory.  Journal of Consciousness Studies 26, 1-2 (2019), 133–162(30).
  * [22] Nilsen, A. S., Juel, B. E., and Marshall, W. Evaluating approximations and heuristic measures of integrated information.  Entropy 21, 5 (2019).
  * [23] Oizumi, M., Albantakis, L., and Tononi, G. From the Phenomenology to the Mechanisms of Consciousness: Integrated Information Theory 3.0.  PLoS Computational Biology 10, 5 (may 2014), e1003588.
  * [24] Pautz, A. What is the integrated information theory of consciousness? A catalogue of questions.  Journal of Consciousness Studies 26, 1-2 (2019), 188–215.
  * [25] Sánchez-Cañizares, J. Formal Causation in Integrated Information Theory: An Answer to the Intrinsicality Problem.  Foundations of Science 27, 1 (mar 2022), 77–94.
  * [26] Schulz, E., Speekenbrink, M., and Krause, A. A tutorial on gaussian process regression: Modelling, exploring, and exploiting functions.  Journal of Mathematical Psychology 85 (2018), 1–16.
  * [27] Searle, J. R. Can Information Theory Explain Consciousness?, 2013.
  * [28] Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and De Freitas, N. Taking the human out of the loop: A review of bayesian optimization.  Proceedings of the IEEE 104, 1 (2015), 148–175.
  * [29] Toker, D., and Sommer, F. T. Information integration in large brain networks.  PLOS Computational Biology 15, 2 (feb 2019), e1006807.
  * [30] Tononi, G. Consciousness as Integrated Information: a Provisional Manifesto.  Biological Bulletin, December (2008), 216–242.
  * [31] Tononi, G. The Integrated Information Theory of Consciousness.  In The Blackwell Companion to Consciousness. Wiley, Chichester, UK, apr 2017, pp. 243–256.
  * [32] Tononi, G., Boly, M., Massimini, M., and Koch, C. Integrated information theory: from consciousness to its physical substrate.  Nature Reviews Neuroscience 17, 7 (2016), 450–461.