Linear Programming and the Most Affordable Healthy Diet — Part 1

Optimization is by far one of the richest ways to apply computer science and mathematics to the real world. Everybody is looking to optimize something: companies want to maximize profits, factories want to maximize efficiency, investors want to minimize risk, the list just goes on and on. The mathematical tools for optimization are also some of the richest mathematical techniques. They form the cornerstone of an entire industry known as operations research, and advances in this field literally change the world.

The mathematical field is called combinatorial optimization, and the name comes from the goal of finding optimal solutions more efficiently than an exhaustive search through every possibility. This post will introduce the most central problem in all of combinatorial optimization, known as the linear program. Even better, we know how to efficiently solve linear programs, so in future posts we’ll write a program that computes the most affordable diet while meeting the recommended health standard.

Generalizing a Specific Linear Program

Most optimization problems have two parts: an objective function, the thing we want to maximize or minimize, and constraints, rules we must abide by to ensure we get a valid solution. As a simple example you may want to minimize the amount of time you spend doing your taxes (objective function), but you certainly can’t spend a negative amount of time on them (a constraint).

The following more complicated example is the centerpiece of this post. Most people want to minimize the amount of money spent on food. At the same time, one needs to maintain a certain level of nutrition. For males ages 19-30, the United States National Institute for Health recommends 3.7 liters of water per day, 1,000 milligrams of calcium per day, 90 milligrams of vitamin C per day, etc.

We can set up this nutrition problem mathematically, just using a few toy variables. Say we had the option to buy some combination of oranges, milk, and broccoli. Some rough estimates [1] give the following content/costs of these foods. For 0.272 USD you can get 100 grams of orange, containing a total of 53.2mg of calcium, 40mg of vitamin C, and 87g of water. For 0.100 USD you can get 100 grams of whole milk, containing 276mg of calcium, 0mg of vitamin C, and 87g of water. Finally, for 0.381 USD you can get 100 grams of broccoli containing 47mg of calcium, 89.2mg of vitamin C, and 91g of water. Here’s a table summarizing this information:

Nutritional content and prices for 100g of three foods

Food         calcium(mg)     vitamin C(mg)      water(g)   price(USD/100g)
Broccoli     47              89.2               91         0.381
Whole milk   276             0                  87         0.100
Oranges      40              53.2               87         0.272

Some observations: broccoli is more expensive but gets the most of all three nutrients, whole milk doesn’t have any vitamin C but gets a ton of calcium for really cheap, and oranges are a somewhere in between. So you could probably tinker with the quantities and figure out what the cheapest healthy diet is. The problem is what happens when we incorporate hundreds or thousands of food items and tens of nutrient recommendations. This simple example is just to help us build up a nice formality.

So let’s continue doing that. If we denote by b the number of 100g units of broccoli we decide to buy, and m the amount of milk and r the amount of oranges, then we can write the daily cost of food as

\displaystyle \text{cost}(b,m,r) = 0.381 b + 0.1 m + 0.272 r

In the interest of being compact (and again, building toward the general linear programming formulation) we can extract the price information into a single cost vector c = (0.381, 0.1, 0.272), and likewise write our variables as a vector x = (b,m,r). We’re implicitly fixing an ordering on the variables that is maintained throughout the problem, but the choice of ordering doesn’t matter. Now the cost function is just the inner product (dot product) of the cost vector and the variable vector \left \langle c,x \right \rangle. For some reason lots of people like to write this as c^Tx, where c^T denotes the transpose of a matrix, and we imagine that c and x are matrices of size 3 \times 1. I’ll stick to using the inner product bracket notation.

Now for each type of food we get a specific amount of each nutrient, and the sum of those nutrients needs to be bigger than the minimum recommendation. For example, we want at least 1,000 mg of calcium per day, so we require that 1000 \leq 47b + 276m + 40r. Likewise, we can write out a table of the constraints by looking at the columns of our table above.

\displaystyle \begin{matrix} 91b & + & 87m & + & 87r & \geq & 3700 & \text{(water)}\\ 47b & + & 276m & + & 40r & \geq & 1000 & \text{(calcium)} \\ 89.2b & + & 0m & + & 53.2r & \geq & 90 & \text{(vitamin C)} \end{matrix}

In the same way that we extracted the cost data into a vector to separate it from the variables, we can extract all of the nutrient data into a matrix A, and the recommended minimums into a vector v. Traditionally the letter b is used for the minimums vector, but for now we’re using b for broccoli.

A = \begin{pmatrix} 91 & 87 & 87 \\ 47 & 276 & 40 \\ 89.2 & 0 & 53.2 \end{pmatrix}

v = \begin{pmatrix} 3700 \\ 1000 \\ 90 \end{pmatrix}

And now the constraint is that Ax \geq v, where the \geq means “greater than or equal to in every coordinate.” So now we can write down the more general form of the problem for our specific matrices and vectors. That is, our problem is to minimize \left \langle c,x \right \rangle subject to the constraint that Ax \geq v. This is often written in offset form to contrast it with variations we’ll see in a bit:

\displaystyle \text{minimize} \left \langle c,x \right \rangle \\ \text{subject to the constraint } Ax \geq v

In general there’s no reason you can’t have a “negative” amount of one variable. In this problem you can’t buy negative broccoli, so we’ll add the constraints to ensure the variables are nonnegative. So our final form is

\displaystyle \text{minimize} \left \langle c,x \right \rangle \\ \text{subject to } Ax \geq v \\ \text{and } x \geq 0

In general, if you have an m \times n matrix A, a “minimums” vector v \in \mathbb{R}^m, and a cost vector c \in \mathbb{R}^n, the problem of finding the vector x that minimizes the cost function while meeting the constraints is called a linear programming problem or simply a linear program.

To satiate the reader’s burning curiosity, the solution for our calcium/vitamin C problem is roughly x = (1.01, 41.47, 0). That is, you should have about 100g of broccoli and 4.2kg of milk (like 4 liters), and skip the oranges entirely. The daily cost is about 4.53 USD. If this seems awkwardly large, it’s because there are cheaper ways to get water than milk.

100-grams-broccoli

100g of broccoli (image source: 100-grams.blogspot.com)

[1] Water content of fruits and veggiesFood costs in March 2014 in the midwest, and basic known facts about the water density/nutritional content of various foods.

Duality

Now that we’ve seen the general form a linear program and a cute example, we can ask the real meaty question: is there an efficient algorithm that solves arbitrary linear programs? Despite how widely applicable these problems seem, the answer is yes!

But before we can describe the algorithm we need to know more about linear programs. For example, say you have some vector x which satisfies your constraints. How can you tell if it’s optimal? Without such a test we’d have no way to know when to terminate our algorithm. Another problem is that we’ve phrased the problem in terms of minimization, but what about problems where we want to maximize things? Can we use the same algorithm that finds minima to find maxima as well?

Both of these problems are neatly answered by the theory of duality. In mathematics in general, the best way to understand what people mean by “duality” is that one mathematical object uniquely determines two different perspectives, each useful in its own way. And typically a duality theorem provides one with an efficient way to transform one perspective into the other, and relate the information you get from both perspectives. A theory of duality is considered beautiful because it gives you truly deep insight into the mathematical object you care about.

In linear programming duality is between maximization and minimization. In particular, every maximization problem has a unique “dual” minimization problem, and vice versa. The really interesting thing is that the variables you’re trying to optimize in one form correspond to the contraints in the other form! Here’s how one might discover such a beautiful correspondence. We’ll use a made up example with small numbers to make things easy.

So you have this optimization problem

\displaystyle \begin{matrix}  \text{minimize} & 4x_1+3x_2+9x_3 & \\  \text{subject to} & x_1+x_2+x_3 & \geq 6 \\  & 2x_1+x_3 & \geq 2 \\  & x_2+x_3 & \geq 1 & \\ & x_1,x_2,x_3 & \geq 0 \end{matrix}

Just for giggles let’s write out what A and c are.

\displaystyle A = \begin{pmatrix} 1 & 1 & 1 \\ 2 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix}, c = (4,3,9), v = (6,2,1)

Say you want to come up with a lower bound on the optimal solution to your problem. That is, you want to know that you can’t make 4x_1 + 3x_2 + 9x_3 smaller than some number m. The constraints can help us derive such lower bounds. In particular, every variable has to be nonnegative, so we know that 4x_1 + 3x_2 + 9x_3 \geq x_1 + x_2 + x_3 \geq 6, and so 6 is a lower bound on our optimum. Likewise,

\displaystyle \begin{aligned}4x_1+3x_2+9x_3 & \geq 4x_1+4x_3+3x_2+3x_3 \\ &=2(2x_1 + x_3)+3(x_2+x_3) \\ & \geq 2 \cdot 2 + 3 \cdot 1 \\ &=7\end{aligned}

and that’s an even better lower bound than 6. We could try to write this approach down in general: find some numbers y_1, y_2, y_3 that we’ll use for each constraint to form

\displaystyle y_1(\text{constraint 1}) + y_2(\text{constraint 2}) + y_3(\text{constraint 3})

To make it a valid lower bound we need to ensure that the coefficients of each of the x_i are smaller than the coefficients in the objective function (i.e. that the coefficient of x_1 ends up less than 4). And to make it the best lower bound possible we want to maximize what the right-hand-size of the inequality would be: y_1 6 + y_2 2 + y_3 1. If you write out these equations and the constraints you get our “lower bound” problem written as

\displaystyle \begin{matrix} \text{maximize} & 6y_1 + 2y_2 + y_3 & \\ \text{subject to} & y_1 + 2y_2 & \leq 4 \\ & y_1 + y_3 & \leq 3 \\ & y_1+y_2 + y_3 & \leq 9 \\ & y_1,y_2,y_3 & \geq 0 \end{matrix}

And wouldn’t you know, the matrix providing the constraints is A^T, and the vectors c and v switched places.

\displaystyle A^T = \begin{pmatrix} 1 & 2 & 0 \\ 1 & 0 & 1 \\ 1 & 1 & 1 \end{pmatrix}

This is no coincidence. All linear programs can be transformed in this way, and it would be a useful exercise for the reader to turn the above maximization problem back into a minimization problem by the same technique (computing linear combinations of the constraints to make upper bounds). You’ll be surprised to find that you get back to the original minimization problem! This is part of what makes it “duality,” because the dual of the dual is the original thing again. Often, when we fix the “original” problem, we call it the primal form to distinguish it from the dual form. Usually the primal problem is the one that is easy to interpret.

(Note: because we’re done with broccoli for now, we’re going to use b to denote the constraint vector that used to be v.)

Now say you’re given the data of a linear program for minimization, that is the vectors c, b and matrix A for the problem, “minimize \left \langle c, x \right \rangle subject to Ax \geq b; x \geq 0.” We can make a general definition: the dual linear program is the maximization problem “maximize \left \langle b, y \right \rangle subject to A^T y \leq c, y \geq 0.” Here y is the new set of variables and the superscript T denotes the transpose of the matrix. The constraint for the dual is often written y^T A \leq c^T, again identifying vectors with a single-column matrices, but I find the swamp of transposes pointless and annoying (why do things need to be columns?).

Now we can actually prove that the objective function for the dual provides a bound on the objective function for the original problem. It’s obvious from the work we’ve done, which is why it’s called the weak duality theorem.

Weak Duality Theorem: Let c, A, b be the data of a linear program in the primal form (the minimization problem) whose objective function is \left \langle c, x \right \rangle. Recall that the objective function of the dual (maximization) problem is \left \langle b, y \right \rangle. If x,y are feasible solutions (satisfy the constraints of their respective problems), then

\left \langle b, y \right \rangle \leq \left \langle c, x \right \rangle

In other words, the maximum of the dual is a lower bound on the minimum of the primal problem and vice versa. Moreover, any feasible solution for one provides a bound on the other.

Proof. The proof is pleasingly simple. Just inspect the quantity \left \langle A^T y, x \right \rangle = \left \langle y, Ax \right \rangle. The constraints from the definitions of the primal and dual give us that

\left \langle y, b \right \rangle \leq \left \langle y, Ax \right \rangle = \left \langle A^Ty, x \right \rangle \leq \left \langle c,x \right \rangle

The inequalities follow from the linear algebra fact that if the u in \left \langle u,v \right \rangle is nonnegative, then you can only increase the size of the product by increasing the components of v. This is why we need the nonnegativity constraints.

In fact, the world is much more pleasing. There is a theorem that says the two optimums are equal!

Strong Duality Theorem: If there are any solutions x,y to the primal (minimization) problem and the dual (maximization) problem, respectively, then the two problems also have optimal solutions x^*, y^*, and two candidate solutions x^*, y^* are optimal if and only if they produce equal objective values \left \langle c, x^* \right \rangle = \left \langle y^*, b \right \rangle.

The proof of this theorem is a bit more convoluted than the weak duality theorem, and the key technique is a lemma of Farkas and its variations. See the second half of these notes for a full proof. The nice thing is that this theorem gives us a way to tell if an algorithm to solve linear programs is done: maintain a pair of feasible solutions to the primal and dual problems, improve them by some rule, and stop when the two solutions give equal objective values. The hard part, then, is finding a principled and guaranteed way to improve a given pair of solutions.

On the other hand, you can also prove the strong duality theorem by inventing an algorithm that provably terminates. We’ll see such an algorithm, known as the simplex algorithm in the next post. Sneak peek: it’s a lot like Gaussian elimination. Then we’ll use the algorithm (or an equivalent industry-strength version) to solve a much bigger nutrition problem.

In fact, you can do a bit better than the strong duality theorem, in terms of coming up with a stopping condition for a linear programming algorithm. You can observe that an optimal solution implies further constraints on the relationship between the primal and the dual problems. In particular, this is called the complementary slackness conditions, and they essentially say that if an optimal solution to the primal has a positive variable then the corresponding constraint in the dual problem must be tight (is an equality) to get an optimal solution to the dual. The contrapositive says that if some constraint is slack, or a strict inequality, then either the corresponding variable is zero or else the solution is not optimal. More formally,

Theorem (Complementary Slackness Conditions): Let A, c, b be the data of the primal form of a linear program, “minimize \left \langle c, x \right \rangle subject to Ax \geq b, x \geq 0.” Then x^*, y^* are optimal solutions to the primal and dual problems if any only if all of the following conditions hold.

  • x^*, y^* are both feasible for their respective problems.
  • Whenever x^*_i > 0 the corresponding constraint A^T_i y^* = c_i is an equality.
  • Whenever y^*_j > 0 the corresponding constraint A_j x^* = b_j is an equality.

Here we denote by M_i the i-th row of the matrix M and v_i to denote the i-th entry of a vector. Another way to write the condition using vectors instead of English is

\left \langle x^*, A^T y^* - c \right \rangle = 0
\left \langle y^*, Ax^* - b \right \rangle

The proof follows from the duality theorems, and just involves pushing around some vector algebra. See section 6.2 of these notes.

One can interpret complementary slackness in linear programs in a lot of different ways. For us, it will simply be a termination condition for an algorithm: one can efficiently check all of these conditions for the nonzero variables and stop if they’re all satisfied or if we find a variable that violates a slackness condition. Indeed, in more mature optimization analyses, the slackness condition that is more egregiously violated can provide evidence for where a candidate solution can best be improved. For a more intricate and detailed story about how to interpret the complementary slackness conditions, see Section 4 of these notes by Joel Sobel.

Finally, before we close we should note there are geometric ways to think about linear programming. I have my preferred visualization in my head, but I have yet to find a suitable animation on the web that replicates it. Here’s one example in two dimensions. The set of constraints define a convex geometric region in the plane

The constraints define a convex area of "feasible solutions." Image source: Wikipedia.

The constraints define a convex area of “feasible solutions.” Image source: Wikipedia.

Now the optimization function f(x) = \left \langle c,x \right \rangle is also a linear function, and if you fix some output value y = f(x) this defines a line in the plane. As y changes, the line moves along its normal vector (that is, all these fixed lines are parallel). Now to geometrically optimize the target function, we can imagine starting with the line f(x) = 0, and sliding it along its normal vector in the direction that keeps it in the feasible region. We can keep sliding it in this direction, and the maximum of the function is just the last instant that this line intersects the feasible region. If none of the constraints are parallel to the family of lines defined by f, then this is guaranteed to occur at a vertex of the feasible region. Otherwise, there will be a family of optima lying anywhere on the line segment of last intersection.

In higher dimensions, the only change is that the lines become affine subspaces of dimension n-1. That means in three dimensions you’re sliding planes, in four dimensions you’re sliding 3-dimensional hyperplanes, etc. The facts about the last intersection being a vertex or a “line segment” still hold. So as we’ll see next time, successful algorithms for linear programming in practice take advantage of this observation by efficiently traversing the vertices of this convex region. We’ll see this in much more detail in the next post.

Until then!

About these ads

(Finite) Fields — A Primer

So far on this blog we’ve given some introductory notes on a few kinds of algebraic structures in mathematics (most notably groups and rings, but also monoids). Fields are the next natural step in the progression.

If the reader is comfortable with rings, then a field is extremely simple to describe: they’re just commutative rings with 0 and 1, where every nonzero element has a multiplicative inverse. We’ll give a list of all of the properties that go into this “simple” definition in a moment, but an even more simple way to describe a field is as a place where “arithmetic makes sense.” That is, you get operations for +,-, \cdot , / which satisfy the expected properties of addition, subtraction, multiplication, and division. So whatever the objects in your field are (and sometimes they are quite weird objects), they behave like usual numbers in a very concrete sense.

So here’s the official definition of a field. We call a set F a field if it is endowed with two binary operations addition (+) and multiplication (\cdot, or just symbol juxtaposition) that have the following properties:

  • There is an element we call 0 which is the identity for addition.
  • Addition is commutative and associative.
  • Every element a \in F has a corresponding additive inverse b (which may equal a) for which a + b = 0.

These three properties are just the axioms of a (commutative) group, so we continue:

  • There is an element we call 1 (distinct from 0) which is the identity for multiplication.
  • Multiplication is commutative and associative.
  • Every nonzero element a \in F has a corresponding multiplicative inverse b (which may equal a) for which ab = 1.
  • Addition and multiplication distribute across each other as we expect.

If we exclude the existence of multiplicative inverses, these properties make F a commutative ring, and so we have the following chain of inclusions that describes it all

\displaystyle \textup{Fields} \subset \textup{Commutative Rings} \subset \textup{Rings} \subset \textup{Commutative Groups} \subset \textup{Groups}

The standard examples of fields are the real numbers \mathbb{R}, the rationals \mathbb{Q}, and the complex numbers \mathbb{C}. But of course there are many many more. The first natural question to ask about fields is: what can they look like?

For example, can there be any finite fields? A field F which as a set has only finitely many elements?

As we saw in our studies of groups and rings, the answer is yes! The simplest example is the set of integers modulo some prime p. We call them \mathbb{Z} / p \mathbb{Z}, or sometimes just \mathbb{Z}/p for short, and let’s rederive what we know about them now.

As a set, \mathbb{Z}/p consists of the integers \left \{ 0, 1, \dots, p-1 \right \}. The addition and multiplication operations are easy to define, they’re just usual addition and multiplication followed by a modulus. That is, we add by a + b \mod p and multiply with ab \mod p. This thing is clearly a commutative ring (because the integers form a commutative ring), so to show this is a field we need to show that everything has a multiplicative inverse.

There is a nice fact that allows us to do this: an element a has an inverse if and only if the only way for it to divide zero is the trivial way 0a = 0. Here’s a proof. For one direction, suppose a divides zero nontrivially, that is there is some c \neq 0 with ac = 0. Then if a had an inverse b, then 0 = b(ac) = (ba)c = c, but that’s very embarrassing for c because it claimed to be nonzero. Now suppose a only divides zero in the trivial way. Then look at all possible ways to multiply a by other nonzero elements of F. No two can give you the same result because if ax = ay then (without using multiplicative inverses) a(x-y) = 0, but we know that a can only divide zero in the trivial way so x=y. In other words, the map “multiplication by a” is injective. Because the set of nonzero elements of F is finite you have to hit everything (the map is in fact a bijection), and some x will give you ax = 1.

Now let’s use this fact on \mathbb{Z}/p in the obvious way. Since p is a prime, there are no two smaller numbers a, b < p so that ab = p. But in \mathbb{Z}/p the number p is equivalent to zero (mod p)! So \mathbb{Z}/p has no nontrivial zero divisors, and so every element has an inverse, and so it’s a finite field with p elements.

The next question is obvious: can we get finite fields of other sizes? The answer turns out to be yes, but you can’t get finite fields of any size. Let’s see why.

Characteristics and Vector Spaces

Say you have a finite field k (lower-case k is the standard letter for a field, so let’s forget about F). Beacuse the field is finite, if you take 1 and keep adding it to itself you’ll eventually run out of field elements. That is, n = 1 + 1 + \dots + 1 = 0 at some point. How do I know it’s zero and doesn’t keep cycling never hitting zero? Well if at two points n = m \neq 0, then n-m = 0 is a time where you hit zero, contradicting the claim.

Now we define \textup{char}(k), the characteristic of k, to be the smallest n (sums of 1 with itself) for which n = 0. If there is no such n (this can happen if k is infinite, but doesn’t always happen for infinite fields), then we say the characteristic is zero. It would probably make more sense to say the characteristic is infinite, but that’s just the way it is. Of course, for finite fields the characteristic is always positive. So what can we say about this number? We have seen lots of example where it’s prime, but is it always prime? It turns out the answer is yes!

For if ab = n = \textup{char}(k) is composite, then by the minimality of n we get a,b \neq 0, but ab = n = 0. This can’t happen by our above observation, because being a zero divisor means you have no inverse! Contradiction, sucker.

But it might happen that there are elements of k that can’t be written as 1 + 1 + \dots + 1 for any number of terms. We’ll construct examples in a minute (in fact, we’ll classify all finite fields), but we already have a lot of information about what those fields might look like. Indeed, since every field has 1 in it, we just showed that every finite field contains a smaller field (a subfield) of all the ways to add 1 to itself. Since the characteristic is prime, the subfield is a copy of \mathbb{Z}/p for p = \textup{char}(k). We call this special subfield the prime subfield of k.

The relationship between the possible other elements of k and the prime subfield is very neat. Because think about it: if k is your field and F is your prime subfield, then the elements of k can interact with F just like any other field elements. But if we separate k from F (make a separate copy of F), and just think of k as having addition, then the relationship with F is that of a vector space! In fact, whenever you have two fields k \subset k', the latter has the structure of a vector space over the former.

Back to finite fields, k is a vector space over its prime subfield, and now we can impose all the power and might of linear algebra against it. What’s it’s dimension? Finite because k is a finite set! Call the dimension m, then we get a basis v_1, \dots, v_m. Then the crucial part: every element of k has a unique representation in terms of the basis. So they are expanded in the form

\displaystyle f_1v_1 + \dots + f_mv_m

where the f_i come from F. But now, since these are all just field operations, every possible choice for the f_i has to give you a different field element. And how many choices are there for the f_i? Each one has exactly |F| = \textup{char}(k) = p. And so by counting we get that k has p^m many elements.

This is getting exciting quickly, but we have to pace ourselves! This is a constraint on the possible size of a finite field, but can we realize it for all choices of p, m? The answer is again yes, and in the next section we’ll see how.  But reader be warned: the formal way to do it requires a little bit of familiarity with ideals in rings to understand the construction. I’ll try to avoid too much technical stuff, but if you don’t know what an ideal is, you should expect to get lost (it’s okay, that’s the nature of learning new math!).

Constructing All Finite Fields

Let’s describe a construction. Take a finite field k of characteristic p, and say you want to make a field of size p^m. What we need to do is construct a field extension, that is, find a bigger field containing k so that the vector space dimension of our new field over k is exactly m.

What you can do is first form the ring of polynomials with coefficients in k. This ring is usually denoted k[x], and it’s easy to check it’s a ring (polynomial addition and multiplication are defined in the usual way). Now if I were speaking to a mathematician I would say, “From here you take an irreducible monic polynomial p(x) of degree m, and quotient your ring by the principal ideal generated by p. The result is the field we want!”

In less compact terms, the idea is exactly the same as modular arithmetic on integers. Instead of doing arithmetic with integers modulo some prime (an irreducible integer), we’re doing arithmetic with polynomials modulo some irreducible polynomial p(x). Now you see the reason I used p for a polynomial, to highlight the parallel thought process. What I mean by “modulo a polynomial” is that you divide some element f in your ring by p as much as you can, until the degree of the remainder is smaller than the degree of p(x), and that’s the element of your quotient. The Euclidean algorithm guarantees that we can do this no matter what k is (in the formal parlance, k[x] is called a Euclidean domain for this very reason). In still other words, the “quotient structure” tells us that two polynomials f, g \in k[x] are considered to be the same in k[x] / p if and only if f - g is divisible by p. This is actually the same definition for \mathbb{Z}/p, with polynomials replacing numbers, and if you haven’t already you can start to imagine why people decided to study rings in general.

Let’s do a specific example to see what’s going on. Say we’re working with k = \mathbb{Z}/3 and we want to compute a field of size 27 = 3^3. First we need to find a monic irreducible polynomial of degree 3. For now, I just happen to know one: p(x) = x^3 - x + 1. In fact, we can check it’s irreducible, because to be reducible it would have to have a linear factor and hence a root in \mathbb{Z}/3. But it’s easy to see that if you compute p(0), p(1), p(2) and take (mod 3) you never get zero.

So I’m calling this new ring

\displaystyle \frac{\mathbb{Z}/3[x]}{(x^3 - x + 1)}

It happens to be a field, and we can argue it with a whole lot of ring theory. First, we know an irreducible element of this ring is also prime (because the ring is a unique factorization domain), and prime elements generate maximal ideals (because it’s a principal ideal domain), and if you quotient by a maximal ideal you get a field (true of all rings).

But if we want to avoid that kind of argument and just focus on this ring, we can explicitly construct inverses. Say you have a polynomial f(x), and for illustration purposes we’ll choose f(x) = x^4 + x^2 - 1. Now in the quotient ring we could do polynomial long division to find remainders, but another trick is just to notice that the quotient is equivalent to the condition that x^3 = x - 1. So we can reduce f(x) by applying this rule to x^4 = x^3 x to get

\displaystyle f(x) = x^2 + x(x-1) - 1 = 2x^2 - x - 1

Now what’s the inverse of f(x)? Well we need a polynomial g(x) = ax^2 + bx + c whose product with f gives us something which is equivalent to 1, after you reduce by x^3 - x + 1. A few minutes of algebra later and you’ll discover that this is equivalent to the following polynomial being identically 1

\displaystyle (a-b+2c)x^2 + (-3a+b-c)x + (a - 2b - 2c) = 1

In other words, we get a system of linear equations which we need to solve:

\displaystyle \begin{aligned} a & - & b & + & 2c & = 0 \\ -3a & + & b & - & c &= 0 \\ a & - & 2b & - & 2c &= 1 \end{aligned}

And from here you can solve with your favorite linear algebra techniques. This is a good exercise for working in fields, because you get to abuse the prime subfield being characteristic 3 to say terrifying things like -1 = 2 and 6b = 0. The end result is that the inverse polynomial is 2x^2 + x + 1, and if you were really determined you could write a program to compute these linear systems for any input polynomial and ensure they’re all solvable. We prefer the ring theoretic proof.

In any case, it’s clear that taking a polynomial ring like this and quotienting by a monic irreducible polynomial gives you a field. We just control the size of that field by choosing the degree of the irreducible polynomial to our satisfaction. And that’s how we get all finite fields!

One Last Word on Irreducible Polynomials

One thing we’ve avoided is the question of why irreducible monic polynomials exist of all possible degrees m over any \mathbb{Z}/p (and as a consequence we can actually construct finite fields of all possible sizes).

The answer requires a bit of group theory to prove this, but it turns out that the polynomial x^{p^m} - x has all degree m monic irreducible polynomials as factors. But perhaps a better question (for computer scientists) is how do we work over a finite field in practice? One way is to work with polynomial arithmetic as we described above, but this has some downsides: it requires us to compute these irreducible monic polynomials (which doesn’t sound so hard, maybe), to do polynomial long division every time we add, subtract, or multiply, and to compute inverses by solving a linear system.

But we can do better for some special finite fields, say where the characteristic is 2 (smells like binary) or we’re only looking at F_{p^2}. The benefit there is that we aren’t forced to use polynomials. We can come up with some other kind of structure (say, matrices of a special form) which happens to have the same field structure and makes computing operations relatively painless. We’ll see how this is done in the future, and see it applied to cryptography when we continue with our series on elliptic curve cryptography.

Until then!

Fixing Bugs in “Computing Homology”

A few awesome readers have posted comments in Computing Homology to the effect of, “Your code is not quite correct!” And they’re right! Despite the almost year since that post’s publication, I haven’t bothered to test it for more complicated simplicial complexes, or even the basic edge cases! When I posted it the mathematics just felt so solid to me that it had to be right (the irony is rich, I know).

As such I’m apologizing for my lack of rigor and explaining what went wrong, the fix, and giving some test cases. As of the publishing of this post, the Github repository for Computing Homology has been updated with the correct code, and some more examples.

The main subroutine was the simultaneousReduce function which I’ll post in its incorrectness below

def simultaneousReduce(A, B):
   if A.shape[1] != B.shape[0]:
      raise Exception("Matrices have the wrong shape.")

   numRows, numCols = A.shape # col reduce A

   i,j = 0,0
   while True:
      if i >= numRows or j >= numCols:
         break

      if A[i][j] == 0:
         nonzeroCol = j
         while nonzeroCol < numCols and A[i,nonzeroCol] == 0:
            nonzeroCol += 1

         if nonzeroCol == numCols:
            j += 1
            continue

         colSwap(A, j, nonzeroCol)
         rowSwap(B, j, nonzeroCol)

      pivot = A[i,j]
      scaleCol(A, j, 1.0 / pivot)
      scaleRow(B, j, 1.0 / pivot)

      for otherCol in range(0, numCols):
         if otherCol == j:
            continue
         if A[i, otherCol] != 0:
            scaleAmt = -A[i, otherCol]
            colCombine(A, otherCol, j, scaleAmt)
            rowCombine(B, j, otherCol, -scaleAmt)

      i += 1; j+= 1

   return A,B

It’s a beast of a function, and the persnickety detail was just as beastly: this snippet should have an i += 1 instead of a j.

if nonzeroCol == numCols:
   j += 1
   continue

This is simply what happens when we’re looking for a nonzero entry in a row to use as a pivot for the corresponding column, but we can’t find one and have to move to the next row. A stupid error on my part that would be easily caught by proper test cases.

The next mistake is a mathematical misunderstanding. In short, the simultaneous column/row reduction process is not enough to get the \partial_{k+1} matrix into the right form! Let’s see this with a nice example, a triangulation of the Mobius band. There are a number of triangulations we could use, many of which are seen in these slides. The one we’ll use is the following.

mobius-triangulation

It’s first and second boundary maps are as follows (in code, because latex takes too much time to type out)

mobiusD1 = numpy.array([
   [-1,-1,-1,-1, 0, 0, 0, 0, 0, 0],
   [ 1, 0, 0, 0,-1,-1,-1, 0, 0, 0],
   [ 0, 1, 0, 0, 1, 0, 0,-1,-1, 0],
   [ 0, 0, 0, 1, 0, 0, 1, 0, 1, 1],
])

mobiusD2 = numpy.array([
   [ 1, 0, 0, 0, 1],
   [ 0, 0, 0, 1, 0],
   [-1, 0, 0, 0, 0],
   [ 0, 0, 0,-1,-1],
   [ 0, 1, 0, 0, 0],
   [ 1,-1, 0, 0, 0],
   [ 0, 0, 0, 0, 1],
   [ 0, 1, 1, 0, 0],
   [ 0, 0,-1, 1, 0],
   [ 0, 0, 1, 0, 0],
])

And if we were to run the above code on it we’d get a first Betti number of zero (which is incorrect, it’s first homology group has rank 1). Here are the reduced matrices.

>>> A1, B1 = simultaneousReduce(mobiusD1, mobiusD2)
>>> A1
array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]])
>>> B1
array([[ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  1,  0,  0,  0],
       [ 1, -1,  0,  0,  0],
       [ 0,  0,  0,  0,  1],
       [ 0,  1,  1,  0,  0],
       [ 0,  0, -1,  1,  0],
       [ 0,  0,  1,  0,  0]])

The first reduced matrix looks fine; there’s nothing we can do to improve it. But the second one is not quite fully reduced! Notice that rows 5, 8 and 10 are not linearly independent. So we need to further row-reduce the nonzero part of this matrix before we can read off the true rank in the way we described last time. This isn’t so hard (we just need to reuse the old row-reduce function we’ve been using), but why is this allowed? It’s just because the corresponding column operations for those row operations are operating on columns of all zeros! So we need not worry about screwing up the work we did in column reducing the first matrix, as long as we only work with the nonzero rows of the second.

Of course, nothing is stopping us from ignoring the “corresponding” column operations, since we know we’re already done there. So we just have to finish row reducing this matrix.

This changes our bettiNumber function by adding a single call to a row-reduce function which we name so as to be clear what’s happening. The resulting function is

def bettiNumber(d_k, d_kplus1):
   A, B = numpy.copy(d_k), numpy.copy(d_kplus1)
   simultaneousReduce(A, B)
   finishRowReducing(B)

   dimKChains = A.shape[1]
   kernelDim = dimKChains - numPivotCols(A)
   imageDim = numPivotRows(B)

   return kernelDim - imageDim

And running this on our Mobius band example gives:

>>> bettiNumber(mobiusD1, mobiusD2))
1

As desired. Just to make sure things are going swimmingly under the hood, we can check to see how finishRowReducing does after calling simultaneousReduce

>>> simultaneousReduce(mobiusD1, mobiusD2)
(array([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 1, 0, 0, 0, 0, 0, 0]]), array([[ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  0,  0,  0,  0],
       [ 0,  1,  0,  0,  0],
       [ 1, -1,  0,  0,  0],
       [ 0,  0,  0,  0,  1],
       [ 0,  1,  1,  0,  0],
       [ 0,  0, -1,  1,  0],
       [ 0,  0,  1,  0,  0]]))
>>> finishRowReducing(mobiusD2)
array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0]])

Indeed, finishRowReducing finishes row reducing the second boundary matrix. Note that it doesn’t preserve how the rows of zeros lined up with the pivot columns of the reduced version of \partial_1 as it did in the previous post, but since in the end we’re only counting pivots it doesn’t matter how we switch rows. The “zeros lining up” part is just for a conceptual understanding of how the image lines up with the kernel for a valid simplicial complex.

In fixing this issue we’ve also fixed an issue another commenter mentioned, that you couldn’t blindly plug in the zero matrix for \partial_0 and get zeroth homology (which is the same thing as connected components). After our fix you can.

Of course there still might be bugs, but I have so many drafts lined up on this blog (and research papers to write, experiments to run, theorems to prove), that I’m going to put off writing a full test suite. I’ll just have to update this post with new bug fixes as they come. There’s just so much math and so little time :) But extra kudos to my amazing readers who were diligent enough to run examples and spot my error. I’m truly blessed to have you on my side.

Also note that this isn’t the most efficient way to represent the simplicial complex data, or the most efficient row reduction algorithm. If you’re going to run the code on big inputs, I suggest you take advantage of sparse matrix algorithms for doing this sort of stuff. You can represent the simplices as entries in a dictionary and do all sorts of clever optimizations to make the algorithm effectively linear time in the number of simplices.

Until next time!

How to Conquer Tensorphobia

A professor at Stanford once said,

If you really want to impress your friends and confound your enemies, you can invoke tensor products… People run in terror from the \otimes symbol.

He was explaining some aspects of multidimensional Fourier transforms, but this comment is only half in jest; people get confused by tensor products. It’s often for good reason. People who really understand tensors feel obligated to explain it using abstract language (specifically, universal properties). And the people who explain it in elementary terms don’t really understand tensors.

This post is an attempt to bridge the gap between the elementary and advanced understandings of tensors. We’ll start with the elementary (axiomatic) approach, just to get a good feel for the objects we’re working with and their essential properties. Then we’ll transition to the “universal” mode of thought, with the express purpose of enlightening us as to why the properties are both necessary and natural.

But above all, we intend to be sufficiently friendly so as to not make anybody run in fear. This means lots of examples and preferring words over symbols. Unfortunately, we simply can’t get by without the reader knowing the very basics of linear algebra (the content of our first two primers on linear algebra (1) (2), though the only important part of the second is the definition of an inner product).

So let’s begin.

Tensors as a Bunch of Axioms

Before we get into the thick of things I should clarify some basic terminology. Tensors are just vectors in a special vector space. We’ll see that such a vector space comes about by combining two smaller vector spaces via a tensor product. So the tensor product is an operation combining vector spaces, and tensors are the elements of the resulting vector space.

Now the use of the word product is quite suggestive, and it may lead one to think that a tensor product is similar or related to the usual direct product of vector spaces. In fact they are related (in very precise sense), but they are far from similar. If you were pressed, however, you could start with the direct product of two vector spaces and take a mathematical machete to it until it’s so disfigured that you have to give it a new name (the tensor product).

With that image in mind let’s see how that is done. For the sake of generality we’ll talk about two arbitrary finite-dimensional vector spaces V, W of dimensions n, m. Recall that the direct product  V \times W is the vector space of pairs (v,w) where v comes from V and w from W. Recall that addition in this vector space is defined componentwise ((v_1,w_1) + (v_2, w_2) = (v_1 + v_2, w_1 + w_2)) and scalar multiplication scales both components \lambda (v,w) = (\lambda v, \lambda w).

To get the tensor product space V \otimes W, we make the following modifications. First, we redefine what it means to do scalar multiplication. In this brave new tensor world, scalar multiplication of the whole vector-pair is declared to be the same as scalar multiplication of any component you want. In symbols,

\displaystyle \lambda (v, w) = (\lambda v, w) = (v, \lambda w)

for all choices of scalars \lambda and vectors v, w. Second, we change the addition operation so that it only works if one of the two components are the same. In symbols, we declare that

(v, w) + (v', w) = (v + v', w)

only works because w is the same in both pieces, and with the same rule applying if we switch the positions of v,w above. All other additions are simply declared to be new vectors. I.e. (x,y) + (z,w) is simply itself. It’s a valid addition — we need to be able to add stuff to be a vector space — but you just can’t combine it any further unless you can use the scalar multiplication to factor out some things so that y=w or x=z. To say it still one more time, a general element of the tensor V \otimes W is a sum of these pairs that can or can’t be combined by addition (in general things can’t always be combined).

Finally, we rename the pair (v,w) to v \otimes w, to distinguish it from the old vector space V \times W that we’ve totally butchered and reanimated, and we call the tensor product space as a whole V \otimes W. Those familiar with this kind of abstract algebra will recognize quotient spaces at work here, but we won’t use that language except to note that we cover quotients and free spaces elsewhere on this blog, and that’s the formality we’re ignoring.

As an example, say we’re taking the tensor product of two copies of \mathbb{R}. This means that our space \mathbb{R} \otimes \mathbb{R} is comprised of vectors like 3 \otimes 5, and moreover that the following operations are completely legitimate.

3 \otimes 5 + 1 \otimes (-5) = 3 \otimes 5 + (-1) \otimes 5 = 2 \otimes 5

6 \otimes 1 + 3\pi \otimes \pi = 3 \otimes 2 + 3 \otimes \pi^2 = 3 \otimes (2 + \pi^2)

Cool. This seemingly innocuous change clearly has huge implications on the structure of the space. We’ll get to specifics about how different tensors are from regular products later in this post, but for now we haven’t even proved this thing is a vector space. It might not be obvious, but if you go and do the formalities and write the thing as a quotient of a free vector space (as we mentioned we wouldn’t do) then you know that quotients of vector spaces are again vector spaces. So we get that one for free. But even without that it should be pretty obvious: we’re essentially just declaring that all the axioms of a vector space hold when we want them to. So if you were wondering whether

\lambda (a \otimes b + c \otimes d) = \lambda(a \otimes b) + \lambda(c \otimes d)

The answer is yes, by force of will.

So just to recall, the axioms of a tensor space V \otimes W are

  1. The “basic” vectors are v \otimes w for v \in V, w \in W, and they’re used to build up all other vectors.
  2. Addition is symbolic, unless one of the components is the same in both addends, in which case (v_1, w) + (v_2, w) = (v_1+ v_2, w) and (v, w_1) + (v,w_2) = (v, w_1 + w_2).
  3. You can freely move scalar multiples around the components of v \otimes w.
  4. The rest of the vector space axioms (distributivity, additive inverses, etc) are assumed with extreme prejudice.

Naturally, one can extend this definition to n-fold tensor products, like V_1 \otimes V_2 \otimes \dots \otimes V_d. Here we write the vectors as sums of things like v_1 \otimes \dots \otimes v_d, and we enforce that addition can only be combined if all but one coordinates are the same in the addends, and scalar multiples move around to all coordinates equally freely.

So where does it come from?!

By now we have this definition and we can play with tensors, but any sane mathematically minded person would protest, “What the hell would cause anyone to come up with such a definition? I thought mathematics was supposed to be elegant!”

It’s an understandable position, but let me now try to convince you that tensor products are very natural. The main intrinsic motivation for the rest of this section will be this:

We have all these interesting mathematical objects, but over the years we have discovered that the maps between objects are the truly interesting things.

A fair warning: although we’ll maintain a gradual pace and informal language in what follows, by the end of this section you’ll be reading more or less mature 20th-century mathematics. It’s quite alright to stop with the elementary understanding (and skip to the last section for some cool notes about computing), but we trust that the intrepid readers will push on.

So with that understanding we turn to multilinear maps. Of course, the first substantive thing we study in linear algebra is the notion of a linear map between vector spaces. That is, a map f: V \to W that factors through addition and scalar multiplication (i.e. f(v + v') = f(v) + f(v') and f(\lambda v) = \lambda f(v)).

But it turns out that lots of maps we work with have much stronger properties worth studying. For example, if we think of matrix multiplication as an operation, call it m, then m takes in two matrices and spits out their product

m(A,B) = AB

Now what would be an appropriate notion of linearity for this map? Certainly it is linear in the first coordinate, because if we fix B then

m(A+C, B) = (A+C)B = AB + CB = m(A,B) + m(C,B)

And for the same reason it’s linear in the second coordinate. But it is most definitely not linear in both coordinates simultaneously. In other words,

m(A+B, C+D) = (A+B)(C+D) = AC + AD + BC + BD \neq AC + BD = m(A,C) + m(B,D)

In fact, if the only operation satisfying linearity in its two coordinates separately and also this kind of linearity is the zero map! (Try to prove this as an exercise.) So the strongest kind of linearity we could reasonably impose is that m is linear in each coordinate when all else is fixed. Note that this property allows us to shift around scalar multiples, too. For example,

\displaystyle m(\lambda A, B) = \lambda AB = A (\lambda B) = m(A, \lambda B) = \lambda m(A,B)

Starting to see the wispy strands of a connection to tensors? Good, but hold it in for a bit longer. This single-coordinate-wise-linear property is called bilinearity when we only have two coordinates, and multilinearity when we have more.

Here are some examples of nice multilinear maps that show up everywhere:

  • If V is an inner product space over \mathbb{R}, then the inner product is bilinear.
  • The determinant of a matrix is a multilinear map if we view the columns of the matrix as vector arguments.
  • The cross product of vectors in \mathbb{R}^3 is bilinear.

There are many other examples, but you should have at least passing familiarity with these notions, and it’s enough to convince us that multilinearity is worth studying abstractly.

And so what tensors do is give a sort of classification of multilinear maps. The idea is that every multilinear map f from a product vector space U_1 \times \dots \times U_d to any vector space Y can be written first as a multilinear map to the tensor space

\displaystyle \alpha : U_1 \times \dots \times U_d \to U_1 \otimes \dots \otimes U_d

Followed by a linear map to Y,

\displaystyle \hat{f} : U_1 \otimes \dots \otimes U_d \to Y

And the important part is that \alpha doesn’t depend on the original f (but \hat{f} does). One usually draws this as a single diagram:

comm-diagram-tensor

And to say this diagram commutes is to say that all possible ways to get from one point to another are equivalent (the compositions of the corresponding maps you follow are equal, i.e. f = \hat{f} \alpha).

In fuzzy words, the tensor product is like the gatekeeper of all multilinear maps, and \alpha is the gate. Yet another way to say this is that \alpha is the most general possible multilinear map that can be constructed from U_1 \times \dots \times U_d. Moreover, the tensor product itself is uniquely defined by having a “most-general” \alpha (up to isomorphism). This notion is often referred to by mathematicians as the “universal property” of the tensor product. And they might say something like “the tensor product is initial with respect to multilinear mappings from the standard product.” We discuss language like this in detail in this blog’s series on category theory, but it’s essentially a super-compact (and almost too vague) way of saying what the diagram says.

Let’s explore this definition when we specialize to a tensor of two vector spaces, and it will give us a good understanding of \alpha (which is really incredibly simple, but people like to muck it up with choices of coordinates and summations). So fix V, W as vector spaces and look at the diagram

comm-diagram-tensor-2

What is \alpha in this case? Well it just sends (v,w) \mapsto v \otimes w. Is this map multilinear? Well if we fix w then

\displaystyle \alpha(v_1 + v_2, w) = (v_1 + v_2) \otimes w = v_1 \otimes w + v_2 \otimes w = \alpha(v_1, w) + \alpha (v_2, w)

and

\displaystyle \alpha(\lambda v, w) = (\lambda v) \otimes w = (\lambda) (v \otimes w) = \lambda \alpha(v,w)

And our familiarity with tensors now tells us that the other side holds too. Actually, rather than say this is a result of our “familiarity with tensors,” the truth is that this is how we know that we need to define the properties of tensors as we did. It’s all because we designed tensors to be the gatekeepers of multilinear maps!

So now let’s prove that all maps f : V \times W \to Y can be decomposed into an \alpha part and a \hat{f} part. To do this we need to know what data uniquely defines a multilinear map. For usual linear maps, all we had to do was define the effect of the map on each element of a basis (the rest was uniquely determined by the linearity property). We know what a basis of V \times W is, it’s just the union of the bases of the pieces. Say that V has a basis v_1, \dots, v_n and W has w_1, \dots, w_m, then a basis for the product is just ((v_1, 0), \dots, (v_n,0), (0,w_1), \dots, (0,w_m)).

But multilinear maps are more nuanced, because they have two arguments. In order to say “what they do on a basis” we really need to know how they act on all possible pairs of basis elements. For how else could we determine f(v_1 + v_2, w_1)? If there are n of the v_i‘s and m of the w_i‘s, then there are nm such pairs f(v_i, w_j).

Uncoincidentally, as V \otimes W is a vector space, its basis can also be constructed in terms of the bases of V and W. You simply take all possible tensors v_i \otimes w_j. Since every v \in V, w \in W can be written in terms of their bases, it’s clear than any tensor \sum_{k} a_k \otimes b_k can also be written in terms of the basis tensors v_i \otimes w_j (by simply expanding each a_k, b_k in terms of their respective bases, and getting a larger sum of more basic tensors).

Just to drive this point home, if (e_1, e_2, e_3) is a basis for \mathbb{R}^3, and (g_1, g_2) a basis for \mathbb{R}^2, then the tensor space \mathbb{R}^3 \otimes \mathbb{R}^2 has basis

(e_1 \otimes g_1, e_1 \otimes g_2, e_2 \otimes g_1, e_2 \otimes g_2, e_3 \otimes g_1, e_3 \otimes g_2)

It’s a theorem that finite-dimensional vector spaces of equal dimension are isomorphic, so the length of this basis (6) tells us that \mathbb{R}^3 \otimes \mathbb{R}^2 \cong \mathbb{R}^6.

So fine, back to decomposing f. All we have left to do is use the data given by f (the effect on pairs of basis elements) to define \hat{f} : V \otimes W \to Y. The definition is rather straightforward, as we have already made the suggestive move of showing that the basis for the tensor space (v_i \otimes w_j) and the definition of f(v_i, w_j) are essentially the same.

That is, just take \hat{f}(v_i \otimes w_j) = f(v_i, w_j). Note that this is just defined on the basis elements, and so we extend to all other vectors in the tensor space by imposing linearity (defining \hat{f} to split across sums of tensors as needed). Is this well defined? Well, multilinearity of f forces it to be so. For if we had two equal tensors, say, \lambda v \otimes w = v \otimes \lambda w, then we know that f has to respect their equality, because f(\lambda v_i, w_j) = f(v_i, \lambda w_j), so \hat{f} will take the same value on equal tensors regardless of which representative we pick (where we decide to put the \lambda). The same idea works for sums, so everything checks out, and f(v,w) is equal to \hat{f} \alpha, as desired. Moreover, we didn’t make any choices in constructing \hat{f}. If you retrace our steps in the argument, you’ll see that everything was essentially decided for us once we fixed a choice of a basis (by our wise decisions in defining V \otimes W). Since the construction would be isomorphic if we changed the basis, our choice of \hat{f} is unique.

There is a lot more to say about tensors, and indeed there are some other useful ways to think about tensors that we’ve completely ignored. But this discussion should make it clear why we define tensors the way we do. Hopefully it eliminates most of the mystery in tensors, although there is still a lot of mystery in trying to compute stuff using tensors. So we’ll wrap up this post with a short discussion about that.

Computability and Stuff

It should be clear by now that plain product spaces V \times W and tensor product spaces V \otimes W are extremely different. In fact, they’re only related in that their underlying sets of vectors are built from pairs of vectors in V and W. Avid readers of this blog will also know that operations involving matrices (like row reduction, eigenvalue computations, etc.) are generally efficient, or at least they run in polynomial time so they’re not crazy impractically slow for modest inputs.

On the other hand, it turns out that almost every question you might want to ask about tensors is difficult to answer computationally. As with the definition of the tensor product, this is no mere coincidence. There is something deep going on with tensors, and it has serious implications regarding quantum computing. More on that in a future post, but for now let’s just focus on one hard problem to answer for tensors.

As you know, the most general way to write an element of a tensor space U_1 \otimes \dots \otimes U_d is as a sum of the basic-looking tensors.

\displaystyle \sum_{k} a_{1,k} \otimes a_{2,k} \otimes \dots \otimes a_{d,k}

where the a_{i,k} may be sums of vectors from U_i themselves. But as we saw with our examples over \mathbb{R}, there can be lots of different ways to write a tensor. If you’re lucky, you can write the entire tensor as a one-term sum, that is just a tensor a \otimes b. If you can do this we call the tensor a pure tensor, or a rank 1 tensor. We then have the following natural definition and problem:

Definition: The rank of a tensor x \in U_1 \otimes \dots \otimes U_d is the minimum number of terms in any representation of x as a sum of pure tensors. The one exception is the zero element, which has rank zero by convention.

Problem: Given a tensor x \in k^{n_1} \otimes k^{n_2} \otimes k^{n_3} where k is a field, compute its rank.

Of course this isn’t possible in standard computing models unless you can represent the elements of the field (and hence the elements of the vector space in question) in a computer program. So we restrict k to be either the rational numbers \mathbb{Q} or a finite field \mathbb{F}_{q}.

Even though the problem is simple to state, it was proved in 1990 (a result of Johan Håstad) that tensor rank is hard to compute. Specifically, the theorem is that

Theorem: Computing tensor rank is NP-hard when k = \mathbb{Q} and NP-complete when k is a finite field.

The details are given in Håstad’s paper, but the important work that followed essentially showed that most problems involving tensors are hard to compute (many of them by reduction from computing rank). This is unfortunate, but also displays the power of tensors. In fact, tensors are so powerful that many believe understanding them better will lead to insight in some very important problems, like finding faster matrix multiplication algorithms or proving circuit lower bounds (which is closely related to P vs NP). Finding low-rank tensor approximations is also a key technique in a lot of recent machine learning and data mining algorithms.

With this in mind, the enterprising reader will probably agree that understanding tensors is both valuable and useful. In the future of this blog we’ll hope to see some of these techniques, but at the very least we’ll see the return of tensors when we delve into quantum computing.

Until next time!