10. Extracting Important Relationships in High Dimensions

For one-dimensional data, we wrote This expression is meaningful for vectors, too, because we can add vectors and divide by scalars. We write and call this the mean of the data. Notice that each component of is the mean of that component of the data. There is not an easy analogue of the median, however (how do you order high dimensional data?) and this is a nuisance. Notice that, just as for the one-dimensional mean, we have (i.e. if you subtract the mean from a data set, the resulting data set has zero mean).

Plotting high dimensional data is tricky. If there are relatively few dimensions, you could just choose two (or three) of them and produce a 2D (or 3D) scatterplot. Figure 10.1 shows such a scatterplot, for data that was originally four dimensional. This is the famous Iris dataset, collected by Edgar Anderson in 1936, and made popular amongst statisticians by Ronald Fisher in that year. I found a copy at the UC Irvine repository of datasets that are important in machine learning. You can find the repository at http://​archive.​ics.​uci.​edu/​ml/​index.​html.

Another simple but useful plotting mechanism is the stem plot. This is can be a useful way to plot a few high dimensional data points. One plots each component of the vector as a vertical line, typically with a circle on the end (easier seen than said; look at Fig. 10.2). The dataset I used for this is the wine dataset, again from the UC Irvine machine learning data repository (you can find this dataset at http://​archive.​ics.​uci.​edu/​ml/​datasets/​Wine). For each of three types of wine, the data records the values of 13 different attributes. In the figure, I show the overall mean of the dataset, and also the mean of each type of wine (also known as the class means, or class conditional means). A natural way to compare class means is to plot them on top of one another in a stem plot (Fig. 10.2).

Another strategy that is very useful when there aren’t too many dimensions is to use a scatterplot matrix. To build one, you lay out scatterplots for each pair of variables in a matrix. On the diagonal, you name the variable that is the vertical axis for each plot in the row, and the horizontal axis in the column. This sounds more complicated than it is; look at the example of Fig. 10.3, which shows both a 3D scatter plot and a scatterplot matrix for the same dataset.

Figure 10.4 shows a scatter plot matrix for four of the variables in the height-weight dataset of http://​www2.​stetson.​edu/​~jrasp/​data.​htm (look for bodyfat.xls at that URL). This is originally a 16-dimensional dataset, but a 16 by 16 scatterplot matrix is squashed and hard to interpret. For Fig. 10.4, you can see that weight and adiposity appear to show quite strong correlations, but weight and age are pretty weakly correlated. Height and age seem to have a low correlation. It is also easy to visualize unusual data points. Usually one has an interactive process to do so—you can move a “brush” over the plot to change the color of data points under the brush.

Variance, standard deviation and correlation can each be obtained by performing a more general operation on data. We have a dataset of N vectors x _i , and we are interested in relationships between the j’th and the k’th components. As with correlation, we would like to know whether one component tends to be large (resp. small) when the other is large. Remember that I write x _i ^(j) for the j’th component of the i’th vector.

Just like mean, standard deviation and variance, covariance can refer either to a property of a dataset (as in the definition here) or a particular expectation (as in Chap. 4). From the expression, it should be clear we have seen examples of covariance already. Notice that

which you can prove by substituting the expressions. Recall that variance measures the tendency of a dataset to be different from the mean, so the covariance of a dataset with itself is a measure of its tendency not to be constant. More important is the relationship between covariance and correlation, in the box below.

This is occasionally a useful way to think about correlation. It says that the correlation measures the tendency of {x} and {y} to be larger (resp. smaller) than their means for the same data points, compared to how much they change on their own.

Covariance matrices are often written as , whatever the dataset (you get to figure out precisely which dataset is intended, from context). Generally, when we want to refer to the j, k’th entry of a matrix , we will write , so is the covariance between the j’th and k’th components of the data.

We have a d dimensional dataset {x}. An affine transformation of this data is obtained by choosing some matrix and vector b, then forming a new dataset {m}, where . Here doesn’t have to be square, or symmetric, or anything else; it just has to have second dimension d.

It is easy to compute the mean and covariance of {m}. We have so you get the new mean by multiplying the original mean by and adding b.

The new covariance matrix is easy to compute as well. We have:

All this means that we can try and choose affine transformations that yield “good” means and covariance matrices. It is natural to choose b so that the mean of the new dataset is zero. An appropriate choice of can reveal a lot of information about the dataset.

Recall a matrix is symmetric if . A symmetric matrix is necessarily square. Assume is a d × d symmetric matrix, u is a d × 1 vector, and λ is a scalar. If we have then u is referred to as an eigenvector of and λ is the corresponding eigenvalue. Matrices don’t have to be symmetric to have eigenvectors and eigenvalues, but the symmetric case is the only one of interest to us.

In the case of a symmetric matrix, the eigenvalues are real numbers, and there are d distinct eigenvectors that are normal to one another, and can be scaled to have unit length. They can be stacked into a matrix . This matrix is orthonormal, meaning that .

This means that there is a diagonal matrix and an orthonormal matrix such that In fact, there is a large number of such matrices, because we can reorder the eigenvectors in the matrix , and the equation still holds with a new , obtained by reordering the diagonal elements of the original . There is no reason to keep track of this complexity. Instead, we adopt the convention that the elements of are always ordered so that the elements of are sorted along the diagonal, with the largest value coming first. This gives us a particularly important procedure.

We start with a dataset of N d-dimensional vectors . We can translate this dataset to have zero mean, forming a new dataset {m} where . Now recall that, if we were to form a new dataset {a} where

the covariance matrix of {a} would be Recall also we can diagonalize to get

But this means we could form the dataset {r}, using the rule The mean of this new dataset is clearly 0. The covariance of this dataset is where is a diagonal matrix of eigenvalues of that we obtained by diagonalization. We now have a very useful fact about {r}: its covariance matrix is diagonal. This means that every pair of distinct components has covariance zero, and so has correlation zero. Remember that, in describing diagonalization, we adopted the convention that the eigenvectors of the matrix being diagonalized were ordered so that the eigenvalues are sorted in descending order along the diagonal of . Our choice of ordering means that the first component of r has the highest variance, the second component has the second highest variance, and so on.

The transformation from {x} to {r} is a translation followed by a rotation (remember is orthonormal, and so a rotation). So this transformation is a high dimensional version of what I showed in Figs. 10.5 and 10.6.

The mean of {r} is zero, and the covariance matrix of {r} is diagonal. It is quite usual for high dimensional datasets to have a small number of large values on the diagonal, and a lot of small values. These values give the variance of the corresponding components of {r}. Now imagine choosing a component of {r} that has a small variance, and replacing that with zero. Because this component has zero mean (like every other one), and small variance, replacing it with zero will not result in much error.

If we can replace components with zero without causing much error, the blob of data is really a low dimensional blob in a high dimensional space. For example, think about a blob lying along, but very close to the x-axis in 3D. Replacing each data items y and z values with zero will not change the shape of the blob very much. As a visual example, look at Fig. 10.3; the scatterplot matrix strongly suggests that the blob of data is flattened (eg look at the petal width vs petal length plot).

The data set {r} is d-dimensional. We will try to represent it with an s dimensional dataset, and see what error we incur. Choose some s < d. Now take each data point r _i and replace the last d − s components with 0. Call the resulting data item p _i . We should like to know the average error in representing r _i with p _i .

This error is Write r _i ^(j) for the j′ component of r _i , and so on. Remember that p _i is zero in the last d − s components. The error is then Because {r} has zero mean, we have that is the variance of the j’th component of {r}. So the error is which is the sum of the diagonal elements of the covariance matrix from s + 1, s + 1 to d, d. If this sum is small compared to the sum of the first s components, then dropping the last d − s components results in a small error. In that case, we could think about the data as being s dimensional. Figure 10.7 shows the result of using this approach to represent the blob we’ve used as a running example as a 1D dataset.

This is an observation of great practical importance. As a matter of experimental fact, a great deal of high dimensional data produces relatively low dimensional blobs. We can identify the main directions of variation in these blobs, and use them to understand and to represent the dataset.

Translating a blob of data doesn’t change the scatterplot matrix in any interesting way (the axes change, but the picture doesn’t). Rotating a blob produces really interesting results, however. Figure 10.8 shows the dataset of Fig. 10.4, translated to the origin and rotated to diagonalize it. Now we do not have names for each component of the data (they’re linear combinations of the original components), but each pair is now not correlated. This blob has some interesting shape features. Figure 10.8 shows the gross shape of the blob best. Each panel of this figure has the same scale in each direction. You can see the blob extends about 80 units in direction 1, but only about 15 units in direction 2, and much less in the other two directions. You should think of this blob as being rather cigar-shaped; it’s long in one direction, but there isn’t much in the others. The cigar metaphor isn’t perfect because there aren’t any four dimensional cigars, but it’s helpful. You can think of each panel of this figure as showing views down each of the four axes of the cigar.

Now look at Fig. 10.9. This shows the same rotation of the same blob of data, but now the scales on the axis have changed to get the best look at the detailed shape of the blob. First, you can see that blob is a little curved (look at the projection onto direction 2 and direction 4). There might be some effect here worth studying. Second, you can see that some points seem to lie away from the main blob. I have plotted each data point with a dot, and the interesting points with a number. These points are clearly special in some way.

The problem with these figures is that the axes are meaningless. The components are weighted combinations of components of the original data, so they don’t have any units, etc. This is annoying, and often inconvenient. But I obtained Fig. 10.8 by translating, rotating and projecting data. It’s straightforward to undo the rotation and the translation—this takes the projected blob (which we know to be a good approximation of the rotated and translated blob) back to where the original blob was. Rotation and translation don’t change distances, so the result is a good approximation of the original blob, but now in the original blob’s coordinates. Figure 10.10 shows what happens to the data of Fig. 10.4. This is a two dimensional version of the original dataset, embedded like a thin pancake of data in a four dimensional space. Crucially, it represents the original dataset quite accurately.

We start with a dataset of N d-dimensional vectors . We translate this dataset to have zero mean, forming a new dataset {m} where . We diagonalize to get and form the dataset {r}, using the rule We saw the mean of this dataset is zero, and the covariance is diagonal. We then represented the d-dimensional data set {r} with an s dimensional dataset, by choosing some r < d, then taking each data point r _i and replacing the last d − s components with 0. We call the resulting data item p _i .

Now consider undoing the rotation and translation. We would form a new dataset , with the i’th element given by (you should check this expression). But this expression says that is constructed by forming a weighted sum of the first s columns of (because all the other components of p _i are zero), then adding . If we write u _j for the j’th column of , we have What is important about this sum is that s is usually a lot less than d. The u _j are known as principal components of the dataset. You can easily derive an expression for r _i ^(j) from all this (exercises), but for reference, here it is:

We can easily determine the error in approximating with . The error in representing {r} by {p _s } was easy to compute. We had This was the sum of the diagonal elements of the covariance matrix of {r} from s, s to d, d. If this sum is small compared to the sum of the first s components, then dropping the last d − s components results in a small error.

The error in representing {x} with is now easy to get. Rotations and translations do not change lengths. This means that which is the sum of the diagonal elements of the covariance matrix of {r} from s, s to d, d which is easy to evaluate, because these are the values of the d − s eigenvalues that we decided to ignore. Now we could choose s by identifying how much error we can tolerate. More usual is to plot the eigenvalues of the covariance matrix, and look for a “knee”, like that in Fig. 10.11. You can see that the sum of remaining eigenvalues is small.

Diffuse surfaces reflect light uniformly in all directions. Examples of diffuse surfaces include matte paint, many styles of cloth, many rough materials (bark, cement, stone, etc.). One way to tell a diffuse surface is that it does not look brighter (or darker) when you look at it along different directions. Diffuse surfaces can be colored, because the surface reflects different fractions of the light falling on it at different wavelengths. This effect can be represented by measuring the spectral reflectance of a surface, which is the fraction of light the surface reflects as a function of wavelength. This is usually measured in the visual range of wavelengths (about 380 nm to about 770 nm). Typical measurements are every few nm, depending on the measurement device. I obtained data for 1995 different surfaces from http://​www.​cs.​sfu.​ca/​~colour/​data/​ (there are a variety of great datasets here, from Kobus Barnard).

Each spectrum has 101 measurements, which are spaced 4 nm apart. This represents surface properties to far greater precision than is really useful. Physical properties of surfaces suggest that the reflectance can’t change too fast from wavelength to wavelength. It turns out that very few principal components are sufficient to describe almost any spectral reflectance function. Figure 10.11 shows the mean spectral reflectance of this dataset, and Fig. 10.11 shows the eigenvalues of the covariance matrix.

This is tremendously useful in practice. One should think of a spectral reflectance as a function, usually written ρ(λ). What the principal components analysis tells us is that we can represent this function rather accurately on a (really small) finite dimensional basis. This basis is shown in Fig. 10.11. This means that there is a mean function r(λ) and k functions ϕ _m (λ) such that, for any ρ(λ), where e(λ) is the error of the representation, which we know is small (because it consists of all the other principal components, which have tiny variance). In the case of spectral reflectances, using a value of k around 3–5 works fine for most applications (Fig. 10.12). This is useful, because when we want to predict what a particular object will look like under a particular light, we don’t need to use a detailed spectral reflectance model; instead, it’s enough to know the c _i for that object. This comes in useful in a variety of rendering applications in computer graphics. It is also the key step in an important computer vision problem, called color constancy. In this problem, we see a picture of a world of colored objects under unknown colored lights, and must determine what color the objects are. Modern color constancy systems are quite accurate, even though the problem sounds underconstrained. This is because they are able to exploit the fact that relatively few c _i are enough to accurately describe a surface reflectance.

An image is usually represented as an array of values. We will consider intensity images, so there is a single intensity value in each cell. You can turn the image into a vector by rearranging it, for example stacking the columns onto one another. This means you can take the principal components of a set of images. Doing so was something of a fashionable pastime in computer vision for a while, though there are some reasons that this is not a great representation of pictures. However, the representation yields pictures that can give great intuition into a dataset.

Figure 10.14 shows the mean of a set of face images encoding facial expressions of Japanese women (available at http://​www.​kasrl.​org/​jaffe.​html; there are tons of face datasets at http://​www.​face-rec.​org/​databases/​). I reduced the images to 64 × 64, which gives a 4096 dimensional vector. The eigenvalues of the covariance of this dataset are shown in Fig. 10.13; there are 4096 of them, so it’s hard to see a trend, but the zoomed figure suggests that the first couple of dozen contain most of the variance. Once we have constructed the principal components, they can be rearranged into images; these images are shown in Fig. 10.14. Principal components give quite good approximations to real images (Fig. 10.15).

The principal components sketch out the main kinds of variation in facial expression. Notice how the mean face in Fig. 10.14 looks like a relaxed face, but with fuzzy boundaries. This is because the faces can’t be precisely aligned, because each face has a slightly different shape. The way to interpret the components is to remember one adjusts the mean towards a data point by adding (or subtracting) some scale times the component. So the first few principal components have to do with the shape of the haircut; by the fourth, we are dealing with taller/shorter faces; then several components have to do with the height of the eyebrows, the shape of the chin, and the position of the mouth; and so on. These are all images of women who are not wearing spectacles. In face pictures taken from a wider set of models, moustaches, beards and spectacles all typically appear in the first couple of dozen principal components.

We will plot the high dimensional point x _i at v _i , which is a two-dimensional vector. Now the squared distance between points i and j in the high dimensional space is (where the superscript is to remind you that this is a squared distance). We could build an N × N matrix of squared distances, which we write . The i, j’th entry in this matrix is D _ij ⁽²⁾(x), and the x argument means that the distances are between points in the high-dimensional space. Now we could choose the v _i to make as small as possible. Doing so should mean that points that are far apart in the high dimensional space are far apart in the plot, and that points that are close in the high dimensional space are close in the plot.

In its current form, the expression is difficult to deal with, but we can refine it. Because translation does not change the distances between points, it cannot change either of the matrices. So it is enough to solve the case when the mean of the points x _i is zero. We can assume that Now write 1 for the n-dimensional vector containing all ones, and for the identity matrix. Notice that Now write Using this expression, you can show that the matrix , defined below, has i, jth entry x _i ⋅ x _j (exercises). I now argue that, to make is close to , it is enough to make close to . Proving this will take us out of our way unnecessarily, so I omit a proof.

We need some notation. Take the dataset of N d-dimensional column vectors x _i , and form a matrix by stacking the vectors, so In this notation, we have Notice is symmetric, and it is positive semidefinite. It can’t be positive definite, because the data is zero mean, so .

We must now choose a set of v _i that makes close to . We do so by choosing a that is close to . But this means we must choose so that is close to . We are computing an approximate factorization of the matrix .

We seek a set of k dimensional v that can be stacked into a matrix . This must produce a that must (a) be as close as possible to and (b) have rank at most k. It can’t have rank larger than k because there must be some which is N × k so that . The rows of this are our v _i ^T .

We can obtain the best factorization of from a diagonalization. Write write for the matrix of eigenvectors of and for the diagonal matrix of eigenvalues sorted in descending order, so we have and write for the matrix of positive square roots of the eigenvalues. Now we have which allows us to write

Now think about approximating by the matrix . The error is a sum of squares of the entries, Because is a rotation, it is straightforward to show that But which means that we could find from the best rank k approximation to . This is obtained by setting all but the k largest entries of to zero. Call the resulting matrix . Then we have and The first k columns of are non-zero. We drop the remaining N − k columns of zeros. The rows of the resulting matrix are our v _i , and we can plot these. This method for constructing a plot is known as principal coordinate analysis.

This plot might not be perfect, because reducing the dimension of the data points should cause some distortions. In many cases, the distortions are tolerable. In other cases, we might need to use a more sophisticated scoring system that penalizes some kinds of distortion more strongly than others. There are many ways to do this; the general problem is known as multidimensional scaling.

Multidimensional scaling gets positions (the of Sect. 10.4.1) from distances (the of Sect. 10.4.1). This means we can use the method to build maps from distances alone. I collected distance information from the web (I used http://​www.​distancefromto.​net, but a google search on “city distances” yields a wide range of possible sources), then applied multidimensional scaling. I obtained distances between the South African provincial capitals, in kilometers. I then used principal coordinate analysis to find positions for each capital, and rotated, translated and scaled the resulting plot to check it against a real map (Fig. 10.16).

One natural use of principal coordinate analysis is to see if one can spot any structure in a dataset. Does the dataset form a blob, or is it clumpy? This isn’t a perfect test, but it’s a good way to look and see if anything interesting is happening. In Fig. 10.17, I show a 3D plot of the spectral data, reduced to three dimensions using principal coordinate analysis. The plot is quite interesting. You should notice that the data points are spread out in 3D, but actually seem to lie on a complicated curved surface—they very clearly don’t form a uniform blob. To me, the structure looks somewhat like a butterfly. I don’t know why this occurs (perhaps the universe is doodling), but it certainly suggests that something worth investigating is going on. Perhaps the choice of samples that were measured is funny; perhaps the measuring instrument doesn’t make certain kinds of measurement; or perhaps there are physical processes that prevent the data from spreading out over the space.

Our algorithm has one really interesting property. In some cases, we do not actually know the datapoints as vectors. Instead, we just know distances between the datapoints. This happens often in the social sciences, but there are important cases in computer science as well. As a rather contrived example, one could survey people about breakfast foods (say, eggs, bacon, cereal, oatmeal, pancakes, toast, muffins, kippers and sausages for a total of 9 items). We ask each person to rate the similarity of each pair of distinct items on some scale. We advise people that similar items are ones where, if they were offered both, they would have no particular preference; but, for dissimilar items, they would have a strong preference for one over the other. The scale might be “very similar”, “quite similar”, “similar”, “quite dissimilar”, and “very dissimilar” (scales like this are often called Likert scales). We collect these similarities from many people for each pair of distinct items, and then average the similarity over all respondents. We compute distances from the similarities in a way that makes very similar items close and very dissimilar items distant. Now we have a table of distances between items, and can compute a and produce a scatter plot. This plot is quite revealing, because items that most people think are easily substituted appear close together, and items that are hard to substitute are far apart. The neat trick here is that we did not start with a , but with just a set of distances; but we were able to associate a vector with “eggs”, and produce a meaningful plot.