Good news, everybody! Isaac Goldbring has brought out a book dedicated to ultrafilters, in nonstandard analysis and many other subjects: Ultrafilters Throughout Mathematics. Go convince someone to buy it for you.

I have decided that the best way for me to regularly contribute to this blog is to send it in a certain direction. Since I love nonstandard analysis and want to get better at it, that will be the general theme. To start with, I thought I would take a look at those weird and wonderful things, ultrafilters. Although the focus will mainly be the use of these in nonstandard analysis (NSA, from now on), we will also look at some other uses.

First, a definition:

Definition 1.1. A set of subsets of some set is a filter on if

1. 
2. implies .
3. If and , then .

We are going to consider filters mainly when .

The idea of a filter was first formulated by Cartan, based on Moore’s nets. Let us recap what a net is. Firstly, we recall that a directed set is a set with a reflexive and transitive relation in which every pair of elements has an upper bound.

Definition 1.2. A function is a net on if is a directed set.

Since a sequence on a space is a function , we see that a net is clearly a generalisation of a sequence, and we can use it to capture ideas of convergence, as we will do with filters. A fundamental idea of filters is that each element is “big” in some sense, and that we will consider two things to be equivalent if they agree on an element of the filter. Accordingly, the conditions on a filter reflect this. The empty set is clearly not big, while the whole set must be. The intersections of two big sets must contain some big thing we are interested in, and a set that is bigger than an already big one must itself be considered big.

If we are going to consider two things the same if they agree on an element of the filter (and only then), they can disagree on at most the complement of that set. By implication, the complement of that set can then not also be part of the filter.

Definition 1.3. A filter on a set is an ultrafilter on if for each , either or .

The example usually used to illustrate an ultrafilter is the collection of cofinite subsets of , that is, subsets of with finite complement.

We are not yet done, because there are ultrafilters that can have “small” elements and that will not be suitable for our purposes. Consider the set of subsets of that contain the number 3. Since , any two sets that contain 3 have an intersection that contains 3, and any set containing a set containing 3 must contain 3, this collection is a filter. Furthermore, if I give you a subset of , and its complement, one, and only one of them can contain 3. Therefore, we have an ultrafilter. However, this is not a very useful filter (for our purposes). For instance, if we were using the ultrafilter to compare subsets of , they would be deemed equivalent if both contained 3.

Any ultrafilter of this kind is called a principal ultrafilter. More generally, any ultrafilter generated by a finite set in this way will be deemed such. Since they are clearly no good, we will only deal with non-principal ultrafilters, that is, an ultrafilter that cannot be generated by any finite set. You may also hear people use the term free ultrafilter at times. This can be categorised as a filter that does not contain any finite sets, or as one for which the intersection of all elements is empty. These are equivalent to the condition of being non-principal, and we can use the terms “free” and “non-principal” interchangeably.

You might ask me now to provide an example of a free ultrafilter, which I sadly cannot do. There are all sorts of complicated reasons which I might get into in a later post, but suffice it to say that ZF guarantees that such things exist. However, I cannot write one down. Nevertheless, let us now fix a non-principal ultrafilter and look at some things we can do with it.

Ultraproducts. This is a fundamental concept in Robinsonian nonstandard analysis, and our constructs will rest on these. (We will be going into Nelson’s nonstandard analysis later, which does not require ultraproducts.) Suppose that we have an index set (which we take to be the natural numbers) and a structure for each . These all need to have the same signature, which is something you need not concern yourself overly much with in our constructions. We form the quotient set

,

where is some fixed free ultrafilter. Let (so is an ultrafilter on and for each . If we consider two sequences of reals, they will be identified as the same in the ultraproduct if they agree on an element of . Although we cannot give specific elements of , we can immediately say that two sequences of reals will be equivalent if they differ on only a finite set, since no finite set can be a member of the ultrafilter in question. Also consider the sequence where if is even and zero otherwise as well as the sequence where if is even and zero otherwise. Are these two sequences equivalent? Of course not – the set of indices on which they agree is empty. Now, can we say if one is equivalent to the sequence consisting of 1’s only, and the other to a constant sequence of zeroes? By the definition of our ultrafilter, either the set of even numbers or the set of odd numbers (but not both!) must be a member. If the set of even numbers is, then will be equivalent to a sequence of ones, and will be equivalent to a sequence of zeroes, and vice versa if the set of odd numbers is. However, that is all we can say – in general, we have no idea which of the two sets will be in our ultrafilter.

Now, what about infinitesimals? In our construction, a sequence tending to zero can be considered an equivalence class. Thus, the sequence can be considered as one infinitesimal, and the sequence as another. They clearly won’t be equivalent, since they agree nowhere. Indeed, we can speculate that the first one represents, in some sense, a smaller quantity than the second, since it tends to zero quicker. In this way, we can define the relation < by saying that (both representing equivalence classes here) if the set belongs to . Of course, there is still the same ambiguity as we faced with equality. Suppose that is a sequence with the value 1 on odd values of , and with value for even values of , and let be the sequence with value 1 on odd values of and with value on even values of . Now, we cannot definitely say that , nor can we definitely say that , since we do not know whether our ultrafilter contains the even or the odd numbers. However, we can say for certain that , at least…

In a previous post I wrote about hypergraph complexity, which I encountered in the paper by Conlon, Fox and Zhao [1]. In that same paper, the hypergraph techniques were used to prove a theorem on the existence of arithmetic progressions. In this post, I want to discuss what arithmetic progressions have to do with graphs, which will perhaps give an indication of why the hypergraph approach is useful. To recap, here is the definition of a hypergraph:

Definition 1. Let be a finite set and a positive integer. Define to be the set of all -element subsets of . An –uniform hypergraph on is defined to be any subset .

Note that the usual definition of a graph corresponds to the case in the above definition. One could picture an -uniform hypergraph as a normal graph, if we require that an edge between exists if and only if there is an edge (in the usual sense of a graph) between and for , plus an edge between and .

Using graphs in the context of arithmetic progressions goes back (at least) to the paper of Ruzsa and Szemeredi [2]. In this, they show that the triangle removal lemma implies Roth’s theorem. For this, I am going to follow the discussion in the book [3]. Consider the following theorem of Ajtai and Szemeredi [3, p23]:

Theorem 2. Let . For there exists an such that if and , then there exist , , for some integer .

Getting from Theorem 2 to Roth’s theorem is not hard. Define by letting if , where we assume is some set of integers with density. We assume that is large enough so that has a density of in . We can see this as a graph with the vertices given by , or we can consider as a grid. In the latter case, the “lower left” (under the diagonal) part of the grid contains all possible differences a subset of could generate, and hence each element of must be represented. We will now try to minimise the occurrences of elements of , in order to be able to apply Theorem 2. Looking at one half of the grid (the easiest thing is to draw a small one of your own), we see that 1 occurs times, 2 occurs times, and so on. Clearly then, if we want to minimise the elements of , they must occur as the largest numbers available in the grid. The number of occurrences is then at least (approximately), which means we have occurrences of in the worst-case scenario. This means that we can apply the Ajtai-Szemeredi theorem.

Another way of using graphs for 3-arithmetic progressions is to consider tripartite graphs. For that, we will need the following:

Theorem 3. [2] (Triangle removal lemma for tripartite graphs) Let be finite non-empty sets of vertices, and let be a tripartite graph on these sets of vertices, thus for . Suppose that the number of triangles in this graph does not exceed for some . Then there exists a graph which contains no triangle whatsoever, and such that for .

In this case, denotes a quantity that is bounded in magnitude by , where as . The following discussion is similar to that in https://yufeizhao.wordpress.com/2012/09/07/graph-regularity/.

Now, suppose that is a subset of so that . Suppose we have three elements of . Is there a way of describing these elements that will only be possible if they are in arithmetic progression, in terms of general integers? Of course, one way to do so would be to say that if are all in , we have a 3AP. However, if we want to use the Triangle Removal Lemma and specifically tripartite graphs, we should use 3 quantities to describe a progressions. It is quick to check that if we have integers such that

are all in , we have a 3AP. Suppose now we pick sets

Consider three disjoint vertex sets , each of cardinality , with the vertices of each labelled as . The important thing about these sets is that the difference (as follows) will let us recover all elements of . Let us connect a vertex and if , and if , and and if . (Check that this representation captures all three APs, and no more.)

We now show how the triangle lemma implies Roth’s theorem. We let our graph be constructed from for some , and assume that . Also assume that contains no non-trivial arithmetic progressions. What does this mean for the tripartite graph? Especially, how many triangles does our graph contain? Since there are no non-trivial arithmetic progressions, we must have

.

From this, we can easily see that if we pick , such that , we can set to satisfy the above equations and obtain a triangle. How many ways are there to do this? Given some , we can pick , and . We therefore look at all triples of the form for , . Given that , we should get at least triangles. But how do we apply the triangle removal lemma here? The form it is in does not immediately lend itself to our problem. Let us rephrase it as follows.

Triangle Removal Lemma (2). For each there exists some such that if removing edges does not make the graph triangle-free, then the graph must have more than triangles.

This gives us enough for Roth’s theorem. We have constructed our triangles in a very specific way, so that we know they are the only triangles in the graph, and they do not share any edges. Given the lower bound on the number of triangles in the graph from above, this means that we will need to remove at least edges to make the graph triangle-free. Therefore, removing edges will not make the graph triangle-free, where is slightly smaller than . We must then have some such that the graph contains at least triangles. This is impossible however, since it is easily shown that the upper bound on the number of triangles is also . For large enough, we have a contradiction.

What is the real essence of this proof? It lies in the fact that to have a lot of triangles in such a graph, we would expect many of them to share edgea. This is why we can remove fewer edges than there are triangles in order to make the graph triangle-free. Since the set has density, we know that there will be a significant number of triangles, but not too many, since they have to be disjoint.

References:

1. David Conlon, Jacob Fox and Yufei Zhao: A relative Szemeredi theorem, GAFA Vol.25 (2015)
2. Ruzsa and Szemeredi: Triple systems with no six points carrying three triangles, Combinatorics (Keszthely, 1976), Coll. Math. Soc. J. Bolyai 18 (1978): 939-945.
3. Terence Tao: A variant of the hypergraph removal lemma, Journal of combinatorial theory, Series A 113.7 (2006): 1257-1280.
4. Fan Chung Graham: The Combinatorics of Patterns in Subsets and Graphs, 2004. https://www.math.ucsd.edu/~fan/teach/262/notes/steve/262_book.pdf

Caveat: I am not a Python expert. I’m just messing around with it.

In a previous post, I used Numpy to calculate some eigenvectors of a four-by-four random matrix . Of course, these values were not exact, since we weren’t using symbolic computation. I have been experimenting with Cython recently in order to speed up some of my (very long) Python calculations, and I was wondering if I could use it to write code which gets close to Numpy speed and accuracy in the computation of eigenvalues. Numpy libraries are C-optimised already, so we shouldn’t be surprised if I can’t easily surpass them. However, this is a good exercise and introduction to using Cython for matrices.

Using Cython is slightly more complicated than simply slotting new commands into Python. You have to install Cython, create a .pyx file with your code, and create a setup.py file. This tutorial has everything you need to get started:

https://pythonprogramming.net/introduction-and-basics-cython-tutorial/

We will now implement an algorithm in Numpy to calculate the largest eigenvalue of a matrix, and compare it to standard Numpy functions. This is based on the power method, which works as follows: Given an initial vector , we form , This will gives us an approximation to the eigenvector belonging to the dominant eigenvalue. To get the eigenvalue corresponding to the vector, we observe that the eigenvalue belonging to an eigenvector is given by the Rayleigh quotient:

First, I am going to give the Python code, applied to a simple matrix for which it is easy to calculate the eigenvalues and eigenvectors by hand.

import numpy as np
import time

# This function takes iterates of the matrix multiplied by an initial vector 
# It also gives the Rayleigh quotient at the end 
# We give inputs the matrix x, the initial vector i, and the number of iterations

def largesteig(x,i,n):
    for j in range(n):
        y = np.matmul(x,i)
        i = y
    a = i    
    b = np.matmul(x,a)
    bb = b.transpose()
    bbb = np.matmul(bb,a)
    #print(bbb)
    c = np.matmul(a.transpose(),a)
    d = bbb/c
    return d

# We now choose the matrix to which we'll apply the function       
    
B = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
I = np.array([[1,1,1]])
K = I.transpose()

t1 = time.time()

def eiglist(A,J,n):
    for i in range(n):
        print(largesteig(A,J,i))       

t2 = time.time() 

print('The algorithm took {} seconds'.format(t2-t1))

Note that we really weren’t careful in choosing the initial vector. The output of the iterations is as follows:

[[0.66666667]]
[[2.]]
[[3.33333333]]
[[3.41176471]]
[[3.41414141]]
[[3.41421144]]
[[3.4142135]]
[[3.41421356]]
[[3.41421356]]
[[3.41421356]]
[[3.41421356]]
[[3.41421356]]
[[3.41421356]]
[[3.41421356]]
[[3.41421356]]
The algorithm took 0.003749370574951172 seconds

It seems like this clearly converges, and quite quickly. (Of course, things can still go awry. I have had experiences where the output seems to stabilise, but with more iterations goes completely wrong.) We can compare the speed with Numpy’s usual way of solving for eigenvalues:

import numpy as np
import time

B = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
I = np.array([[1,1,1]])
K = I.transpose()


t1 = time.time()

A = np.linalg.eigvals(B)
a = max(A)

t2 = time.time()

print(a)
print("Numpy took {} seconds.".format(t2-t1))

This give us the output

3.4142135623730914
Numpy took 0.0005006790161132812 seconds.

This is a lot faster, which should not surprise us greatly, since the Numpy code is optimised. We can increase the speed of our power method code by specifying fewer iterations (or preferably, by introducing some convergence criteria), but we still won’t get close to Numpy. Now let us see if we can close the gap by using some Cython. I am not under the illusion that we are going to beat Numpy to the largest eigenvalue, but even getting close should demonstrate the worth of using Cython in your Python code. All we are really going to do is say what types the different variables are. Usually, this is not necessary in Python, because it is “dynamically typed”. While this makes Python easy to read, it unfortunately slows it down quite a lot. The following script will not be the usual .py file, but .pyx instead.

import numpy as np
import time
cimport numpy as np

def largesteig(np.ndarray[np.int_t, ndim = 2] x, np.ndarray[np.int_t, ndim=2] i,int n):
    for j in range(n):
        y = np.matmul(x,i)
        i = y
    b = np.matmul(x,i)
    bb = b.transpose()
    bbb = np.matmul(bb,i)
    #print(bbb)
    c = np.matmul(i.transpose(),i)
    d = bbb/c
    return d

All we have done is specify the matrix properties (integer entries and dimension). Notice that we had to cimport numpy apart from importing numpy as usual. Before we can use it, we have to create the following .py script. Just call it setup.py:

from distutils.core import setup 
from Cython.Build import cythonize
import numpy

setup(ext_modules = cythonize('Power_method_with_Cython.pyx'),include_dirs=[numpy.get_include()])

In the terminal, we then run the command

home@home:pathtoyourfiles$ python setup.py build_ext --inplace

Finally, we have the following Python script to use what we have built on an example:

import Power_method_with_Cython as PmC
import numpy as np
import time

B = np.array([[2,-1,0],[-1,2,-1],[0,-1,2]])
I = np.array([[1,1,1]])
K = I.transpose()

t1 = time.time()

print(PmC.largesteig(B,K,10))

t2 = time.time() 

print('The algorithm took {} seconds'.format(t2-t1))

Once we have run the above, we get the output

[[3.41421356]]
The algorithm took 0.00029778480529785156 seconds

This is actually pretty unbelievable – we have only added a few type declarations in our functions, specifically to do with matrices, and we are already surpassing Numpy speed! The algorithm may still take longer than Numpy though, depending on the number of iterations used. Then there are further issues about convergence of the method, since a solution is not necessarily stable under more iterations. Still, pretty good for a small change in the code. There are even more ways to speed things up, some of which I might add to this post (once I understand them). Interestingly, I’ve discovered that some of the Cython commands which are supposed to speed things up could actually slow them down, but more one that later.

Some useful links:

https://stackoverflow.com/questions/28727717/how-are-numpy-eigenvalues-and-eigenvectors-computed

https://blog.paperspace.com/faster-numpy-array-processing-ndarray-cython/

Duality is an extraordinarily important concept in both linear and nonlinear programming. However, the equivalency of the primal and dual problems is not completely obvious, and difficult to give a geometric interpretation to. However, we can get a pretty good idea of how it works by examining the equivalency of the two problems. The following is mostly based on the notes from AA Ahmadi, but I have expanded them a little. Specifically, I show how you can also solve the dual problem from the primal. (I also recommend these notes for a really intuitive explanation of Farkas’ Lemma.)

Let us consider the following problem as the primal:

where if each component is .

We contend that this problem will have a solution if the dual problem does:

In both the above problems, we let , and an real matrix. When necessary, we shall refer to the primal problem as (P) and the dual as (D). Note that it does not really matter which one you pick as primal or dual, since the dual of the dual has to be the primal again. What we want to do is to show that the primal problem will have a solution as long as the dual one does.

We can start be rephrasing the problem (P) as follows:

The strategy here can be seen as maximising the lower bound for , instead of minimising directly. We can rephrase this as an ostensibly stronger problem. The above will have a solution as long as there exist , with such that

The stronger problem can now be reformulated as

Now, affine functions of will be equal everywhere if and only if their coefficients are equal, hence the problem can be written as

But this is the same as solving the problem

which is the dual problem. Thus, being able to solve the dual problem means we can solve the primal problem. Conversely, suppose that we can solve the primal problem. We phrase the dual problem as

The dual problem can be solved if there exists some such that:

This just means that, since is nonnegative, implies that is also nonnegative, meaning that we have a nonnegative , and we can minimise to be as close to as possible. Once again, we look at the coefficients in the above affine function of , and can set them to be equal. We therefore get the problem

which is of course just the primal problem again.

The usefulness of primality and duality is not limited to spaces such as . In fact, these basic concepts will hold for any vector lattice, which allows us to do things like semidefinite programming.

This post will again not contain anything very advanced, but try to explain a relatively advanced concept by breaking it down into the ideas that led to its formulation. Once again, the star is that wonderful thing, linear algebra. I have gone from really disliking it upon first encounter (as a naive undergrad), to having immense appreciation for the way in which it can illustrate fundamental mathematical concepts. Here, I will first present some computations that could lead to the formulation of the spectral theorem in the finite dimensional case, and in the next installment I will consider instances of the generalisation. This is similar to how I would approach any new mathematical concept, which is first to reduce it to something I have a decent understanding of, before trying to understand it in any generality.

One version of the spectral theorem can be stated as follows:

Theorem. Suppose is a compact self-adjoint operator on a Hilbert space . Then there is an orthonormal basis of which consists of eigenvectors of , and for which each eigenvalue is real.

Recall that a compact operator is a linear operator such that the image of any bounded subset has compact closure. To get a good idea of the Spectral Theorem, we must first get an idea of what a compact self-adjoint operator would be in a finite dimensional setting. For the moment, let us only consider real matrices. Since the adjoint of a matrix is its conjugate transpose, in this case we consider only symmetric matrices. The question of compactness is trivial in this case. Our Hilbert space will just be with the usual inner (dot) product. In this case, the spectral theorem states that we should be able to find an orthonormal basis of consisting of eigenvectors of our matrix , with corresponding real eigenvalues. Let us consider some simple cases.

Since I don’t feel like working this out by hand, I’ll use some Python:

import numpy as np
A = np.array([[3,2,1],[2,2,3],[1,3,5]])
a = np.linalg.eig(A)
print(a) 

This gives us the following output:

(array([ 7.66203005,  2.67899221, -0.34102227]),
 array([[-0.3902192 , -0.83831649,  0.38072882],
        [-0.53490509, -0.13015656, -0.83482681],
        [-0.74940344,  0.52941923,  0.39763018]])

Note that the eigenvectors are given as the columns of the array. Now, it would not be too much to expect that the eigenvectors will give us a basis for . It would be a stretch to expect these to be orthonormal, however. We can test this out quickly:

b = [np.dot(a[1][i],a[1][j]) for i in range(3) for j in range(3)]
print(b)  

(If you look at the code, you’ll see that I actually took the dot product of the rows, not the columns. Why would this have the same result?)

[0.9999999999999998, -5.551115123125783e-17, 2.7755575615628914e-16, -5.551115123125783e-17, 0.9999999999999998, -4.440892098500626e-16, 2.7755575615628914e-16, -4.440892098500626e-16, 1.0000000000000004]

Those values are awfully close to either 0 or 1, and we can suspect that a precise calculation will probably give us those exact answers. But was this just luck, or would it work on any matrix? Let’s generate a random 4 by 4 matrix, which we multiply by its transpose to make it symmetric (where else do we see this product being used to generate Hermitian operators?):

B = np.random.rand(4,4)
BB = np.matmul(B, B.transpose())
bb = np.linalg.eig(BB)
c = [np.dot(bb[1][i],bb[1][j]) for i in range(3) for j in range(3)]
print(c)

[1.0000000000000007, -8.881784197001252e-16, 5.689893001203927e-16, -8.881784197001252e-16, 1.000000000000001, 1.6653345369377348e-16, 5.689893001203927e-16, 1.6653345369377348e-16, 0.9999999999999989]

Once again, we get entries that are practically 0 or 1, and it seems as if NumPy actually gives us the vectors in the desired form. You can get exact results using symbolic computation, for instance using SymPy:

import sympy as sym
A = sym.Matrix([[3,2,1],[2,2,3],[1,3,5]])
A.eigenvects()

This takes surprisingly long however, and doesn’t tell us much we don’t know. We can now formulate a possible finite dimensional version of the spectral theorem:

Conjecture 1. The eigenvectors of a symmetric real matrix form an orthonormal basis of .

Let us first check this in the case where the matrix has full rank, and we have linearly independent eigenvectors. Suppose that are linearly independent eigenvectors of with eigenvalues . For :

and

Since the dot product is commutative, we must have either that or . The latter part is clearly what we desire, but what if ? If we take a look at the wording of the Spectral Theorem, we can see that it does not state that any set of eigenvectors form a basis, but that there is a basis which consists of eigenvectors. So if we have, say, that and , we can form a new vector as follows:

This new vector is also an eigenvector with eigenvalue :

Furthermore, will be orthogonal to . To get orthonormality, we just normalise . (If there are more than two vectors with the same eigenvalue, we just continue applying the Gram-Schmidt process.) We have therefore used the eigenvectors and with the same eigenvalue to find two orthonormal vectors ( and ) with the same eigenvalue as before.

We can see that we must now replace Conjecture 1 with one that is more in line with the statement of the Spectral Theorem:

Conjecture 2. Given a real symmetric matrix , there is a basis of consisting of eigenvectors of .

Now, where could this conjecture go wrong? In the above examples, we had nice, non-zero eigenvalues. So let us consider a matrix in which the columns are linearly dependent, and which has an eigenvalue of 0 with multiplicity 2:

A = np.array([[1,1,1],[1,1,1],[1,1,1]])
np.linalg.eig(A)

(array([-2.22044605e-16,  3.00000000e+00,  0.00000000e+00]),
 array([[-0.81649658,  0.57735027,  0.        ],
        [ 0.40824829,  0.57735027, -0.70710678],
        [ 0.40824829,  0.57735027,  0.70710678]]))

As expected, we only have one non-zero eigenvalue, but we still have three orthonormal eigenvectors. NumPy has kindle normalised them for us, but it is easy to see that our eigenvectors will be multiples of , and . Clearly, these are orthogonal, and once again we have an orthonormal basis. (Think of a geometric reason why this case is actually simpler than the case of 3 distinct non-zero eigenvalues.)

All of this seems to support our conjecture (but is, of course, not yet a proof). Perhaps our conjecture can be made more general? Will it not perhaps hold for arbitrary matrices? Let’s pick an arbitrary integer matrix, generated for me by Python, and get its eigenvalues and eigenvectors:

[[ 1467.   418. -1015.]
 [  418.  1590. -1565.]
 [-1015. -1565. -1560.]]

D = np.linalg.eig(C)
print(D)

(array([-2369.06809548,  1115.59533457,  2750.47276091]), array([[ 0.20557476, -0.80791041,  0.55228596],
       [ 0.34091587,  0.58811029,  0.73341847],
       [ 0.91734148, -0.03751073, -0.39633011]]))

We check to see if we have an orthonormal basis:

E = [np.dot(D[1][i],D[1][j]) for i in range(3) for j in range(3)]
print(E)

[1.0000000000000002, 5.551115123125783e-17, 0.0, 5.551115123125783e-17, 0.9999999999999998, -2.220446049250313e-16, 0.0, -2.220446049250313e-16, 1.0000000000000002]

Within computational error, this is definitely what we want. One can have some more fun with this by generating a matrix and solving for the eigenvectors, but for now, the conjecture looks pretty solid, and we can try to prove it properly (but not here). It is also valuable at this point to look at what the spectral theorem means for the diagonalisability of the matrix.

This has been fairly simple, but the goal is to study the spectral theorem in its full form. I an upcoming post, I will consider some more interesting operators on infinite dimensional spaces.

This is of course a ridiculously simple thing, that many encounter in their first year of university mathematics, and perhaps even in high school. But, like many things I encountered in those days, I kind-of just accepted its existence and moved on. The calculations were easy enough, so I never really had to think about it.

These days I am having a lot of fun examining the undergraduate maths I haven’t touched in years, ever since being forced to do a lot of linear algebra for a research project. From a more mature perspective (I guess), there is so much I missed on in my undergraduate studies by not questioning the things I was learning more. I have a simple question I ask myself now when encountering a mathematical concept: “What the hell is this and why should it be this way?” It might seem a bit belligerent, but I think a bit of aggression is not entirely unwarranted when encountering mathematics. To me, a mathematical concept or entity is not satisfying until I can say where it came from and what it is good for, and I am slightly ashamed of the years I spent just coasting on the definitions of others, never questioning why things should be this way. Perhaps we too easily accept the work that has gone before, in all its abstract glory. Sure, things were done that way for a reason, but you would be so much better off if you knew that reason…

Hence, the cross or vector product in Euclidean space. For now I am going to accept the existence of the dot/scalar product as is, although if you would like a really nice explanation of it, I would suggest 3Blue1Brown‘s very excellent video. (The entire series of videos is worth watching, and should make you appreciate the wonderful world of linear transformations.) The video on the cross product did not entirely satisfy me though, which is why I’ve written this post.

Firstly, why should there be such a thing as a vector product? In other words, how could such a concept arise naturally? I mean, the definition is certainly not obvious. For the moment, let us leave aside the size of the product and simply ask: when will I need to get a vector perpendicular to two other (linearly independent) vectors? As with so much of mathematics, a physical intuition is the best answer. Suppose we have some kind of flow (for now, we can suppose we are dealing with a liquid) through some surface. The flow is not completely uniform, and the surface is not necessarily flat. Supposing we have a function describing the flow at each point in space (in other words, a vector field), can we determine how much is flowing through the surface?

As with so much of calculus, we are going to approximate the surface linearly (assuming the surface is nice enough to do this, of course – let us assume that it is given by a smooth enough function). Therefore, we consider the an approximation to the surface given by tiny parallelograms. By making these small enough, we can estimate the flow through the surface arbitrarily well. Thus, in order to find the total flow, we need some way of determining how much is flowing out of each of our parallelograms. Any part of the flow that is parallel to the plane created by our parallelogram is not relevant, because it is not flowing out. Therefore, we only need to consider the flow with components perpendicular to our parallelogram, in the direction of what we would call a normal vector.

Suppose we know how to determine our tiny parallelograms (for instance, by using partial derivatives of the surface at certain points), and that we know our parallelograms is defined by two vectors, say, and (supposing that these vectors are not parallel). To get a normal vector, we have to get a vector which is perpendicular to both of these. Fortunately, we can assume that we already know about the dot product, and can use it to test for perpendicularity. If there is a vector perpendicular to both of the above, it will have to satisfy the equations

We now have two equations with three unknowns. But if we decide not to care about the length of the vector, we can set any of its components equal to , and solve for the other two. We therefore get the equations

Now it becomes a simple case of solving two equations in two unknowns, and we get

So far it’s not very pretty, but we can make it look nicer by multiplying every vector by :

This is now starting to look very familiar. In order to get the unit normal vector that we require, we simply need to take

(Exercise: Show that the cross product corresponds to a unique linear map from to .)

We’re not quite done with the original form of , though. Supposing once again that we are trying to calculate the flow through some surface by approximating the surface with small parallelograms. To find the flow through the surface then, we have to find the part of the flow perpendicular to the parallelogram and multiply with the area. As if by magic though, we already have the surface area of the parallelogram! Specifically, it is given by the magnitude of the normal vector that we created earlier.

This is something that we are usually taught, and just assume. But why is that magnitude the area of the parallelogram? It is easy enough to get the surface area of a parallelogram in two dimensions; supposing that the parallelogram is determined by the vectors and , the area is given by (show this!) . This looks familiar, at least, and conforms to the determinant method we usually use to get the area.

In three dimensions, the calculation is not very hard either, but let us forget about the cross product for a moment. To get the area of any parallelogram we, of course, only need to multiply the base times the height. In this case, we can take one of the vectors as the base, and find its length in the usual way. To find the height, we just need the dot product and Pythagoras’ theorem.

Specifically, suppose we have the vectors and , since I’m tired of using subscripts. Taking as the base, we have the first term, . To get the height, we first take the length of projected onto :

Using Pythagoras to evaluate the height , we get

The square of the area is then given be

We can see from this that the usual, determinant form of evaluation the area agrees with the more elementary version above. At least we know that the form of the cross product is justified. There is still, to my mind, one deeper issue that needs to be resolved. We now know, formally, that the area of the parallelogram formed by two independent vectors is given by the size of the cross product, but we can also notice that the determinant form of the cross product consists of evaluating the size of three two-dimensional parallelograms. Why should this be? Left to the reader…

I am no Python expert, and have only recently encountered SymPy, for symbolic calculations. I have been trying to do some (relatively simple) matrix calculations, and it has taking me an embarrassingly long time to figure out how to do this. So, I am sharing what I have learned here to help someone else avoid the rather large number of internet searches I had to do to piece it together.

There is a lot out there on how to use SymPy to solve matrix equations of the form . What I am interested in is taking a bunch of given matrices (with numerical values) and constants, performing some operations with an unknown matrix, and setting each entry of the final matrix equal to zero and solving. In other words, suppose we are given matrices and , which are determined beforehand. There is some constant which can be varied (this forms part of an iterative scheme), and an unknown matrix , which is represented purely symbolically, as such:

There is a function , and we want to find the entries of for which

where indicates the zero matrix.

In order to solve an equation in SymPy, you have to declare the “symbols” that you are solving for. Now, defining a matrix symbol in SymPy is easy, but this did not help me in solving for the matrix, and I kept getting an empty output. I needed a way to iteratively declare each entry of the matrix as a symbol, whilst putting them together as a single matrix. This turned out to be the key to the whole thing. First, let us state the preamble:

import sympy as sym
from sympy.interactive import printing
printing.init_printing(use_latex=True)

The last two lines have no effect on the calculations, but they do give the option of displaying your matrices very nicely. Next, let us define some function with which to work:

def matrixfunction(A,B,C,t):
    mf = A - t*C + B.inv()
    return mf

(The final part of the last line is simply how we compute the inverse of .) Here is one of the first things that tripped me up. In SymPy, you should distinguish between operations involving symbolic matrices and usual operations between matrices. For instance, if I were to declare and to be two arbitrary matrices and wanted, for instance, to multiply them, I would use

D = sym.MatrixSymbol('D', 5, 5)
E = sym.MatrixSymbol('E', 5, 5)
sym.MatMul(D,E)

The output to this would be

D*E

and we would be able to see the symbolic entries of this matrix by using

X = sym.MatMul(D,E)
X.as_explicit()

The same holds for MatAdd. However, if you have defined the matrix by declaring all of its entries to be symbols, there does not seem to be a need to use this method, and a simple * can be used for multiplication and for addition.

Our objective is now to set each entry in the matrix obtained from the function “matrixfunction” equal to zero and solve for the unknown matrix . To do so, we define

def solvefor(A,B,t):
    C = Matrix(n,n,sym.symbols('D0:n(0:n)'))
    sol = sym.solve(matrixfunction(A,B,C,t))
    display sol

Of course, in the above, needs to be replaced by an actual numerical value. Once the functions are defined, we can assign values to , and and run “solvefor”. The output will be a set of values assigning the solution value to each entry of . In this case, has to be invertible. Note that SymPy automatically sets the argument of “sym.solve” equal to zero unless otherwise instructed – in this case, it is set equal to the zero matrix. We are also not specifying which symbols to solve for, since SymPy will automatically solve for the free variables here.

This isn’t an ad, and I am not getting money for this. It just really, really works:

I have taken an interest in regularity lemmas recently, starting with the weak regularity lemma of Frieze and Kannan (in order to understand Szemerédi’s theorem on arithmetic progressions), and decided that I needed to understand the proof of the full regularity lemma. The reference I’m using is Reinhard Diestel‘s book Graph Theory. The regularity lemma seems almost magically powerful to me, and has many applications (and generalisations, but we’ll get to that), but the proof is not all that hard. This post won’t present the full proof, just an overview to get the feel of it. (For another perspective on the regularity lemma by someone who actually knows what he is doing, see Yufei Zhao’s blog entry https://yufeizhao.wordpress.com/2012/09/07/graph-regularity/.)

Much like certain proofs of Roth’s theorem, the idea behind the proof of the regularity lemma is to increase a bounded positive quantity by a constant in each iteration. There are only so many times you can do this, and when you stop, what is left must satisfy the conditions of the lemma. Firstly, some definitions.

Definition 1. Let be a graph. We say that a pair of disjoint subsets is -regular if

for any subsets and for which and .

In the above, as usual denotes edge density

between and .

Definition. Let be given and let , , be a graph. We say that a partition of the vertices of the graph is -regular if

1. 
2. 
3. Each pair , is -regular for all but at most of the possible pairs.

The regularity lemma (in one form; there are several equivalent ones) can be stated as follows.

Regularity Lemma. Give some and some , there exists an integer such that any graph of order at least admits an -regular partition where .

The quantity that we will bound is a kind of norm on the graph. If and disjoint, set

where . It should be clear that . In effect, we consider as a bipartite graph. If the graph is complete, will be , and when empty. It does not seem to me that the quantity has a very intuitive meaning. However, it does satisfy the conditions of being bounded and increasing under successive partitions, which is what we’ll need. Suppose that is a partition of and is a partition of . We define

.

 The fundamental fact which we can prove (quite easily) is that, with all objects as above,

If we have a partition of the graph , we define

We can now prove from the above that if is a refinement of the partition , then

We can already start to spy the strategy of the proof here – we have a bounded quantity that increases with further partitions. What we have to prove is that such refinement of partitions will increase by a fixed amount, which means that it can be done only a bounded number of times, given that is bounded by .

Suppose we are now given a partition of , namely , where all the members except have the same size. We regard as a set of members of the partition itself, each a singleton. We set

and define the quantity for our partition as . With our definitions in place, we can start assembling the required lemmas.

Lemma 0. From the Cauchy-Schwarz inequality, we can conclude that

We won’t be using this lemma in this post (hence calling it Lemma 0), but it is essential in working out the details of the following lemmas.

Lemma 1. If are disjoint and not -regular, there are partitions of and of such that

For the proof, we pick the sets and to be two sets of sizes no less than and , respectively, witnessing the irregularity of , and and to be what is left over in and .

Using Lemma 1, we can prove the following crucial lemma.

Lemma 2. Let be such that . If is an -irregular partition of where and , then there is a partition , , of such that and such that

The crucial part here is that the constant increased by does not depend on any of the parameters except . The idea behind the proof of Lemma 2 is to compare all pairs of sets in the first partition, that might be irregular, then use Lemma 1 to partition them. We take a common refinement of all of these partitions, then partition further to obtain sets of equal cardinality, shunting all the vertices left over into .

What astounds me about Lemma 2 is how wasteful it seems to be. We just partition according to Lemma 1 until we get the desired refinement, and yet this is still powerful enough to lead directly to the result we’re after.

We can now carefully pick our constants and iterate Lemma 2 until we obtain the regularity lemma. Pick such that . We will need to determine the constants and , depending on but not on . These depend on how many times we have to implement Lemma 2, and how many further partitions each use of Lemma 2 creates. Because the -value of a partition cannot be greater than , we know that the number of iterations has to be bounded by . Now, at each stage the cardinality of the exceptional set may increase by , if we start with a partition . We want to prevent the size of the exceptional set at any stage from exceeding , and since we start with some with , we have to make sure that iterations of adding do not exceed either. (This is not precise, since the exceptional set will not actually increase by the same factor at each stage, but this will certainly do.) Now must be chosen large enough so that for the initial exceptional set , we have a large so that and . The exceptional set should be allowed up to members, to ensure we have equal cardinality for all the other sets. If is large enough to allow , we have and

for .

What is left is to choose . We choose larger than the maximum of and the number of sets (non-exceptional) that can be generated by applications of the previous lemma, that is, iterations of elements of the partition increasing to elements.

As long as our graph is relatively small, it is easy to show that a suitable partition exists. If the graph has order , where , we choose and let . For , we want to start with a partition which satisfies for Lemma 2 to be valid. Let be as we established above. Now, let be a minimal subset of so that divides . We then partition into sets of equal cardinality. Because , Lemma 2 is applicable, and we simply apply it repeatedly to obtain the partition we are after.