Principal Components Analysis in 2D

Miguel A. Lerma

November 11, 2019

Abstract

Here we study 2-dimensional PCA and discuss an application to

the location of a set of points in the plane in an elliptical region.

1 Introduction

Principal component analysis (PCA) is a mathematical procedure intended

to replace a number of correlated variables with a new set of variables that
are linearly uncorrelated. This can be achieved with a change of variables
given by an orthogonal transformation. The technique is well known, how-
ever the case of two variables is particularly simple to analyze, and has
some interesting application to image analysis. Here we analyze the math
behind the method, which becomes particularly simple to approach in two
dimensions, and use it to solve the problem of finding a elliptical region
approximately encompassing a given set of points. In our example we will
use the method to image (1) consisting of a rotated “HELLO” word.

2 The Math behind PCA

Assume that we have a collection of 2-dimensional points (x1 , y1 ), (x2 , y2 ),. . . ,

(xn , yn ). The correlation between the xi and the
P yi can be measured by their
variance σxy = (x − x̄)(y − ȳ), where u = n1 ni=1 ui represents the average
value of u. In order to simplify computations we will assume x̄ = ȳ = 0,
which can be accomplished with a translation (x, y) 7→ (x − x̄, y − ȳ).

1
 Figure 1: A rotated “HELLO”

Our purpose is to find a change of coordinates (x, y) 7→ (x0 , y 0 ) given by a

rotation:
 0     
x x cos θ − sin θ
=R , R= ,
y0 y sin θ cos θ

such that x0 and y 0 are uncorrelated, i.e., σx0 y0 = 0.

Let X be the matrix whose columns are the original coordinates of the given
points:  
x1 x2 · · · xn
X=
y1 y2 · · · yn
and X 0 is the analogous matrix with the new coordinates:
 0
x1 x02 · · · x0n

0
X = = RX .
y10 y20 · · · yn0
The correlation matrix of X is
σx2 σxy
 
1
C = XX T = ,
n σxy σy2

and that of X 0 is
1 0 0T 1
C0 = X X = RXX T RT = RCRT ,
n n
Since R is a rotation it is orthogonal and we have RT = R−1 . On the other
hand we want σx0 y0 = 0, hence C 0 must be a diagonal matrix:
 2 
σx0 0
C0 = .
0 σy20

2
So the problem amounts to diagonalizing C. This can be done in the usual
way, by finding the eigenvalues and eigenvectors of C.

The eigenvalues of C are the roots of its characteristic polynomial:

σx2 − λ σxy
= λ2 − (σx2 + σy2 )λ + σx2 σy2 − σxy
2
σxy σy2 − λ

The roots can be found with the usual formula to solve second degree equa-
tions: q
σx2 + σy2 ± (σx2 + σy2 )2 − 4σx2 σy2 + 4σxy
2
λ=
q 2
2 2
σx + σy ± (σx − σy2 )2 + 4σxy
2 2
= .
2
Let λ+ be the root obtained using the plus sign in the formula, and λ− the
one obtained with the minus sign. Then it can be easily checked that the
following are eigenvectors verifying CV± = λ± V± :
 2   2 
σx + σxy − λ− σx + σxy − λ+
V+ = V− =
σy2 + σxy − λ− σy2 + σxy − λ+

These vectors can now be normalized U+ = k+ V+ , U− = k− V− , with k+

and k+ chosen so that ||U+ || = ||U− || = 1, and U+ × U− = 1. If we choose
the rotation matrix R so that
 
T U+,x U−,x
R = (U+ |U− ) = ,
U+,y U−,y

then λ+ = σx20 , and λ− = σy20 . Since λ+ ≥ λ− we will have that x0 will be

the variable with the maximum dispersion, and y 0 the one with minimum
dispersion.

3 Application to optimal location of a set of points

in an elliptical region

We now apply the method to locating the position, orientation, and ap-
proximate size of a set of points {(x1 , y1 ), (x2 , y2 ), . . . , (xn , yn )} in the plane,
and illustrate it for the example shown in figure (1). The set of points in

3
this example is the one given by the black points of the word “HELLO”
in the image. Each has some coordinates (xi , yi ), and its barycenter (x̄, ȳ)
is the first reference we use for the location of the set on the plane. The
eigenvectors V+ and V− determine the directions of maximum and minimum
dispersion, and their directions can be used for the axes of an ellipse centered
in the barycenter. It remains the determine the values of the semi axes p of
the ellipse, pwhich will be proportional to the standard deviations σx0 = λ+
and σy0 = λ− respectively. In this particular example the sizes has been
chosen to be double of the standard deviations, and the result can be seen
in figure (2). The ellipse has been drawn in green, and the eigenvectors,
scaled so that their lengths coincide with double the standard deviations of
x0 and y 0 , are drawn in red and yellow respectively.

Figure 2: The rotated “HELLO” and its associated ellipse

See the appendix for specific code designed to perform the computations.
The program used in the example shown here was actually implemented in
Matlab, and has not been included here, but the computations are basically
identical to the ones performed by the code included below.

4 Conclusion

We have reviewed the mathematical details of PCA in two dimensions, and

used it to solve a problem in image analysis, namely locating a set of points
in the plane in an elliptical region.

4
Appendix: sample code implementation

The following is C computer code illustrating how to perform the computa-

tions (based on tested code but modified and simplified to be included here,
it may contain typos or bugs):

// some variable declarations

long int sumx=0, sumy=0, sumxx=0, sumyy=0, sumxy=0;
double xbar, ybar, varx, vary, covxy, sumvars, diffvars;
double discriminant, sqrtdiscr, lambdaplus, lambdaminus;
double aplus, bplus, aminus, bminus;
int npix=0;

// scan the image (nLength and nHeight are the dimensions of the image)
for(x=0; x<nLength; x++)
{
for(y=0; y<nHeight; y++)
{
if (region(x,y)) // true if we are in the region of interest
{
sumx+=x; sumy+=y; sumxx+=x*x;
sumyy+=y*y; sumxy+=x*y; npix++;
}
}
}

// baricenter
xbar = ((double)sumx)/npix;
ybar = ((double)sumy)/npix;

// variances and covariance

varx = sumxx/npix - xbar*xbar;
vary = sumyy/npix - ybar*ybar;
covxy = sumxy/npix - xbar*ybar;

sumvars = varx + vary;

diffvars = varx - vary;
discriminant = diffvars*diffvars + 4*covxy*covxy;
sqrtdiscr = sqrt(discriminant);

5
// eigenvalues
lambdaplus = (sumvars + sqrtdiscr)/2;
lambdaminus = (sumvars - sqrtdiscr)/2;

//eigenvectors - these are the components of the two vectors

aplus = varx + covxy - lambdaminus;
bplus = vary + covxy - lambdaminus;

aminus = varx + covxy - lambdaplus;

bminus = vary + covxy - lambdaplus;

// normalizing the vectors

double aParallel; double bParallel;
double aNormal; double bNormal;

double denomPlus = sqrtf(aplus*aplus + bplus*bplus);

double denomMinus= sqrtf(aminus*aminus + bminus*bminus);

aParallel = aplus/denomPlus;
bParallel = bplus/denomPlus;
aNormal = aminus/denomMinus;
bNormal = bminus/denomMinus;

// semi axes
double k = 2 // scale factor
double majoraxis = k*sqrtf(lambdaplus);
double minoraxis = k*sqrtf(lambdaminus);