The True Centroid of a Polygon

 My apologies for taking a break from the quadratic programming thread. I will be coming back to it. But for now I thought I would record a method for finding the true centroid of a planar polygon, given as a sequence of vertices that traverse the boundary in order. A quick and dirty approximation is to simply take the barycenter of the vertices. This is often good enough, and for convex polygons, is guaranteed to lie within the interior of the polygon, but consider the following example. Take any polygon and then split one of its vertices into two vertices connected by a very short edge. This new polygon should have approximately the same centroid, but the barycenter of its vertices changes quite a bit since the split vertex now exerts twice the pull on the centroid than it did before.

Luckily there is a very simple way to get the centroid, at least for planar polygons. I believe the formulas should generalize to higher dimensional embeddings of polygons, but have not given the matter sufficient thought. We start with Green's Theorem 

$$\iint_R \frac{\partial M}{\partial x}-\frac{\partial L}{\partial y}\,dx\,dy = \oint_{\partial R} L\,dx+M\,dy$$

for integrating a 1-form on the boundary of a region. 

As a warm-up, let $M=x$ and $L=0$. Then we discover that the area of the polygon is the integral

$$\oint_{\partial R} x\,dy$$ along its boundary. This integral can be evaluated on each line segment comprising the boundary to get a simple combinatorial formula.

Let $p_0=\left(\begin{matrix} x_0\\y_0\end{matrix}\right)$ and $p_1=\left(\begin{matrix} x_1\\y_1\end{matrix}\right)$. Then the path connecting them is $(1-t)p_0+tp_1$ for $0\leq t \leq 1$.  

Then $dy = (y_1-y_0)\,dt$ and $x= (1-t)x_0+tx_1$, so we get

$$\int_0^1 ((1-t)x_0+tx_1)(y_1-y_0)\,dt= \frac{1}{2}(x_1+x_0)(y_1-y_0).$$

This gives us a beautiful formula for area. If the polygon has $n$ points $p_1,\ldots,p_n$ we have

$$A=\frac{1}{2}\sum_{i=1}^n(x_i+x_{i+1})(y_{i+1}-y_i)$$

where the indices are read modulo $n$. $x_{n+1}=x_1, y_{n+1}=y_1$. As an exercise see if you can find two more similar formulas for area with different choices of $L$ and $M$ above. These types of formulas are the basis for mechanical planimeters which calculate area by tracing around the perimeter.

Returning to the centroid, we need to calculate the average position of both $x$ and $y$:

$$\overline{x} = \frac{1}{A}\iint_R x\,dx\,dy,\,\,\,\overline{y} = \frac{1}{A}\iint_R y\,dx\,dy$$

Tackling the $x$ coordinate, let $M=\frac{x^2}{2}$ and $L=0$. Then we have

$$\iint_R x\,dx = \sum_{i=1}^n\int_0^1 \frac{((1-t)x_i)^2+t x_{i+1}}{2}(y_{i+1}-y_i)\,dt$$

$$\iint_R x\,dx = \frac{1}{6}\sum_{i=1}^n (x_i^2+x_ix_{i+1}+x_{i+1}^2)(y_{i+1}-y_i)$$

In a similar way we can derive the formula

$$\iint_R y\,dy = \frac{1}{6}\sum_{i=1}^n (y_i^2+y_iy_{i+1}+y_{i+1}^2)(x_{i}-x_{i+1}).$$

It does not feel to me like these formulas are unique to planar embeddings of polygons, but the technique of proof, using Green's theorem, will not work for higher ambient dimensions.

Minimizing a convex quadratic with equality constraints

In our last few posts we tackled the question of finding a global minimum for a function $f(x)=c^Tx+\frac{1}{2} x^TCx.$ The heart of the algorithm was finding an orthonormal basis for $\mathbb R^n$ with respect to the inner product $\langle x,y\rangle =x^TCy$.

Now suppose we have a set of equality constraints $Ax=b$. This defines an affine subspace of $\mathbb R^n$ and if we can find an orthonormal basis of the same plane translated to the origin $Ax=0$, then we can just use the same algorithm as before.

There is a neat trick for doing this. First, assume that the rows of $A$ $a_1\ldots,a_m$ are linearly independent. (We will remove that restriction later.) Now add extra rows to $A$ so that its rows form a basis of $\mathbb R^n$, to get a matrix we call $D$. 

Claim: The $m+1,\ldots, n$ columns of $D^{-1}$ form a basis for the subspace $Ax=0$.

Proof:

Suppose $Ax=0$. Then $Dx=(0,\ldots,0,*,\ldots,*)^T$ is $0$ in the first $m$ coordinates. Applying $D^{-1}$ to both sides, we have $x=D^{-1}(0,\ldots,0,*,\ldots,*)^T$, which is in the column space of $D^{-1}$ spanned by the last $n-m$ columns. $\Box$

We proceed by doing Gram-Schmidt orthogonalization to the last $n-m$ columns. So to begin, normalize the $m+1$st column and subtract off projections of this column from all other columns, including the initial $m$. These column operations correspond to similar row operations on $D$, so that we are adding a multiple of the $m+1$st row of $D$ to the vectors spanning the subspace $Ax=b$ (see below for a simple explanation of this phenomenon). Thus by the above claim, the last $n-m-1$ vectors of the new $D^{-1}$ are orthogonal to the subspace spanned by the first $m+1$ columns.

 Repeating this process with the remaining columns as we move to the right, we get a collection of orthonormal vectors spanning our subspace.

 Thus we can use the algorithm developed in the unconstrained case with very little modification, but with a little more work needed in the initialization. 

• First, you need to be able to identify a point $x_0$ satisfying the constraints: $Ax_0=b$.
• Second, you need to add rows to $A$ to get an invertible matrix $D$, and start the algorithm off with $GS=D^{-1}$ rather than just the identity.
• Also as a note, you have to make sure to apply the matrix update rules to the first $m$ columns at each step, not just the columns to the right of the currently active one.
  

 I hope to write this out as a detailed algorithm in my next post, but for now, let us return to the question of how row/column operations affect the inverse.

Lemma: An elementary row operation on a matrix $M$ corresponds to an elementary column operation on $M^{-1}$.

Proof: An elementary row operation is realized by multiplying by an elementary matrix on the left $EM$. Taking the inverse we have $M^{-1}E^{-1}$. The inverse of an elementary matrix is an elementary matrix, and right multiplication is an elementary column operation. $\Box$

Unconstrained quadratic optimization: the algorithm

 We finally have all the tools in place to describe a step algorithm for minimizing a (convex) quadratic objective function without any constraints.

• Minimize $f(x)= c\cdot x+\frac{1}{2}x^TCx$ where $C$ is a positive semidefinite matrix.

We start with any seed value $x_0\in \mathbb R^n$. Recall that the modified Gram-Schmidt process starts with a set of independent vectors and gradually replaces them with an orthonormal set. We will start with the standard basis $e_1,\ldots,e_n$, recorded as the columns of a matrix $GS_0=I$. We use a vector $J\in\{Y,N\}^n$ to keep track of which columns have been orthonormalized. So to begin we have $J=\{N,\ldots,N\}$, and as we work on a column to create a conjugate direction, we add $Y$ to that index of $J$. Set $g_0=\nabla f(x_0)=c+Cx_0$ and $j=0$ to initialize, where $j$ represents the iteration of the algorithm.

• Step 1: Computation of search direction $s_j$: Let $GS_j=[c_{1j},\ldots, c_{nj}]$ and $J_j=\{\alpha_{1j},\ldots,\alpha_{nj}]$. If each component of $J_j$ is equal to $Y$, then stop with optimal solution $x_j$. Otherwise, compute the smallest index $k$ such that $$|g_j^Tc_{kj}|= \max\{|g_j^Tc_{ij}|\,|\, \forall i \alpha_{ij}=N\}.$$ If $g^T_jc_{kj}=0$, stop with optimal solution $x_j$. Otherwise set $s_j=\operatorname{sgn}(g^T_jc_{kj})c_{kj}.$ Go to step 2.
• Step 2: Computation of step size $t_j$: Compute $\langle s_j,s_j\rangle$. If it is $0$, then stop. The objective function is unbounded from below. Otherwise set: $$t_j=\frac{g_j^Ts_j}{\langle s_j,s_j\rangle} $$ Go to step 3.
• Step 3: Updating data: Set
• $x_{j+1}:=x_j-t_js_j$.
• $g_{j+1} := \nabla f(x_{j+1}) = c+Cx_{j+1}$
• $f(x_{j+1}) = c^Tx_{j+1} +\frac{1}{2}x_{j+1}^TCx_{j+1}$
• $d_j = s_j/\sqrt{\langle s_j,s_j\rangle}$
• $GS_{j+1}$ is formed from $GS_j$ by 
• replacing the $k$th column $c_k$ by $c_k/(d_j^Tc_k)$ 
• replacing all columns $c_i$ where $\alpha_{ij}=N$ by $c_i-\frac{d_j^Tc_i}{d_j^Tc_k}c_k$
• leaving columns $c_i$ where $\alpha_{ij}=Y$ unchanged.
• $J_{j+1}$ is formed from $J_j$ be setting $\alpha_{k}$ to $Y$ and leaving everything else unchanged.
• Increment $j:=j+1$ and go to step 1!

To summarize the algorithm works as follows. Pick some point $x_0$. Choose a search direction $s_1$ parallel to the axes that maximizes the gradient slope. Locate the minimum of $f$ on the line through $x_0$ parallel to $s_1$ to get $x_1$. Comb the remaining basis vectors of $\mathbb R^n$ to be orthogonal to $s_1$ with respect to the inner product defined by $C$. Next find a search direction $s_2$ among these combed vectors that maximizes the gradient slope, locate the minimum of $f$ on the line through $x_1$ parallel to $s_2$ to get $x_3$, etc. After at most $n$ steps the algorithm will terminate.

I went ahead and programmed this into C# using the MathNet.Numerics.LinearAlgebra library and it seems to work quite well.

First we have a helper function that updates the matrix $GS$, followed by the main minimization routine.

public Matrix<double> GramSchmidt(Matrix<double> matrix, Vector<double> d, int k, bool[] J)

        {

            Matrix<double> A = matrix.Clone();

            double lambda = d.DotProduct(matrix.Column(k));

            A.SetColumn(k, matrix.Column(k) / lambda);

            for(int i = 0; i < A.ColumnCount; i++)

            {

                if (!J[i]&&i!=k)

                {

                    double mu = d.DotProduct(matrix.Column(i));

                    A.SetColumn(i, matrix.Column(i)-matrix.Column(k) * (mu / lambda));

                }

            }

            return A;

        }

 public Vector<double> Minimize(Vector<double> c,Matrix<double> C,Vector<double> x0, out bool unbounded)

        {

            unbounded = false;

            Vector<double> x = x0;

            int n = c.Count;

            bool[] J = new bool[n];//builds array of falses

            Matrix<double> GS = Matrix<double>.Build.DenseIdentity(n);

            Vector<double> g = c + C * x;

            for(int j = 0; j < n; j++)

            {

                //Step 1:

                bool init = true;

                foreach(bool index in J) { init = init && index; }

                if (init) { return x; }

                double max = 0;

                int k = -1;

                for(int i =0; i < n; i++)

                {

                    if (!J[i] && Math.Abs(g.DotProduct(GS.Column(i))) > max)

                          { max = Math.Abs(g.DotProduct(GS.Column(i)));k = i; }

                }

                if (g.DotProduct(GS.Column(k)) == 0) { return x; }

                Vector<double> s = GS.Column(k);

                if (g.DotProduct(GS.Column(k)) < 0) { s *= -1; }

                //Step 2:

                double normSquaredOfs = (s.ToRowMatrix() * C * s.ToColumnMatrix())[0,0];

                if (Math.Abs(normSquaredOfs) < .0000001) { unbounded = true; return x; }

                double t = g.DotProduct(s) / normSquaredOfs;

                //Step 3:

                x -= t * s;

                g = c + C * x;

                Vector<double> d = 1/ Math.Pow(normSquaredOfs, 0.5)*(C * s);

                GS = GramSchmidt(GS, d, k, J);

                J[k] = true;

            }

            return x;

        }

Finding an orthonormal basis with respect to a symmetric positive definite form

Let $C$ be an $n\times n$ symmetric positive definite matrix. Define an inner product on $\mathbb R^n$ by the formula $\langle x,y\rangle =x^TCy$.

 In the last post we wanted to find a set of conjugate directions $s_1,\ldots,s_n$ which were a set of vectors in $\mathbb R^n$ such that $\langle s_i,s_j\rangle=\delta_{ij}$, where $\delta$ is the Kronecker delta function. To do this we start with any basis of $\mathbb R^n$ $b_1,\ldots,b_n$ and slowly improve it. 

At the $k$th stage of the process we will have a list of vectors $s_1^{(k)},\ldots s_n^{(k)}$ where the first $k$ vectors span $S_k$ and the last $n-k$ vectors span $W_k$. The vectors $s_1^{(k)},\ldots s_k^{(k)}$ form an orthonormal basis for $S_k$ and $\mathbb R^n = S_k\oplus W_k$ is an orthogonal decomposition. 

The basic process (called the Modified Gram-Schmidt process) is as follows. Let $v_1,\ldots,v_k$ be a set of independent vectors.

• Let $\hat{v}_1 = v_1/\sqrt{\langle v_1,v_1\rangle}$.
• Let $\hat{v}_i = v_i-\langle v_i,\hat{v}_1\rangle\hat{v}_1$ for all $i\geq 2$.

Then $\hat{v}_1$ is orthogonal to the subspace generated by $\hat{v}_2,\ldots, \hat{v}_k$.

$$\langle \hat{v}_i,\hat{v}_1\rangle = \langle v_i,\hat{v}_1\rangle - \langle v_i,\hat{v}_1\rangle\langle \hat{v}_1,\hat{v}_1\rangle$$

$$\hspace{4em}=\langle v_i,\hat{v}_1\rangle - \langle v_i,\hat{v}_1\rangle\cdot 1=0.$$

This process iterates as follows. Suppose we have a decomposition $S_k\oplus W_k$. Then use this procedure to find a basis for $W_k$ of the form $c_{k+1},\ldots,c_n$ where $c_{k+1}$ has norm $1$ and is orthogonal to the rest of the vectors. Then $c_{k+1}$ is orthogonal to $S_k$ as it lies in $W_k$, so it can be added to $S_k$ to get a new decomposition $S_{k+1}\oplus W_{k+1}$.

In the next post we will use this algorithm to solve the unconstrained quadratic minimization problem. There are a couple things to note. The conjugate directions that are produced this way have a choice of sign, and we will need to pick the right sign when implementing the algorithm. The other thing to note it that we started from the first vector in the list and then steadily moved from left to right in creating the orthonormal basis. However, we have a choice of what order to do this in, and the algorithm will be set up in such a way as to pick certain directions as making the most progress. As far as I can tell this is not important because the unconstrained algorithm will always take $n$ steps to reach a solution (unless you happen to get lucky in your choice of initial point). My guess is that keeping things flexible will help to generalize the argument later.

Quadratic Programming: The unconstrained case

In order to build our understanding of how quadratic programming works, we are first going to look at the unconstrained case:

• Minimize $f(x) = c^Tx+\frac{1}{2} x^TCx$ where $C$ is a positive definite matrix.

Note this is a global optimization problem and can be solved using calculus. The minimum occurs where the gradient vanishes:

$$\nabla f = c^T+Cx = 0.$$

So one could simply solve this linear system in whatever way one wishes and be done with the problem, but following Best's textbook, we will present a basic algorithm that we can build off of in the constrained cases.

The basic idea is as follows. Start with any initial approximation $x_0$, and suppose we have a search direction in the form of a vector $s_0$. We want to walk along this search direction as far as possible to minimize $f(x)$.

Let $g_0=\nabla f(x_0)$. Suppose $g_0^Ts_0\neq 0$. Replace $s_0$ by $-s_0$ if necessary so that $g_0^Ts_0>0$. Looking along the line $x_0-t s_0$, we write the Taylor Series as

$$f(x_0-ts_0) = f(x_0) -tg_0^Ts_0+\frac{1}{2}t^2s_0^TCs_0$$

Note that $s_0^TCs_0>0$ since $C$ is a positive definite matrix. So this function has a minimum, which can be found by taking the derivative with respect to $t$:

$$\frac{d}{dt}f(x_0-ts_0)=-g_0^Ts_0+ts_0^TCs_0=0$$

$$t=\frac{-g_0^Ts_0}{s_0^TCs_0}=\frac{-\nabla f(x_0)^Ts_0}{s_0^TCs_0}.$$

This quantity is called the optimal step size and tells you how far to walk along the ray to minimize the objective function on that ray.

One can iterate this by choosing a new search direction, so the question becomes how to choose these directions. A convenient choice of search directions is given by the concept of conjugate directions.

Definition: A set of vectors $s_1,\ldots,s_k$ is said to be conjugate (with respect to $C$) if $s_i^TCs_j=0$ for $i\neq j$ (and $s_i^TCs_i>0$). 

Theorem: Let $s_0,\ldots,s_{n-1}$ be conjugate directions, and $x_0$ be arbitrary. Construct $x_1,\ldots,x_n$ where $x_i=x_{i-1}-t_{i-1}s_{i-1}$, with $t_{i-1}=\frac{-\nabla f(x_{i-1})^Ts_0}{s_{i-1}^TCs_{i-1}}$. Then

• $\nabla f(x_{j+1})s_i=0$ for $0\leq i\leq j\leq n-1$
• $x_n$ is the optimal solution for the QP.

In other words, at each iteration the gradient becomes orthogonal to more search directions. In the end it is orthogonal to all of them, and since they span $\mathbb R^n$ (exercise) the gradient must be zero.

Proof:

First we want to show that $\nabla f(x_{j+1})^T s_j = 0$. To see this note that

$$\nabla f(x_{j+1})=\nabla f(x_j) -t_jCs_j$$ by Taylor series, and taking the inner product with $s_j$ we get

$$\nabla f(x_{j+1})^Ts_j=\nabla f(x_j)^Ts_j -t_js_j^TCs_j=0$$ by definition of $t_j$.

Now choose $i<j\leq n-1$. Again by Taylor series, we have

$$\nabla f(x_{j+1})=\nabla f(x_{i+1})+C(x_{j+1}-x_{i+1}).$$

Taking the inner product with $s_i$:

$$\nabla f(x_{j+1})^Ts_i=\nabla f(x_{i+1})^Ts_i+(x_{j+1}-x_{i+1})^TCs_i.$$

$$\hspace 4em = (x_{j+1}-x_{i+1})^TCs_i.$$

$$\hspace 4em = -t_{i+1}s^T_{i+1}Cs_i-\cdots-t_js^T_jCs_i.$$

Each term is equal to $0$ by definition of conjugate directions. $\Box$

The only missing ingredient is identifying a set of conjugate directions for a matrix $C$. We will do this in the next post.

Quadratic Programming

 In my next series of posts, I would like to start exploring the quadratic programming problem. This is similar to the linear case, except the objective function is allowed to be a convex quadratic function 

$$f(x) = c\cdot x+\frac{1}{2} x^TCx$$

where $C$ is a positive semidefinite matrix. We are still considering the constraints to be linear, and in fact, following Best's textbook, we will assume they are inequality constraints of the form $Ax\leq b$. 

Definition: Let $x_0\in\mathbb R^n$. We say a constraint $a_i\cdot x\leq b_i$ is active at $x_0$, if $a_i\cdot x_0=b_i$.

In other words, $x_0$ lies on the facet of the feasible region defined by the $i$th constraint. 

Definition: Let $I(x_0)$ be the set of indices of active constraints for $x_0$. We say $x_0$ is quasi-stationary if $x_0$ is the optimal solution for the problem

• Minimize $c\cdot x+\frac{1}{2} x^TCx$
• Subject to $a_i \cdot x=b_i, i\in I(x_0)$

A couple of things to note: any corner of the feasible region is quasi-stationary since the active region is just a point. 

Lemma: An optimal solution for a QP must be quasi-stationary.

Proof: Let $I$ be the set of active constraints for an optimal point $x_0$. Note that $I$ could be empty.

We need to show that $x_0$ is optimal on the entire facet defined by its active constraints. This is obvious, due to convexity of the objective function. If there were a point with smaller objective function, just walk toward it until you hit a facet of higher codimension, meaning that you had not included all active constraints. $\Box$

So we have a very similar situation to the linear programming case. Instead of looking at basic variables, we look at active constraints. Just as in the LP case, a crude but inefficient algorithm is available. Just look at all possible subsets of the constraints, solve them as equalities, and then find the global minimum of the objective function on the hyperplane defined by that facet.

We will conclude this first post on quadratic programming by deriving tests for a point to be quasi-stationary or optimal. 

Suppose we have a quasi-stationary point $x_0$. Recall $\nabla f$ is the vector which gives the direction of fastest increase of $f$. Because $x_0$ is optimal for its set of active constraints, $\nabla f$ must be orthogonal to the hyperplane given by those constraints. Hence $\nabla f$ is in the linear span of $\{a_i\,|\,i\in I(x_0)\}$:

$$-\nabla f = \sum_{i\in I(x_0)} u_i a_i,$$

where we introduce the negative sign for later convenience. If $x_0$ is actually optimal, then the direction of fastest increase must point into the feasible polytope, which means that $u_i\geq 0$ in the above.

A convenient way to encode this is as follows:

Theorem: A point $x_0$ is quasi-stationary for the QP 

$$\min\{f(x)\,|\, a_i\cdot x\leq b_i\}$$ if and only if

• $x_0$ lies in the feasible region $R$: $\forall i\,\,a_i\cdot x_0\leq b_i $
• $-\nabla f(x_0) = \sum_{i=1}^m u_ia_i$ for some scalars $u_i\in\mathbb R$
• $\forall i\,\, u_i(a_i\cdot x_0-b_i)=0$

The last condition gives an alternative: either $u_i=0$ or $i\in I(x_0),$ so is just a reformulation of the statement that the gradient is a linear combination of the active constraint coefficients.

Theorem: A point $x_0$ is optimal if in addition to the above three conditions, each $u_i$ is nonnegative. Writing this in matrix notation:

• $Ax_0\leq b$
• $-\nabla f(x_0) = A^Tu$ for some $u\geq 0$
• $u^T(Ax_0-b)=0$

Linear programming versus row echelon form

 As mentioned in my previous post, my plan was to illustrate the LP algorithm we have developed so far over the last several posts with a real world example. I actually did take an example from a project I am working on with a coefficient matrix $A$ of dimensions $5\times 14$. Obviously this must be a massively overdetermined system, but at the moment, my program simply feeds it through the steps by first constructing a phase 1 and then an auxiliary LP.

This did not seem like a great example to show because several rows were simply zero and other rows were copies or negations. So I decided to row reduce it first. However, one thing I noticed Wikipedia ominously mention, is that Gaussian elimination is numerically unstable, and indeed when I did naive elimination, the resulting system was inconsistent. I was able to get results if I used the CoerceZero() method first, which coerces any elements of a matrix sufficiently close to 0 to actually be 0. However, I also noticed a couple elements of the row echelon form were absurdly big. 

So I actually think that row reducing the matrix is actually a bad idea due to the numerical instability. It's possible a more sophisticated approach using the $LU$ decomposition might be better, but in fact, the program works great at the moment without doing the reduction, so perhaps its okay to leave well enough alone, although I will probably at least remove the all $0$ rows. 

I would also like to note that the algorithm for computing row echelon form which I found at Rosetta Code does not appear to be correct, to add insult to injury.

Linear programming - the auxiliary LP

 Yesterday, I wrote a series of posts explaining how to solve a linear programming problem of the form

• Minimize the objective function $f(x)=c\cdot x$
• subject to the constraints
• $Ax=b$ where $A$ is an $m\times n$ matrix
• $x\geq 0$

The simplex algorithm works great, but you need an initial basic feasible solution to get started. Luckily, one can construct a Phase I program to decide if any feasible solutions exist, and provide such a solution if so. The last remaining piece is to turn this into a basic feasible solution, which we can do with the help of an auxiliary LP.

Recall that the phase I LP is the following

Phase 1 LP (standard)

Introduce new variables $y=(y_1,\ldots,y_n)$.

• Minimize $\sum_{i=1}^n y_i=e\cdot y $ where $e=(1,\ldots,1)$
• subject to
• $Ax+Iy=b$
• $x,y\geq 0$

If we have an optimal solution $(x^*,y^*)$ with $y^*=0$, then $x^*$ is a solution to the original problem. If all of the $y^*$ variables are non-basic, then we can easily construct a bfs for the original LP. However it is quite possible that this is not the case. To remedy this, we introduce the auxiliary LP.

Auxiliary LP

Start with the phase 1 LP and add one new variable $z$.

• Minimize $c\cdot x$
• subject to
• $Ax+Iy=b$
• $e\cdot y+z = 0$ 
• $x,y,z\geq 0$

The auxiliary LP is useful for the following reasons:

• Fact 1: If $(x,y,z)$ is optimal for the auxiliary LP, then $x$ is optimal for the original LP.
• Fact 2: If the auxiliary LP is unbounded, then the original LP is unbounded.
• Fact 3: If $(x^*,0)$ is an optimal solution to the phase 1 LP, then $(x^*,0,0)$ is a bfs for the auxiliary LP. The basic variables are the same as those for $(x^*,0)$ with $z$ added.

Facts 1 and 2 basically come from the fact that any solution to the auxiliary LP must have $y=0$ and $z=0$, because the sum of nonnegative terms $e\cdot y+z$  is equal to $0$. 

The main thing about Fact 3 is to check that the new matrix, after appending the column $\left(\begin{matrix} 0\\ 1\end{matrix}\right)$, is invertible, which I leave as an exercise.

So now we have a complete algorithm to solve LPs of the form indicated at the start of the post. I've implemented this algorithm in C# and it works quite well. The one issue that I kept grappling with is the issue of thresholds, and I essentially just kept tinkering with them until the program stopped giving me error messages. The linear systems I have been dealing with are often highly degenerate, and one thing I am picking up from these posts is that I may want to row reduce the coefficient matrix and throw away extraneous rows before implementing the algorithm. 

So now the question is, what's next? I would like to understand how to do quadratic programming, but I would also like to run through at least one real life (but smallish) example of the LP algorithm we just described both for illustrative purposes, but also so I can see how I might be able to improve it.

Linear Programming - Phase 1 algorithms

 In my previous few posts I discussed solving the LP

• Minimize $c\cdot x$
• Subject to
• $Ax=b$ and 
• $x\geq 0.$

Here $A$ is an $m\times n$ matrix, while $c$ and $x$ are $n$-vectors.

We went through the basic simplex algorithm, but it required an initial basic feasible solution to get started. In this post I want to discuss two algorithms for finding a feasible solution.

We begin with what appears to be the standard Phase 1 LP, taken from Solow's textbook.

Phase 1 LP (standard)

Introduce new variables $y=(y_1,\ldots,y_n)$.

• Minimize $\sum_{i=1}^n y_i=(1,\ldots,1)\cdot y $
• subject to
• $Ax+Iy=b$
• $x,y\geq 0$
• Initial bfs: $B=I$ is formed from the last $m$ columns of  the coefficient matrix $[A \, I]$.

Theorem: Let $(x^*,y^*)$ be an optimal solution to the phase 1 LP. Then the original LP is feasible iff $y^*=0$. In that case, $x^*$ is a solution.

Proof:

Clearly if $y^*=0$, then $b=Ax^*+Iy^*=Ax^*$, so $x^*$ is feasible for the original LP. Moreover, $x^*\geq 0$ since $(x^*,y^*)$ is feasible for the phase 1 LP.

Now suppose that the original LP is feasible with solution $x^*$. Then $(x^*,0)$ is feasible for the phase 1 LP. Moreover the objective function $f(x,y)=\sum y_i$ is zero on this solution. Moreover $f(x,y)>0$ if $y\neq 0$, so that means any optimal solution must have $y^*=0$. $\Box$

This gives us a very pretty solution to the question of when the feasible region is nonempty, but if we are to use it to continue with the simplex algorithm, we still need to have a bfs for the original LP, not just a solution. We will return to this question later when we talk about auxilliary LPs. For now, let us look at another Phase 1 LP from Best's book on Quadratic Programming.

The set up for this algorithm is slightly more general. We are looking at a region expressed by the constraints

• $a_i^Tx\leq b_i$ for $i=1,\ldots, m$
• $a_i^Tx = b_i\geq 0$ for $i=m+1,\ldots,m+r$

The question is whether the region carved out by these constraints is nonempty.

Phase 1 LP (Best's Algorithm)

Let $d=-\sum_{i=m+1}^{m+r}a_i$, and introduce a single additional variable $\alpha$. 

• Minimize $d\cdot x+\alpha$.
• Subject to:
• $a_i^Tx-\alpha \leq b_i$ for $i=1,\ldots m$
• $a_i^Tx \leq b_i$ for $i=m+1,\ldots,m+r$
• $-\alpha\leq 0$

Best says that taking $x=0$ and $\alpha$ sufficiently large is a  feasible solution to the phase I LP, which is true, but he does not say that it is a basic feasible solution. Indeed, this is an inequality constrained LP, so the discussion from our previous posts does not apply. In his book, Best develops an algorithm for optimizing a quadratic objective function subject to both equality and inequality constraints, which he specializes to the case of an linear objective function. I will have to dig into that algorithm more to see what initial data is required. Perhaps simply a feasible solution is enough.

 

Leaving these concerns aside, let us check to see if Best's algorithm really does determine if the original LP is feasible.

Proposition: Best's Phase 1 LP is bounded from below. Hence an optimal solution exists. 

Proof:

The objective function 

$$d\cdot x+\alpha=\alpha-\sum_{i=m+1}^{m+r}a_i\geq \alpha -\sum_{i=m+1}^{m+r}b_i\geq-\sum_{i=m+1}^{m+r}b_i.$$

$\Box$

Theorem: Let $(x^*,\alpha^*)$ be an optimal solution to Best's Phase 1 LP. Then the initial LP is feasible iff $\alpha^*=0$.

Proof: We start with the lower bound on the objective function 

$$d\cdot x^*+\alpha^*\geq -\sum_{i=m+1}^{m+r}b_i.$$

If the initial LP were feasible with solution $x^*$, then $(x^*,0)$ would make the above inequality an equality. Hence if $d\cdot x^*+\alpha^*> -\sum_{i=m+1}^{m+r}b_i$, the initial LP is not feasible. If on the other hand $d\cdot x^*+\alpha^*= -\sum_{i=m+1}^{m+r}b_i$, then $\alpha=0$ and the inequalities $A_ix\leq b$ for $i=m+1,\ldots,m+r$ must actually be equalities, implying $x^*$ is a feasible solution to the original LP. $\Box$

So it looks like Best's Phase 1 algorithm checks out, and is actually quite clever.

Linear Programming: pivoting

 Let us summarize our progress so far. 

• $x=(x_B,x_N)=(B^{-1}b,0)$ is a basic feasible solution to $Ax=b, x\geq 0$, 
• $j^*$ is an index such that $d:=(c_N-c_BB^{-1}N)_{j^*}<0$ and
• $t^*=\min\{-x_i/d_i\,|\, 1\leq i\leq n\text{ and } d_i<0\}$. Let $k^*$ be a choice of index where this minimum is realized.

Then we have concluded that $x+t^*d$ is a feasible solution. In this post we want to show that it is a bfs, and how to calculate the new index sets.

One thing to note: if $t^*=0$ then we aren't actually moving to a different point, but we are moving to a different bfs representation of that same point. In this case the objective function doesn't strictly decrease, which means that if we repeat the basic steps in our algorithm, we might cycle back to the same point. Apparently cycling is quite rare in practice, though one can also preclude its happening through the right choice of pivoting rule - i.e. which $j^*$ and $k^*$ to pick! (One such choice is Bland's Rule.) 

Theorem: With the hypotheses above, $x+t^*d$ is a bfs $(x_{B'},x_N')$ where $B'$ is formed by replacing column $k^*$ of $B$ with column $j^*$ of $N$. 

Proof: By construction, we zeroed out the $x_{k^*}$ variable, so $x_{N'}=0$. If we can show $B'$ is invertible, then automatically $x_{B'}=(B')^{-1}b$ and we are done. To see that $B'$ is invertible, we claim that $B'=BE$ where $E$ is formed by replacing column $k^*$ of the $m\times m$ identity matrix with $-d_B$. 

If $k\neq k*$, $(BE)_{\cdot k}=BE_{\cdot k}=BI_{\cdot k}= B_{\cdot k}$.

On the other hand $(BE)_{\cdot k^*}=BE_{\cdot k^*}=B(-d_B)=B(B^{-1}N_{\cdot j^*})=N_{\cdot j^*}.$ By definition the column vector $N_{\cdot j^*}$ is equal to $B'_{\cdot k^*}$.

Finally we need to show that $E$ is nonsingular. This follows because $(d_B)_{k^*}\neq 0$. So one can use row operations to clear out the rest of the $k^*$ column to get a diagonal matrix with $1$'s on the diagonal, except for an occurrence of $(d_B)_{k^*}$. $\Box$

This operation of updating $B$ to $B'$ is called pivoting.

Now we have all the ingredients for the basic simplex algorithm to minimize an objective function $c\cdot x$ subject to the constraints $Ax=b$ and $x\geq 0$.

1. Start with an initial bfs $x=(x_B,x_N)=(B^{-1}b,0)\geq 0$.
2. Compute $c_N-c_BB^{-1}N$.
3. If  $c_N-c_BB^{-1}N\geq 0$, $x$ is optimal. Terminate program.
4. Otherwise, select $1\leq j^*\leq n-m$ such that $(c_N-c_BB^{-1}N)_j<0.$
5. Compute $d_B=-B^{-1}N_{\cdot j^*}.$
6. If $d_B\geq 0$, the LP is unbounded. Terminate program.
7. Otherwise, select $1\leq k^*\leq m$ so that $-x_{k^*}/d_{k^*} = \min\{-x_i/d_i\,|\, 1\leq i\leq n\text{ and } d_i<0\}$.
8. Create a new bfs by replacing the $k^*$ column of $B$ with the $j^*$ column of $N$. 
9. Go to step 2.

This is all very well, but how do we find that initial bfs? How do we even determine if there is a feasible solution at all? That was actually my primary motivation for looking into this. How can we find a solution to $Ax=b$ such that $x\geq 0$? Let's not even worry about the objective function!

In the next post, I want to introduce what is called the phase I program. This is a linear program formed from the original, which has the property that it has an obvious bfs, and such that the solution can tell us if the initial problem was feasible. In fact, I plan to explore two different constructions. One is given in the book Linear Programming by Daniel Solow. The other is from the book Quadratic Programming with Computer Programs by Michael J. Best. Solow's phase I algorithm looks pretty standard from the various sources I have looked at. Best's algorithm appears to me much more efficient, requiring only one additional variable, and if it works, it seems to be a major improvement. I am a bit skeptical, but hopefully we can sort it out in subsequent posts.

Linear programming: improving a basic feasible solution

 Recall our basic set up. We have a linear program

• Minimize $c\cdot x$ where $x$ and $c$ are $n$-vectors.
• subject to the constraints
• $Ax=b$ ($A$ is an $m\times n$ matrix.
• $x\geq 0$

and we have decided to look for basic feasible solutions, which correspond to partitioning the index set of columns of $A$ into the $B$ and $N$ indices. Somewhat abusing notation, $B$ and $N$ are also the submatrices of $A$ formed by these columns. $x$ is said to be a basic feasible solution if $x_N=0$, and $x_B=B^{-1}b\geq 0$.

In the last post we saw that if $x$ passes the test for optimality $c_N-c_BB^{-1}N\geq 0$, then it is optimal. In this post we consider what happens if our bfs fails the test. Our first step will be to find a direction $d$ along which the objective function decreases. For calculus aficionados, denoting the objective function by $f(x)=c\cdot x$, the directional derivative in the direction of $d$, is just given by $c\cdot d$. Hence we want to find a direction $d$ such that $c\cdot d<0$. (We also need the direction to point into the polytope.)

Theorem: If $x$ is a bfs which fails the test for optimality, let $j$ be an $N$-index  such that $1\leq j\leq n-m$ and $(c_N-c_BB^{-1}N)_j<0$. Let $d=(d_B,d_N)=(-B^{-1}N_{\cdot j},I_{\cdot j})$, where $I$ is the $(n-m)\times(n-m)$ identity matrix. Then $c\cdot d<0$.

Proof:

$$c\cdot d = c_Bd_B+c_Nd_N=-c_BB^{-1}N_{\cdot j}+c_NI_{\cdot j}$$

$$\hspace 4em = (c_N-c_BB^{-1}N)_j<0$$

$\Box$

So we have a direction $d$, possibly several. Let $j^*$ be some choice of index such that $(c_N-c_BB^{-1}N)_j<0$ which gives us this direction. Next we must determine how far to go in the direction $d.$ That is, if we have a bfs $x$, and a direction $d$, how big can we make $t$ so that $x+td$ remains in our feasible polytope. One nice thing is that the equality constraints $Ax=b$ remain true along the entire line if $d$ is of the special form $(-B^{-1}N_{\cdot j},I_{\cdot j})$.

Lemma: Suppose $x=(x_B,x_N)=(B^{-1}b,0)$ is a bfs for $Ax=b$ and $d=(d_B,d_N)=(-B^{-1}N_{\cdot j},I_{\cdot j})$ for some $j$, then for all $t$ $A(x+td)=b$.

Proof:

It suffices to show that $Ad=0$. 

$$Ad=Bd_B+Nd_N=B(-B^{-1}N_{\cdot j})+NI_{\cdot j}=N_{\cdot j}-N_{\cdot j}=0.$$

$\Box$

So we need to find the maximal $t$ so that the nonnegativity constraints $x\geq 0$ are satisfied.

Lemma: Suppose $x\geq 0$ is a solution for $Ax=b$ and $d$ is a direction having at least one negative component, then

$$t^* = \min\{-x_i/d_i\,|\, 1\leq i\leq n\text{ and } d_i<0\}$$ has the property that $x+t^*d\geq 0$. In particular, if $x,d$ are as in the previous lemma, then $x+t^*d$ is feasible.

If $d\geq 0$, the LP is unbounded: the entire ray $x+td, t\geq 0$ lies in the feasible region and $f$ is decreasing along it.

Proof of lemma: 

By the previous lemma, we just need to satisfy the inequality constraints $x\geq 0$. Let us check for each coordinate $x_i$. If $d_i\geq 0$, then $(x+td)_i = x_i+td_i\geq 0$ for all $t>0$. Otherwise suppose $d_i<0$. Then $(x+t^*d)_i = x_i+t^*d_i$, and we know $t^*\leq -x_i/d_i$ by definition, so $t^*d_i\geq -x_i$ and $x_i+t^*d_i\geq x_i-x_i=0$.

$\Box$

In summary, we have chosen a direction along which the objective function decreases that lies in the affine subspace $Ax=b$. Then we have calculated the maximum distance one can travel along this ray and stay in the positive orthant. 

In the next post we will reinterpret this operation combinatorially and show that $x+t^*d$ is also a bfs.

Linear programming and the test for optimality.

[Note: I am using Linear Programming by Daniel Solow as a reference for this and future posts. This is an economical Dover Book which I highly recommend!]

Our goal in this post is to describe an algorithm for solving a particular type of linear program, one which is an optimization problem of the following form:

• Minimize $f(x)=c\cdot x$ 
• Subject to
• $Ax=b$
• $x\geq 0$

Let's just take a step back and think about this problem for a second. $Ax=b$ is some affine subspace of $\mathbb R^n$, and we are considering its intersection with the positive orthant. This gives us some higher dimensional polytope, and we are trying to identify the point that minimizes $f$. Since $f$ is a linear function, such a point, if it exists, will have to be one of the vertices of this polytope. Unfortunately, a polytope can have exponentially many vertices in the number of constraints, so any method that seeks to enumerate all vertices first is going to be very slow.

I am going to outline a method, due to Dantzig, called the simplex method. There are two ingredients of this method that make it work:

• There is a simple test to check if a point is optimal. So if you happen to land on an optimal point, you can easily check if you are done.
• If you are not on an optimal point, you can find a direction where $f$ decreases, and walk along an edge of the polytope to find a new point. 

You do need an initial vertex to get the algorithm started (more on that later), but once you have it, you can keep walking along edges, decreasing $f$ each time, until you find the minimum.

Next lets think about what a vertex of this polytope is combinatorially. To begin, we will assume that the solution set to $Ax=b$ is generic in the sense that the rows of $A$ are linearly independent. That means it carves out a codimension $m$ affine subspace of $\mathbb R^n$. Generically, this will intersect an $m$ dimensional subspace in a point, so the corners of our polytope will come from setting $n-m$ coordinates equal to $0$. And indeed, if we set $n-m$ variable equal to $0$, we are left with the same number of unknowns as equations, and we expect there to be a unique solution generically. This motivates the definition of what is known as a "basic feasible solution" or bfs for short.

Definition: A basic feasible solution (bfs) to the linear system $Ax=b, x\geq 0$, is a vector $x$ such that

• there is a nonsingular $m\times m$ submatrix of $A$ formed by choosing $m$ columns, denoted by $B$, and the remaining columns are denoted by $N$ such that
• $x_B=B^{-1}b$
• $x_B\geq 0$, and
• $x_N=0$

Here $x_B$ and $x_N$ are the projections of $x$ onto the subspace spanned by the $B$ and $N$ columns of $A$ respectively. 

Again, a bfs is generically nothing more than a vertex of the solution polytope in the first orthant. If $Ax=b$ is not generic, one could first row reduce the system and throw out the $0$ rows. There are other options available as well to create a system with an initial bfs.

Let us return now to the two ingredients I mentioned above, a test for optimality, and a method to determine a direction to walk if the current point is not optimal.

Test for optimality: Consider the linear program to minimize $f(x)=c\cdot x$ subject to $Ax=b, x\geq 0$. Let $x^*$ be a bfs with $(x_B,x_N)=(B^{-1}b,0)\geq 0$. If $$c_N-c_BB^{-1}N\geq0$$ then $x^*$ is optimal for this LP.

Proof:

Let $x=(x_B,x_N)$ be some other feasible solution. We want to show $c\cdot x\geq c\cdot x^*$.

$$ c\cdot(x-x^*) = c_B(x_B-x_B^*)+c_N(x_N-x_N^*) $$

$$\hspace{2em}=c_B[B^{-1}(b-Nx_N)-B^{-1}b]+c_N(x_N-0)$$

Here we use that $x$ is feasible, so $Ax=Bx_B+Nx_N=b$, and solving for $x_B$, we get $x_B=B^{-1}(b-Nx_N)$.

Continuing,

$$c_B[B^{-1}(b-Nx_N)-B^{-1}b]+c_N(x_N-0)=c_B[B^{-1}b-B^{-1}Nx_N-B^{-1}b]+c_Nx_N$$

$$\hspace{4em}=-c_BB^{-1}Nx_N+c_Nx_N$$

$$\hspace{4em}=(c_N-c_BB^{-1}N)x_N.$$

So to reiterate, for any feasible solution $x$, and any bfs $x^*$, we have the equation

$$c\cdot(x-x^*)=(c_N-c_BB^{-1}N)x_N.$$

By hypothesis, $(c_N-c_BB^{-1}N)\geq 0$, and by feasibility $x_N\geq 0$. So if $x^*$ passes the optimality test, then $c\cdot x\geq c\cdot x^*$ for any feasible solution. Hence it really is optimal. $\Box$

Now, what about the converse? Can we show that if a feasible solution is optimal, then it is a basic feasible solution satisfying the optimality test? This is a good question which we will return to later. To get some intuition, let's do an example.

Example: Minimize $x+2y+3z,$ subject to $x+y+z=1$ and $x,y,z\geq0$. 

There are $3$ basic feasible solutions. $(1,0,0)$, $(0,1,0)$ and $(0,0,1)$. In each case $B$ is the $1\times 1$ matrix with a $1$ in it. So $B=B^{-1}=(1)$, while $N=(1,1)$. For the first basic feasible solution we have $c_N=(2,3)$ and $c_B=(1)$. So we get

$$c_N-c_BB^{-1}N=(2,3)-(1)*1*(1,1)=(1,2)\geq 0.$$

So indeed $y=z=0$ gives an optimal solution. It is interesting to note that this works even in the degenerate case when we consider $x=y=z=0$. 

This seems like a good place to stop for this post. In the next post, we'll look at how to find an improved bfs if your existing bfs fails the test for optimality.

Solving linear systems, the pseudoinverse, and linear programming.

Suppose we want to solve a system of linear equations $Ax=b$. There are many ways to go about this. If $A$ happens to be a square invertible matrix, you can just write $x=A^{-1}b$. If $A$ is not square, or is not invertible, then we have to look for other methods. One trick that often works is to multiply both sides of the equation $Ax=b$ by the transpose of $A$, to get

$$A^TAx=A^Tb$$

The matrix $A^TA$ is a square matrix, and in good cases, it is actually invertible, so we have

$$x=(A^TA)^{-1}A^Tb.$$

This is actually a special case of the Moore-Penrose pseudoinverse of a matrix $A^+$. When $A^TA$ is invertible, we have that $A^+=(A^TA)^{-1}A^T$, but $A^+$ is defined for any matrix.

It is not my intention to get deep into the theory of pseudoinverses, but I do want to point out some important properties. (See wikipedia for more detail.)

• The system $Ax=b$ has at least one solution if and only if $AA^+b=b$.
• If a solution exists, $x=A^+b$ is the solution with smallest Euclidean norm. (i.e. it is closest to the origin.)
• If a solution exists, one can obtain all solutions by $x=A^+b+(I-A^+A)w$ for an arbitrary vector $w$.  

Calculation of pseudoinverses is widely supported in software libraries. For example, in C#, I use the MathNet numerics package, and it is a convenient way to solve linear systems.

Instead of finding the solution closest to the origin, what if you want to find the solution closest to some other point $x_0$?

Exercise: Let $A$ be an $m\times n$ matrix. What is a method for finding the solution to a system of equations $Ax=b$ which is closest to a given $x_0\in\mathbb R^n$.

Solution: Subtract $Ax_0$ from both sides of the equation to get: $$A(x-x_0)=b-Ax_0.$$

Then $x-x_0=A^+(b-Ax_0)$ is the solution such that $x-x_0$ is closest to $0$. Thus $$x=A^+(b-Ax_0)+x_0$$ is closest to $x_0$. $\Box$

We have just solved a simple quadratic optimization problem, that of minimizing the distance to a point of a solution to a linear system. 

Let us consider a slightly more complex question. 

Problem: Minimize the distance to a point $x_0$ of a solution to $Ax=b$, assuming each coordinate of $x$ is nonnegative.

This is a more difficult question, and even the simpler feasibility question is not so obvious:

Feasibility Problem: Does there exist a solution to $Ax=b$ where each coordinate of $x$ is nonnegative?

The feasibility problem can be solved using the technique of linear programming, which for our purposes means that we want to minimize some linear objective function $f\colon \mathbb R^n\to\mathbb R$ subject to the constraints $Ax=b$ and $x\geq 0$. (Although it looks like a typo, the convention is that $x\geq0$ means each coordinate of $x$ is nonnegative.) At first glance, this seems like a harder problem than determining feasibility, but bear with me. Once we have set up the algorithm for solving linear programs of the above form, we will be able to easily apply it to the feasibility issue.

In order to make these posts as readable as possible, I'll stop here and continue in the next post to talk about linear programs.

Toward the 4D Poincare Conjecture

The Disc Embedding Theorem is a recent book carefully explicating the work of Mike Freedman in proving the topological 4D Poincare Conjecture. I am very excited about this book. I have always wanted to understand the proof of the theorem, but prior to this book, there has not been a feasible accessible approach to it. A dream of mine is to be able to actually see one of the crazy topological disks that the proof creates, for example bounding a topologically slice knot.

The book is based on a series of lectures that Freedman gave, though saying that tends to obscure the amount of work the contributors of the book put into it! One nice feature of the book is its reflection of the structure of Freedman's original presentation, introducing important techniques through historical examples.

One such example is the "shrinking" proof of the Schönflies Theorem. Not only is this a cool and foundational result in topology, the argument introduces the concept of shrinking, which figures into the 4D proof.

I'd actually like to go through this proof of Schönflies. It is super elegant and really showcases the power of point set topology as found in Munkres's standard text that so many of us learned from as undergraduates!

Definition: A subset $X$ of an $n$-manifold $M$ is said to be cellular if it is a nested intersection of countably many closed $n$-cells. More precisely, $X=\cap_{n\geq 1} C_n,$ where $C_n$ is homeomorphic to the closed $n$ dimensional ball $D^n$, and $C_{n+1}\subset \operatorname{int}C_n.$

I invite the reader to come up with examples of cellular sets. The one given in the book is a literal letter $X$, which you can imagine a nested series of disks approaching. You can also imagine more exotic examples like a truncated topologist's sine curve, or simply connected sets with fractal like filaments.

Of course $X$ has to be compact. It seems obvious that $X$ has to be connected and simply connected, though I don't even see a proof of these "easy" facts. In fact, I would bet that $X$ is contractible, though I would have to consider further why that is the case. 

Update: After thinking about it, aside from connectivity, these "obvious" facts are wrong. First, the topologist's sine curve in the plane, defined by $Y=\{(x,\sin(1/x)\,|\,0<x\leq 1\}\cup \{0\}\times[-1,1],$ is a cellular set. I learned this by watching a video of Arunima Ray. You can see it is a cellular set by having each successive closed cell trace out more of the oscillations of the sine curve. Hence, cellular sets do not have to be path connected. Indeed, one can modify this example by considering the topologist's sine curve of revolution in $\mathbb R^3$, formed by rotating $Y$ around the $y$-axis. This space is not simply connected.

As for connectivity, suppose $X$ is separated by sets $U$ and $V$ which are open in the ambient manifold. Then $C_i\subset U\cup V$ for sufficiently large $i$, but then $U\cup V$ would separate $C_i$ which is a contradiction.

Cellular sets seem like they can get pretty pathological, but in fact we have the following theorem.

Theorem: If $X$ is a cellular set in a compact manifold $M$, then the quotient map $\pi\colon M\to M/X$  is a homeomorphism.

The proof of this fact is a beautiful function space argument. One constructs homeomorphisms $h_\epsilon\colon M\to M$ which collapse $X$ to small radius and then argue that the sequence $h_\epsilon$ converges to a limit homeomorphism in the uniform topology on the appropriate function space. 

In a follow-up post, I'll sketch the proof of the topological Schönflies theorem.

 

Prime number theorem

The prime number theorem states that the prime counting function $\pi(x)$ is asymptotically equivalent to $\frac{x}{\log x}$. At some point it might be interesting to delve further into the proof of this fact, but in this post I just wanted to bring about the connection to zeroes of the Riemann zeta function.

To do this, consider the function $\psi(x)=\sum_{p^r<x} \log p,$ originally defined by Chebyshev. It is not that difficult to show that the prime number theorem is equivalent to $\psi(x)\sim x$. It is more difficult to show the following amazing formula of Von Mangoldt

$$\psi(x)=x-\log(2\pi)-\frac{1}{2}\log(1-x^{-2})-\sum_{\zeta(\rho)=0}\frac{x^\rho}{\rho},$$

where the sum is over non-trivial zeroes of the Riemann zeta function. My mind boggles at the beauty and simplicity of this formula!

Now if every zero $\rho =a +ib$ has real part $a<1$, we can see that $|x^\rho|=|x|^a$, so dividing both sides of the Von Mangoldt formula by $x$, each summand of the series approaches $0$ and with a bit more effort one can show that the whole series approaches $0$. Thus PNT follows simply from showing that there are no zeroes of the form $1+ib$.

One reason I find this formula so beautiful is that by plugging in the mysterious zeroes of $\zeta$ one gets closer and closer approximations to $\psi(x).$ See this page for more on that as well as a neat animation.

Chebyshev prime number estimates

 In my previous post, I promised that I would show how to prove that the prime counting function $\pi(n)$ is sandwiched in between two functions of the form $c_1\frac{n}{\log n}$ and $c_2\frac{n}{\log n}$ for some constants $c_1$ and $c_2$. Chebyshev was a pioneer in this area and after a lot of hard work was able to determine in an 1854 paper that this inequality is true for large $n$ and with

$$c_1=0.922\ldots\text{ and }c_2=1.105\ldots.$$

The Prime Number theorem says we should be able to get these constants arbitrarily close to $1$.

Following Lemmermeyer's exposition, we will content ourselves with getting some much cruder estimates, but for far less work! We will show you can take 

$$c_1=\frac{\log 2}{2}\approx 0.347\text{ and }c_2=6\log 2\approx 4.159$$

These are poor estimates but you can get them for some surprisingly little work.

As before the number $N=\binom{2n}{n}$ will play a critical role in the proofs. Let $\nu_p$ be the exponent of $p$ in $N$. In the second proposition of the previous post we showed

$$p^{\nu_p}\leq 2n,$$

which implies $$\nu_p\leq \lfloor\frac{\log 2n}{\log p}\rfloor.$$

Also recall the estimate $\binom{2n}{n}\geq \frac{2^{2n}}{2n}$. Taking logs, we have

$$2n\log 2-\log 2n \leq \log\binom{2n}{n}=\log\left(\prod p^{\nu_p}\right)$$

$$\leq \sum_{p\leq 2n} \left\lfloor\frac{\log 2n}{\log p}\right\rfloor \log p $$

$$\leq \sum_{p\leq 2n} \log 2n=\pi(2n)\log 2n$$

Putting these all together, we get

$$2n\log 2-\log 2n\leq \pi(2n)\log 2n,$$

and rearranging

$$\pi(2n)\geq \log 2\frac{2n}{\log 2n}-1$$

It is not difficult to see that $\log 2\frac{2n}{\log 2n}-1>\frac{\log 2}{2}\frac{2n}{\log 2n}$. This shows $\pi(2n)<\frac{\log 2}{2}\frac{2n}{\log 2n}$. In fact, for general $x$, we need a slightly stronger inequality. Assume $2n\leq x<2n+2$. Then

$$\pi(x)\geq\pi(2n)\geq \log 2\frac{2n}{\log 2n}-1$$

and in order to complete the argument showing $c_1=\frac{\log 2}{2}$, we need to show that

$$\log 2\frac{2n}{\log 2n}-1\geq \frac{(n+1)\log 2}{\log (2n+2)},$$ since this is clearly $\geq \frac{\log 2}{2}\frac{x}{\log x}$.

Lemma: $\log 2\frac{2n}{\log 2n}-1\geq \frac{(n+1)\log 2}{\log (2n+2)}$

Proof:

It suffices to show $\log 2\frac{2n}{\log (2n+2)}-1\geq \frac{(n+1)\log 2}{\log (2n+2)}$. Multiplying both sides by $\log(2n+2)$ we get the following equivalent inequalities:

$$(\log 2)(2n)-\log(2n+2)\geq (n+1)\log 2$$

$$(\log 2)(n-1)\geq \log 2+\log(n+1)$$

$$n\log 2\geq \log(n+1)$$

$$\log 2^n\geq \log (n+1)$$

which is satisfied for all $n\geq 1$. $\Box$

So we have shown $c_1=\frac{\log 2}{2}$. Let us return to the lower bound, which we talked about in the last post. We were zeroing in on Bertrand's Postulate, so the resulting estimate was not the greatest. 

Consider $\prod_{n<p \leq 2n}p.$ This is less than $4^n$ (which follows from the last post.) If we replace each $p$ in this expression with $n$, we get the inequality $n^{\pi(2n)-\pi(n)}\leq\prod_{n<p \leq 2n}p$. So we can conclude

$$n^{\pi(2n)-\pi(n)}<4^n$$

$$\pi(2n)-\pi(n)<\frac{n\log 4}{\log n}$$

Lemma: $\pi(2^k)\leq \frac{3}{k}2^k$.

Proof:

Check it directly for $k\leq 5.$ Now

$$\pi(2^{k+1})-\pi(2^k)<\frac{2^k\log 4}{\log 2^k}=\frac{2^{k+1}}{k},$$

so inductively $$\pi(2^{k+1})=\pi(2^k)+\frac{2^{k+1}}{k}\leq \frac{3}{k}2^k+\frac{2^{k+1}}{k}=\frac{5\cdot 2^k}{k}.$$ It is a simple exercise to show that $\frac{5\cdot 2^k}{k}\leq \frac{3}{k+1}2^{k+1}$, for $k\geq 5$, which will complete the proof. $\Box$.

Now suppose $2^k<x\leq 2^{k+1}$. Then we have

$$\pi(x)\leq \pi(2^{k+1})\leq \frac{6\cdot 2^k}{k+1}\leq \frac{6\cdot 2^k}{k}=6\log 2\frac{2^k}{\log 2^k}\leq 6\log 2\frac{x}{\log x}.$$ So we can take $c_2=6\log 2$.

Bertrand's Postulate

 I have always been fascinated by Bertrand's Postulate, the claim that there is always a prime number between $n$ and $2n$. When I first heard about it as an undergraduate, I tried proving it for myself but was not able to. I'd like to present here what appears to be a now standard proof of the fact due to Erdös. I am basing my exposition on two articles, one by Aigner and Ziegler from the book "Proofs from the book" as well as a PDF found on Franz Lemmermeyer's webpage.  Lemmermeyer has a short section on Bertrand's Postulate which appears to be largely based on the the Aigner-Ziegler article, but he also has a nice exposition on getting some Prime Number Theorem type estimates, which I want to talk about in my next post.

Let $\Theta(x)$ be the product of all primes less than or equal to $x$. More generally, let $\Theta(x,y)$ be the product of all primes in the interval $[x,y]$.

Here's a neat fact:

Proposition: $\Theta(x)\leq 4^{x-1}$ if $x\geq 2$.

Proof: 

Suppose we know the proposition inductively for all values less than $k$. If $k$ is even, then $\Theta(k)=\Theta(k-1)$, so $\Theta(k)<4^{k-1}<4^k$. If $k=2m+1$ is odd, consider the integer $\binom{2m+1}{m}$. By the Binomial theorem, $\binom{2m+1}{m}<2^{2m+1},$ but in fact, because there are two equal terms $\binom{2m+1}{m}=\binom{2m}{m}$, we get the better inequality  

$$\binom{2m+1}{m}<2^{2m}=4^m.$$

Moreover, we claim that $\Theta(m+2,2m+1)$ divides  $\binom{2m+1}{m}$. To see this, notice that $\binom{2m+1}{m}$ can be written as $\frac{(2m+1)(2m)\cdots(m+2)}{m!}.$ Therefore, all of the primes in the range $[m+2,2m+1]$ are factors of the numerator, and they are all larger primes than can be contained in the denominator.

Now $$\Theta(2m+1)=\Theta(m+1)\Theta(m+2,2m+1)$$

$$\leq\Theta(m+1)\binom{2m+1}{m}$$

$$<4^m\cdot 4^m=4^{2m}.$$

This completes the inductive step. $\Box$

In that argument, we used the upper bound $\binom{2m+1}{m}<4^m$. To proceed, we'll want a lower bound on $\binom{2n}{n}$. Consider the set of binomial coefficients $\{\binom{2n}{0},\binom{2n}{1},\ldots\binom{2n}{2n}\}.$ The average value is the sum $2^{2n}$ divided by $2n+1$. Since the middle term is the largest, we get the inequality $\binom{2n}{n}>\frac{2^{2n}}{2n+1}$. In fact, we can improve this to $$\binom{2n}{n}>\frac{2^{2n}}{2n},$$ by considering the set $\{\binom{2n}{0}+\binom{2n}{2n},\binom{2n}{1},\binom{2n}{2},\cdots,\binom{2n}{2n-1}\}$ instead.

Next, we'll look at what we can say about primes appearing in $N=\binom{2n}{n}.$ It is quite easy to check that no prime in the range $(2/3n,n]$ is a divisor of $N$. Think of $n=\frac{(2n)(2n-1)\cdots(n+1)}{n!}$, and notice that $p$ appears once in the denominator, $2p$ appears once in the numerator, while $3p$ is already out of range. A less obvious pair of facts is the following:

Proposition:  

(a) If a prime $p$ satisfies $p>\sqrt{2n}$, then it appears at most once in $\binom{2n}{n}$.

(b) The exponent of $p$ in $\binom{2n}{n}$ is less than or equal to $2n$.

Proof:

The proof  of both these facts uses the Legendre formula that $n!$ contains the prime factor $p$ exactly $\sum_{k\geq 1}\lfloor \frac{n}{p^k}\rfloor$ times. (Incidentally, the Beast Academy 5th Grade Series from AoPS proves this quite nicely, though not using all of the fancy language.)

So the exponent of $p$ in $N$ is given by $$\sum_{k\geq 1}(\lfloor \frac{2n}{p^k}\rfloor-2\lfloor\frac{n}{p^k}\rfloor).$$

It is a simple exercise to check that $\lfloor 2x\rfloor- 2\lfloor x\rfloor\in\{0,1\}$, so each summand is at most $1$. Moreover the summands vanish as soon as $p^k>2n$, so we can say that the exponent of $p$ in $N$ is at most $\operatorname{max}\{k:p^k\leq 2n\}$. This immediately implies statement (b) of the proposition.

To prove (a), notice that if $p>\sqrt{2n}$, the largest power of $p$ less than $2n$ is $1$, so its exponent in $N$ will be at most $1$. $\Box$

Let $\nu_p$ be the exponent of $p$ in $N$. Write

$$N=\prod_{p\leq \sqrt{2n}}p^{\nu_p} \cdot\prod_{\sqrt{2n}<p<\frac{2}{3}n}p^{\nu_p}\cdot\prod_{n<p\leq 2n}p.$$

Notice that we removed the range $[\frac{2}{3}n,n]$ since there are no prime factors of $N$ here. In a similar way to how we argued previously for $\binom{2m+1}{m}$, every prime from $n$ to $2n$ definitely appears with exponent $1$. We can also note that the exponents on the middle term must be either $0$ or $1$. Now we know that $p^{\nu_p}<2n$ from the (b) statement in the above proposition. So we have $$N\leq (2n)^{\sqrt{2n}}\cdot (4^{\frac{2}{3}n})\cdot (2n)^{P(n)},$$

where $P(n)$ denotes the number of primes between $n$ and $2n$. Recall that we are trying to show $P(n)>0$ in order to establish Bertrand's Postulate. On the other hand $N>\frac{4^n}{2n}$, so we get the following inequality:

$$4^{\frac{n}{3}}<(2n)^{\sqrt{2n}+1+P(n)},$$

and taking logarithms, we have

$$P(n)>\frac{2n}{3\log_2(2n)}-(\sqrt{2n}+1).$$

Now we are basically done. The term $\frac{2n}{3\log_2(2n)}$ is going to dominate $\sqrt{2n}+1$ in the long run, and in fact one can check that the right hand side is positive after $n=468$. One can then just verify by hand or computer that $P(n)>0$ for $n<468$. In fact one just needs to find a sequence of primes where each is smaller than twice the previous until you exceed 468. 

It is interesting to notice that this argument establishes a lot more than Bertrand's Postulate. Not only does it show $P(n)>0$, but it give a pretty substantial lower bound of $\frac{2n}{3\log_2(2n)}-(\sqrt{2n}+1)$.

If we let $\pi(n)$ denote the number of primes $\leq n$, clearly $P(n)<\pi(2n)$, so we also get a lower bound on the number of primes.

$$\pi(2n)\geq \frac{2n}{3\log_2(2n)}-(\sqrt{2n}+1).$$

In fact, using the inequality $\sqrt{2n}-1\geq \frac{21}{4}\log_2(2n)$ for $n\geq 2^{11}$, we get

$$\pi(2n)>P(n)>\frac{1}{7}\frac{2n}{\log_2(2n)}.$$

The Prime Number Theorem says that $\pi(n)$ is asymptotically the same as $\frac{n}{\log n}$, so this estimate is in the same ballpark. 

As a step toward the prime number theorem Chebyshev proved that $c_1 \frac{n}{\log n}\leq \pi(n)\leq c_2 \frac{n}{\log n}$ for some positive constants $c_1$ and $c_2$. Our estimate for $P(n)$ gives $c_1=\frac{\ln 2}{7}$. In our next post we will improve $c_1$ and find a $c_2$ that works, following Lemmermeyer's exposition.

A fun argument due to Erdös

 I wanted to share a fun argument that I read in the book Gamma, Exploring Euler's Constant by Julian Havil. Whenever I start reading through proofs of the prime number theorem, this sort of thing crops up and I always find it hard to get my mind to grok it. I understand it line by line, but have difficulty really feeling the essence.

Theorem: The number of primes less than or equal to $n$ is at least $\frac{1}{2}\log_2n.$

The proof is as follows. Let $p_1,p_2,\ldots,p_k$ be the primes less than or equal to $n$. Then any natural number less than or equal to $n$ can be uniquely written as $$p_1^{\epsilon_1}\cdots p_k^{\epsilon_k}\cdot m^2,$$ where the exponents $\epsilon_i\in\{0,1\}$ and $m$ is an integer. In other words, we have decomposed the number into a square free part times a square, and all the exponents of the primes in the square free part will be either $0$ or $1$. Obviously $m^2\leq n$, so $m\leq\sqrt{n}$.

So in order to construct any integer less than or equal to $n$ you must make $k$ binary choices on the exponents of the squarefree part as well as pick an integer $m\leq\sqrt{n}$. Not all such choices will make an integer less than or equal to $n$, but you will definitely construct every integer less than or equal to $n$ in this way! So $$n\leq\sqrt{n}\cdot 2^k.$$ Dividing both sides by $\sqrt{n}$ and taking the base $2$ logarithm, we get the desired result. $\Box$

Well, I am happy I wrote that down, and I have to say I understand the argument a lot better now! I think one of the things that made it a bit hard to grasp at first is that the inequalities are not just about numbers but about numbers of numbers. It is a really nice argument now that I understand it.

Actually, I got to thinking about how this argument fits into a set of "nearby" arguments.  The most straightforward variation is to see what happens when you dispense with the square free factorization. Every number $\leq n$ can be written as $$p_1^{r_1}\cdots p_k^{r_k},$$ but how do we bound the exponents. One try would be to note that the largest any exponent can get is when the base is $2$, where it cannot exceed $\log_2 n$. So, with a very crude estimate, each exponent can be chosen in $\leq \log_2n$ ways, leading to the the inequality $$n\leq (\log_2n)^k.$$ Taking natural logs, we derive 

$$k\geq \frac{\ln n }{\ln(\log_2 n)}. $$

I am kind of surprised this gives as much information as it does. I had expected to get an inequality with no content like $2n>n$! This crude inequality still approaches infinity as $n$ grows, implying the infinitude of primes. When $n=100,000$, we get $k\geq 4.097$ which is a fairly modest statement. On the other hand the original estimate of $1/2 \log_2n$ gives $k\geq 8.304,$ which is not much better. However, one can see that the ratio of the Erdös estimate to the crude estimate goes to infinity, so the Erdös estimate is a substantial improvement.

Another way one could vary the argument is to look at the cube free factorization instead. In that case the exponents all have three choices and there is a factor $m^3$. So we get $$n\leq 3^k\cdot \sqrt[3]{n},$$ leading to the inequality $$k\geq\frac{2}{3}\log_3 n.$$ When $n=100,000$ we get $k\geq 6.986,$ which is again not bad compared to the original, though still worse.  In fact the cubic estimate is always $\frac{4\log2}{3\log 3}\approx 0.84$ times the original estimate. So the original estimate is "only" an improvement of the multiplicative constant.

Chirality and the Conway Polynomial

I'd like to relay a simple conjecture in knot theory, which, despite its simplicity to state, remains unproven. It has to do with the notion of chirality and its opposite amphicheirality. A knot $K\subset S^3$ is said to be amphicheiral if it is isotopic to its mirror image. There are two types. If the strand orientation is preserved, it is said to be $+$ amphicheiral, while if the strand orientation is reversed, it is said to be $-$ amphicheiral. Equivalently in both cases, there exists an orientation-reversing homeomorphism $h\colon S^3\to S^3$ which fixes the knot setwise and which either preserves of reverses the knot's strand orientation. If the knot is hyperbolic, then $h$ can be realized by an isometry, and must be of finite order. One particularly interesting case is when $h$ is an involution. In that case, whether or not $K$ is hyperbolic, we say the knot is strongly $\pm$ amphicheiral. It is possible that a knot can be both $+$ and $-$ amphicheiral simultaneously. The figure $8$ knot has this property.

The Jones Polynomial $V_K(t)\in \mathbb Z[t,t^{-1}]$ is a knot invariant which does not see strand orientation and which satisfies $V_{K^*}(t)=V_{K}(t^{-1})$ where $K^*$ is the mirror image of $K$. Thus the Jones Polynomial of an amphicheiral knot must be symmetric with respect to $t$, and indeed, this is a good way to detect chirality of many knots, e.g. the trefoil.

The Conway Polynomial $C_K(z)\in \mathbb Z[z^2]$ on the other hand is a knot invariant which does not distinguish between $K$ and $K^*$, so it cannot detect chirality in the same way as the Jones. However, we have

Conjecture:
Let $\overline{C_K(z)}\in\mathbb Z_4[z^2]$ be the image of $C_K(z)$ under the natural projection $\mathbb Z[z^2]\twoheadrightarrow \mathbb Z_4[z^2]$. If $K$ is amphicheiral, then
$$\exists f(z)\in \mathbb Z_4[z] \text{ such that }\overline{C_K(z)}=f(z)f(-z)\in \mathbb Z_4[z].$$

(Like a true topologist, I am using $\mathbb Z_k$ to denote the cyclic group of order $k$.)

In this blog post, I would like to explain several confluent lines of evidence for this conjecture. I made an equivalent conjecture in 2006, in a paper with the same name as my blog post. Later, my student Vajira Manathunga showed that my conjecture is equivalent to the one just stated.

The case of strongly amphicheiral knots:
Kawauchi and Hartley proved some time ago that if $K$ is strongly amphicheiral, then $C_K(z)=f(z)f(-z)$ for some $f(z)\in \mathbb Z[z]$. Indeed in the $+$ case, they showed that $C_K(z)=f(z^2)^2$ for some $f$. (They phrased their results in terms of Alexander Polynomials.)

Later Hartley extended this result to the case of all $-$ amphicheiral knots. At this point, the reader may wonder whether the real conjecture ought to be over the integers and not over $\mathbb Z_4$, but it turns out to be false over $\mathbb Z$. Hartley and Kawauchi did not consider the general case in their papers. It may be that they didn't see a reason why it should hold, or perhaps they had a counterexample which they didn't publish. In any event the first published counterexample was due to Ermotti et al, who found a $+$ amphicheiral knot with an orientation reversing symmetry of order $4$, which does not split over $\mathbb Z$.

Independently of the Kawauchi and Hartley story, I was led to conjecture the $\mathbb Z_4$ version after looking at certain properties of Vassiliev invariants related to the Conway polynomial. Gratifyingly the Ermotti et all example does split over $\mathbb Z_4$. Manathunga has since come up with lots of examples which don't split over $\mathbb Z$. See our joint paper on the arXiv. Indeed, in that paper, we prove the conjecture for the case of a positive amphicheiral symmetry that preserves a braid axis, using the Burau representation to calculate the Alexander (and hence Conway) polynomial.

I would like to sketch the argument given by Hartley and Kawauchi for the strongly positive case as well as explain a little bit the Vassiliev invariant background which led to this conjecture.

Review of the Alexander Polynomial:
The Alexander Polynomial is an invariant which is strongly related to the Conway polynomial, and has a very pleasant topological definition. The complement of any knot in $S^3$ has $H_1(S^3\setminus K)=\mathbb Z$, so the covering space $X_K$ of $S^3\setminus K$ associated with the commutator subgroup has deck transformation group $\mathbb Z$. This means that $H_1(S^3\setminus K)$ is a $\mathbb Z[\mathbb Z]$-module. We will think of the group ring $\mathbb Z[\mathbb Z]$ as Laurent polynomials in $t$: $\mathbb Z[t,t^{-1}]$. Moreover, we will switch to rational coefficients, since $\mathbb Q[t^{\pm 1}]$ is a PID. Let $M_K$ be the first homology of $M_K$ with rational coefficients, looked upon as a $\Lambda=\mathbb Q[t^{\pm 1}]$ module. This is called the Alexander Module. One can show, for example by finding an explicit presentation for $M_K$, that $M_K$ is a finitely generated torsion module over $\Lambda$. By the classification of modules over a PID, we have
$$M_K\cong \bigoplus \Lambda/ f_i \Lambda$$
for some $f_i\in \Lambda$.

The Alexander polynomial $\Delta_K(t)$ is defined to be the product of all $p_i$: $\Delta_K(t)=\prod_i f_i$. I.e. it is the order of the Alexander module $M_K$. The polynomials $f_i$ are only well-defined up to units in $\Lambda$. So in addition to multiplying by rational numbers, one can also multiply by powers of $t$. The Alexander polynomial is therefore only well-defined up to these operations. However,  there is a satisfying normal form. It turns out that $\Delta_K(t)$ is equivalent up to multiplying by powers of $t$ to $\Delta_K(t^{-1})$. So one can assume that it is symmetric $\Delta_K(t)=\Delta_K(t^{-1})$. Moreover one can multiply by a constant so that $\Delta_K(1)=1$. With these conventions, though not obviously, $\Delta_K(t)$ has integer coefficients. In fact, Levine showed that each $f_i$ has these properties.

The Conway Polynomial is a repackaging of the Alexander polynomial under the change of variables $z=t^{1/2}-t^{-1/2}$.

Hartley and Kawauchi's Argument for strongly positive amphicheiral knots:
Theorem: If $K$ is strongly positive amphicheiral, then the Alexander module is a direct double $M_k\cong N \oplus N$ for some $\Lambda$-module $N$. This implies that the Conway polynomial is a square $f(z^2)^2$.

Proof: Let $\alpha\colon S^3\to S^3$ be an orientation-reversing involution that preserves the knot's strand orientation. It turns out that it must have two fixed points in the complement of the knot. This implies that $\alpha$'s restriction to $S^3\setminus K$ lifts to an involution $\tilde{\alpha}$ of the universal abelian cover $X_K$. We will use that $M_K$ has a nondegenerate sesquilinear pairing, called the Blanchfield form:
$$B\colon M_K\otimes M_K\to \mathbb Q(t)/\mathbb Q[t,t^{-1}].$$
By definition, $B(x,y)$ is defined as follows. Since $M_K$ is torsion, there is some $\Delta\in \mathbb Q[t,t^{-1}]$ such that $\Delta\cdot y=\partial \sigma$. Indeed one could simply let $\Delta$ be the Alexander polynomial. Then we let
$$B(x,y)=\frac{1}{\Delta}\sum t^{-n}(t^nx)\cdot \sigma,$$
where $\cdot$ denotes the algebraic intersection number.

$B$ is sesquilinear: $B(fx,gy)=f\bar g B(x,y)$ and $B(x,y)=\overline{B(y,x)}$.

Moreover, because $\alpha$ preserves the knot's orientation, it must reverse a meridian. Hence $\tilde{\alpha} t\tilde\alpha=t^{-1}$, which implies that $B(\tilde\alpha x,\tilde\alpha y)=-\overline{B(x,y)}.$ Now the trick is to define $B'(x,y)=B(x,\tilde{\alpha} y)$. It is not hard to see that this is an antisymmetric nondegenerate bilinear form (not sesquilinear) on $M_K$. Now $M_K$ is a finitely generated module over a PID; so it decomposes as a direct sum of pieces $\Lambda/f\Lambda$. Moreover we may assume that the $f$'s are powers of primes $p\in \Lambda$. We can therefore decompose $M_K$ into a direct sum of $p$-primary components: $$M_K=\oplus_p M_K^{(p)}$$ where each $M_K^{(p)}$ is a direct sum of cyclic modules of order $p^k$ for some $k$. It is not hard to see that this must be an orthogonal decomposition with respect to $B'$. Now, each $M_K^{(p)}$ decomposes into a direct sum of cyclic modules of order $p$, plus those of order $p^2$ and so forth. This decomposition is not necessarily orthogonal, but can be fixed up to be so, using an upper triangular argument. So we get summands $\oplus_{i=1}^{j(k,p)} \Lambda/p^k\Lambda$ supporting a nondegenerate antisymmetric bilinear form. Tensoring with the field $F=\Lambda/p\Lambda$, we get a symplectic vector space of dimension $j(k,p)$, which must therefore be even! This implies that $M_K$ is a direct double, as desired. $\Box$



So, now that we understand this classical result of Hartley and Kawauchi, how did Vassiliev invariants enter the picture? For the purposes of this blog post it is not that important to know what exactly a Vassiliev invariant is. It is a type of knot invariant. Some are "Vassiliev" and some are not. The name is used as an adjective to describe some knot invariants in much the way the proper name "Hausdorff" is used as an adjective to describe certain topological spaces. The first thing we want to bring attention to is that the coefficients of the Conway polynomial are Vassiliev invariants of knots:
$$C_K(z)=1+\sum_{i\geq 1}c_{2i}(K)z^{2i},$$
and their Vassiliev degree is $2i$. These invariants are not additive under connected sum, which can be desirable. A standard trick to make them additive is to apply the formal logarithm of power series $\log\colon 1+z\mathbb Z[[z]]\to z\mathbb Z[[z]]$ to $C_K(z)$ and then consider the coefficients. However, this trick does not preserve integrality of the coefficients. This leads to the following definitions:

Definition: Define $$\exp_{\mathbb Z}\colon z\mathbb Z[[z]]\to 1+z\mathbb Z[[z]]$$ by the formula $$\exp_{\mathbb Z}(a_1 z+a_2 z^2+\cdots)= (1-z)^{a_1}(1+z^2)^{a_2}(1-z^3)^{a_3}\cdots=\prod (1+(-z)^n)^{a_n}$$

It is not hard to see that $\exp_{\mathbb Z}$ is bijective and sends addition to multiplication. We define $\log_{\mathbb Z}\colon 1+z\mathbb Z[[z]]\to z\mathbb Z[[z]]$ as its inverse.

Primitive Vassiliev Invariants from the Conway Polynomial: We now define primitive (i.e. additive under connected sum) Vassiliev invariants of degree $2k$ coming from the Conway polynomial, $pc_{2k}\colon \{\text{Knots}\}\to \mathbb Z$, by
$$\log_{\mathbb Z}(C_K(z))=\sum_{i=1}^\infty pc_{2i}(K) z^{2i},$$
where the formal logarithm $\log_{\mathbb Z}$ is taken with respect to the variable $z^2$.

These invariants have an interesting property. Modulo $2$, they drop in degree. $pc_{2n}\pmod 2$ is of degree $2n-1$. On the other hand, no torsion is known among Vassiliev knot invariants. It is thus natural to conjecture

Conjecture: There exist integer-valued Vassiliev knot invariants $v_{2n-1}$ of degree $2n-1$ such that $v_{2n-1}\cong pc_{2n}\pmod 2$.

Such invariants have been explicitly identified for $n=1,2$.

A knot invariant taking values in an abelian group is said to be odd if it changes sign under mirror image. Thus odd $\mathbb Z$-valued knot invariants vanish on amphicheiral knots.

Somewhat unmotivated conjecture: The invariants $v_{4n-1}$ can be chosen to be odd.

This appears to be untrue for general $v_{2n-1}$, but experimental evidence suggested to me that it could be true in the stated half of cases! If this is true, it implies that $pc_{4n}\cong 0\pmod 2$ on amphicheiral knots.

The following proposition then tells us how these considerations on Vassiliev invariants lead to the mod $4$ splitting conjecture.

Proposition: The invariants $pc_{4n}\cong 0\pmod 2$ iff $C_K(z)$ factors as $f(z)f(-z)$ modulo $4$.

So, one approach to proving the mod 4 splitting conjecture is to explicitly identify degree $4n-1$ Vassiliev invariants $v_{4n-1}$ which have the same parity as $pc_{4n}$. This appears to be a difficult problem. Already for  $n=2$ the answer is not known. In degree $7$, the space of all Vassiliev invariants is spanned by polynomial invariants coming from the HOMFLY and Kauffman polynomials, together with a $2$-fold cable. Computer calculations show that $v_{4n-1}$ must involve this cabled invariant, though we have not explicitly identified it.

The Kurosh Problems for PI algebras via Combinatorics

This is a continuing series of posts on some of the highlights of a class at UCSD being given by Efim Zelmanov.

For a field $F$, recall the Kurosh problem as I stated it for algebras.

Kurosh Problem for nilalgebras: Suppose $A$ is a finitely generated $F$-algebra which is nil. That is $\forall a\in A\exists N\in\mathbb Z a^N=0$. Is $A$ finite dimensional (equivalently, nilpotent)?

We saw that the Golod-Shafarevich construction gives a counterexample. However, it turns out that the answer is positive for PI algebras, which we will show. There is another version of the Kurosh problem which I didn't state since it was not relevant for the Burnside problem. It replaces the concept of "nil" with "algebraic."

Definition: An element $a\in A$ is said to be algebraic if there is a polynomial $f(t)$ with $F$-coefficients such that $f(a)=0\in A$. An algebra is said to be algebraic if every element is algebraic.

Kurosh Problem for algebraic algebras: Suppose $A$ is algebraic and finitely generated. Is $A$ finite dimensional?

Since nil implies algebraic, the Golod-Shafarevich construction shows that this is false in general, but we will show that it is indeed true for PI algebras.

Everything follows from a combinatorial lemma about words. But first we need a definition. Let $X$ be a finite alphabet with a fixed order, say $x_1<x_2\ldots x_m$. Order the set of words $X^*$ lexicographically, with the convention that if $v$ is an initial word for $w$, then $v>w$. (This is the opposite of the usual convention.) So for example $x_1x_2<x_1^2$ and $x_1x_2>x_1x_2x_3$.

Definition: A word $w\in X^*$ is $n$-divisible if $w=w'u_1\cdots u_n w''$ for some subwords $u_i$ and every nontrivial permutation $u_{\sigma(1)}\ldots u_{\sigma(n)}$ is a smaller word than $u_1\ldots u_n$. If a word is not $n$-divisible, we say it is $n$-indivisible.

If $A$ satisfies a polynomial identity of degree $n$, then by linearizing we can assume it is of the form $x_1\ldots x_n+\sum_{1\neq\sigma\in \mathcal S_n}\alpha_\sigma x_{\sigma(1)}\cdots x_{\sigma(n)}$. Plugging in $x_i=u_i$, we can rewrite an $n$-divisible word as a linear combination of smaller words. Hence $A$ is spanned by words in its generators which are $n$-indivisible.

Shirshov $\mathcal N(m,k,n)$ lemma: There exists a function $\mathcal N\colon \mathbb N^3\to \mathbb N$ such that any word of length $\geq \mathcal N(m,k,n)$ in the alphabet $X=\{x_1,\ldots, x_m\}$ either has a subword of the form $v^k$ or is $n$-divisible.

Proof of Shirshov lemma:
We proceed by induction on the pair $(n,m)$ ordered lexicographically. For the base cases $m=1$, we can take $\mathcal N(1,k,n)=k$.  So now suppose we know how to define it for all pairs less than $(n,m)$.

Define a set of words
$$T=\{x_m^sx_{i_1}\cdots x_{i_l}\,|\, s\geq 1, \ell\geq 1, \forall j\,i_j<m, s<k, \ell<\mathcal N(m-1,k,n)\}
$$
Note that $T$ is finite. We say a word is a $T$-word if it is a product of elements of $T$. Any word that starts with $x_m$ and ends with some $x_i\neq x_m$, and is both $n$-indivisible and doesn't contain any $k$th powers must be a $T$ word. This is because any such word can be written as a product of $T$-type words except that $s$ and $\ell$ could be large. However, the size of $s$ is bounded by the hypothesis of no $k$th powers and each $\ell$ is inductively bounded by $\mathcal N(m-1,k,n)$ since it does not contain any powers and is $n$-indivisible.

Now we want to show that a long word $v$ in $X$ must contain $k$th powers or be $n$-divisible. So suppose $v$ does not contain $k$th powers and is $n$-indivisible. Any such word $v$ can be written as $v'(x_1,\ldots,x_{m-1})v''x_m^\mu$, where $v''$ is a $T$-word, $\ell(v')<\mathcal N(m-1,k,n)$ and $\mu<k$. So we need to control the length of the $T$-word $v''$.

Now $T$-words can be ordered by looking at the order with respect to $X$ or with respect to $T$. Luckily, or actually by design, these give the same order on $T$-words of the same composition.

Lemma: Given two $T$-words $u,v$ of the same composition (that is they use the same letters with the same multiplicity from the alphabet $X$, or equivalently $T$), then
$$u<_T v\Leftrightarrow u<_X v.$$
Proof of lemma:
If $u$ and $v$ begin with the same letter of $T$, then proceed by induction. Otherwise, suppose that the first $T$ word of $v$ is an initial word of the first $T$ word of $u$. Then, $u<v$ with respect to both orders. $\Box$

Now suppose that $\ell_T(v'')>\mathcal N(|T|,k,n-1)+1$. Then $v''=\underline{\hspace{2em}}t$, where there is an $(n-1)$-division in the initial space. So we need such an $n-1$ division to promote to an $n$-division. Suppose $u_1\cdots u_{n-1}$ is an $n-1$ division, where the $u_i$ are $T$-words.

Lemma: $u_1\cdots u_{n-1} x_m$ is $n$-divisible.
Proof of lemma:
Write $u_i=x_mu'_i$. Divide the word as $$\underbrace{x_m}\underbrace{u_1'x_m}\underbrace{u_2'x_m}\cdots\underbrace{u'_{m-1}x_m}.$$
 Every permutation of these words decreases the order, given our assumption that this is true for the $u_i$'s. $\Box$

So, to finish the proof, we can define $
\begin{multline*}
\mathcal N(m,k,n)=[\mathcal N(m-1,k,n)+k] +\\ [\mathcal N(m-1,k,n)+k]\cdot\mathcal N(m^{\mathcal N(m-1,k,n)},k,n-1)
\end{multline*}
$
Here the first $[\mathcal N(m-1,k,n)+k]$ controls the length of $v'$ and $x_m^\mu$. The second $[\mathcal N(m-1,k,n)+k]$ bounds the length of an element of $T$, while $m^{\mathcal N(m-1,k,n)}$ is an upper bound on the size of $T$.

This completes the proof of the Shirshov Lemma. $\Box$

We actually need a slightly strengthened version of the Shirshov Lemma, which itself requires an elementary lemma.

Lemma: Let $v,w\in X^*$ be words in the alphabet $X$. Then $v,w$ commute iff they are powers of the same word.
Proof of lemma: 
We use induction on the sum of lengths $\ell(v)+\ell(w)$. Assume that $vw=wv$. Then if $\ell(v)=\ell(w)$, it is clear that $v=w$, since they are the initial words on both sides of the equality. On the other hand if $\ell(v)<\ell(w)$, say, then $v$ is an initial word of $w$: $w=vw'$. So $vvw'=vw'v$, and hence $vw'=w'v$. By induction $v=u^k$ and $w'=u^\ell$, and $w=u^{k+\ell}$. $\Box$

Lemma: Let $v$ be a word that is not a nontrivial power. If $\ell(v)>\geq n$, then $v^{2n}$ is $n$-divisible.
Proof of lemma:
Write $v=x_{i_1}\cdots x_{i_\ell}$ for $\ell\geq n$. Consider all cyclic permutations of the word. If two cyclic permutations are equal, then $v'v''=v''v'$ for some decomposition, meaning that $v$ is a proper power by the previous lemma. Now write
$$v=v_1'v_1''=v_2'v_2''=\cdots=v_\ell'v_\ell',$$
which are all the ways of decomposing $v$ into two words. Order them in such a way that
$$v_1''v_1'>v_2''v_2'>\cdots>v_\ell'' v_\ell'.$$
Note that $v^2$ contains all of these as subwords: $v^2=v_i'\underbrace{v_i''v_i'}v_i''$. Consider
$$v^{2n}=v^2\cdots v^2=(v_1'\underbrace{v_1''v_1'v_1'')(v_2'}\underbrace{v_2''v_2'v_2'')(v_3'}v_3''v_3'v_3'')\cdots$$
The bracketed subwords give an $n$-division. $\Box$

This leads to the following strengthening of the Shirshov Lemma:
Theorem: Let $k\geq 2n$. $\exists \mathcal N(m,k,n)$ such that a word $w$ of length $\geq \mathcal N$ is $n$-divisible or contains a subword $v^k$ with $\ell(v)<n$.

We can now deduce some interesting consequences.

Corollary: Let $A=\langle a_1,\cdots, a_m\rangle$ be a finitely generated PI $F$-algebra such that every product $a_{i_1}\ldots a_{i_s}$ for $s<n$ is nilpotent. Then $A$ is a nilpotent algebra, and hence finite dimensional.
Proof of corollary:
Let $k$ be the maximum of all the degrees of nilpotency and $2n$. Then $A^{\mathcal N(m,k,n)}=(0)$. $\Box$
This implies the Kurosh problem for nilalgebras.

Corollary: Suppose that all words $a_{i_1}\cdots a_{i_s}$ $s<n$ are algebraic. Then $A$ is finite dimensional.
Proof of corollary:
$A$ is generated by products of length $<\mathcal N(m,k,n)$. $\Box$

Corollary: Suppose $A$ is a finitely generated nilalgebra where the degrees of nilpotency are bounded. Then $A$ is finite dimensional.
Proof of corollary:
We have $x^n=0$ for all $x\in A$, so it satisfies a polynomial identity. $\Box$

Corollary: Suppose $A$ is a finitely generated algebraic algebra of bounded degree. Then $A$ is finite dimensional.
Proof of corollary:
Say $n$ bounds the degree. Then $S_n([y,x^n],[y,x^{n-1}],\cdots,[y,x])=0$. Why? If $x$ is algebraic, of degree $n$, then $\{x,\ldots, x^n\}$ are linearly dependent, so that the standard polynomial will be $0$. We put the $y'$s in there so the identity is not itself $0$ as a noncomuttatuve polynomial. $\Box$