Subject: Principal components summary

From: "[email-masked]" <[email-masked]>
Date: Mon, 21 Aug 1995 11:10:28 +1

Dear Netters,

Some time ago I had a question about principal components and
plane fitting to points. My original question was:

I am looking for algorithm/numerical methods for the
determining of the principal components of a data set. In the
words of geometry, I have a number of points determined by
y1, y2 and y3, and I am looking for the best fitting plane
to the points, i.e. the plane determined by the two greatest
principal components.

Thanks for everybody who send me the answers. Here is the
summary

Tamas E. Gunda
------------------------------------------------------------
***************************************************

That looks like a rather straightforward least-squares
problem. You want the parameters a,b,c in the equation for
a plane y3 = a y1 + b y2 + c. So you make up a matrix X:

           y1(1)   y2(1)   1  
           y1(2)   y2(2)   1
            ...
           y1(n)   y2(n)   1

and a vector Z = (y3(1), y3(2), ..., y3(n)), plus a
3-component vector of unknowns A = (a,b,c). Then

      A = X^+ Z

gives you the best values for the parameters, where X^+ is
the generalized inverse of X.

Konrad Hinsen   | E-Mail: [email-masked]
Departement de chimie    | Tel.: +1-514-343-6111 ext. 3953
Universite de Montreal   | Fax:  +1-514-343-7586
C.P. 6128, succ. A 
Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7  | Francais (phase experimentale)

*********************************************

There are listings for FORTRAN programs pertaining to all
aspects of factor analyses, including pca in

"Factor Analysis of Data Matrices", P. Horst, Holt, Rinehart
and
Winston, New York, 1965.

I also have written my own code, using the routines TRED2 and
TQLI from "Numerical Recipes in Fortran" to diagonalize the
covariance matrix and locate the eigenvalues and
eigenvectors.

Heather Gordon
Department of Chemistry
Queen's University
Kingston, Ontario

********************************************

        What you say you want sounds a lot like the
transformation of a molecule to the coordinate system in
which the moment-of-inertia tensor is diagonal. If your (y1,
y2, y3) components are *orthogonal*, you ought to be able a
similar approach. This kind of code is buried in all sorts
of programs
.
It should be a simple matter to adapt a suitable code
fragment to you use.

                           FREDERIC A. VAN-CATLEDGE

Scientific Computing Division         ||   Office: (302)
695-1187 or 529-2076
Central Research & Development Dept.  ||          
The DuPont Company                    ||      FAX: (302)
695-9658
P. O. Box 80320                       ||
Wilmington DE 19880-0320              || Internet:
[email-masked]

********************************************

Tamas, 
     A friend sent me your e-mail on the subject of an
algorithm for calculating the first two PC's from a dataset.
I have written PCA (NIPALS) routines in BASIC and C and will
be happy to send them to you, if required. I also have a book
(Numerical Recipes) which details PCA (Singular Value
Decomposition) which does a similar thing. If you still need
some help, send me an e-mail

Steve Gurden,
Bristol Chemometrics Group,
School of Chemistry,
University of Bristol,
Cantock's Close,
BRISTOL BS8 1TS

Tel: +44 (0)117 9289000 extn. 4421
e-mail: [email-masked]

********************************************

  Dear Dr. Gunda,

  One of the most simple and elegant method to perform
Principal Component Analysis is that of the NIPALS algorithm.
It can be formulated in a few lines of programming and it
works really fine.

  You can find it in the following references:

* Analytica Chimica Acta, vol. 185, p. 1-17 (1986)
* Chemometrics and Intelligent Laboratory Systems, vol. 2,
p. 37-52 (1987)

  I hope this helps you.

  Sincerely.

  Jose I. Garcia
--

Dr. Jose Ignacio Garcia-Laureiro              Phone :
34-(9)76-350475
Departamento de Quimica Organica              Fax   :
34-(9)76-567920
Instituto de Ciencia de Materiales de Aragon  e-mail:
[email-masked]
C.S.I.C.-Universidad de Zaragoza                     
[email-masked]
E-50009 ZARAGOZA (SPAIN)

*****************************************************

I thought the principle (no pun intended) was quite simple;
just calculate the variance/covariance matrix of the
variables (or better the correlation matrix) and diagonalize
that. You can get the correlation matrix by autoscaling the
data and calculating the covariance matrix then.

Next, select the number of principal components that you want
by calculating the cumulative sum of the eigenvalues divided
by the sum of all the eigenvaeigenvectors need to be sorted by their size by the
diagonalization procedure or afterwards)

I mean:

T = \sum_i  \labda_i

suppose, as an example:

\labda_1 / T              = 90 %
(\labda_1 + \labda_2) / T = 98 %
 ...
(\labda_1 + \labda_2 + ...) /T = ..

if this is sufficiently accurate for you, use the subspace
spanned by the first two eigenvectors as your new space, and
project all datapoints onto it. otherwise, use more
eigenvectors until (\sum_j \labda_j) / T is sufficiently near
100 % (1.00), and project the datapoints on the subspace
(unfortunately it is no longer a 2-D plane, which would be
most easy for viewing) 

Hopefully you can use this information, I can't think of any
references off-hand..
There is a book by Box (and Hunter?) which describes these
multivariate techniques in great detail.
Sorry I can't help you better,

Frits

Frits Daalmans
OIO Conformational Analysis
Gorlaeus Laboratoria
Leiden, The Netherlands
E-mail: [email-masked]
Tel: [+31] (0)71-274505

********************************************

Check the NIST CD for handwritten characters -- there is a
principal  component algorithm coded there (may be in C, but
you could then recode it in FORTRAN).  There are several
sites that catalog algorithms, too --   here are a few:
http://gams.nist.gov/,   
http://netlib.att.com/netlib/att/cs/home/1127.ht,
http://www.netlib.org/.

Good luck.

Joe McDaniel
[email-masked]

********************************************

Dear Dr. Gunda,
There are many references to Principal Component Analysis.It
is actually a straightforward problem in matrix
diagonalization. However, if your real goal is to simply fit
a plane to a set of points (minimizing the sum of squares of
the distances of the points to the plane) or a line to a set
of points (minimizing the sum of squares of the distances of
the points to the line),
then this is also a simple matrix diagonalization problem
whose solution is given and discussed in the following
references:

"To Fit a Plane or Line to a Set of Points by Least Squares",
V. Schomaker, J. Waser, R.E. Marsh, and G. Bergman, Acta
Cryst.,
vol. 12, pp. 600-604 (1959).

and

"To fit a plane to a set of points by least squares",
D.M. Blow, Acta Cryst., vol. 13, p. 168 (1960).

Regards,

Marvin Waldman, Ph.D.
Director, Rational Drug Design
Biosym Technologies, Inc.
e-mail: [email-masked]

********************************************

Dear Dr. Gunda,

    There are two ways to do what you want (as I understand
it). One is to use a statistical analysis package which
supports principal components analysis as statisticians
define it.  The
plane you are looking for is the one defined by the first two
principal components.

    The other is to perform the equivalent of a moment of
inertia calculation with all masses set to 1.  First find the
mean value of each of the three coordinates (call them x, y,
z).  This is equivlant to your center of mass.  Build the
symmetric 3x3 tensor with the following elements, summed over
all points (distances are relative to the centroid):

    y*y + z*z         -x*y           -x*z
      -x*y          x*x + z*z        -y*z
      -x*z            -y*z         x*x + y*y
    
    Diagonalize the tensor (the Jacobi method works well here
- see any of the "Numerical Recipes" series by Press et all.) 
The transformation matrix gives you your three principal
coordinates.

Hope this helps,
Paul Soper

Paul Soper                        All the usual disclaimers
apply
DuPont Central Research            
[email-masked]  
P.O. Box 80328                                 Tel
(302)-695-1757  
Wilmington, DE 19880-0328                      FAX
(302)-695-8412  

********************************************

   According to the math of the normal least square
procedure, as far as I know, NOT the distances of the points
and the line are minimized (i.e. a perdendicular from a point
to the line), but the  difference of the measured and the
calculated y values,which is represented graphically by a
line between the point and the regression line and it is
parallel with the y axis:

That depends on the error criterion that you use. The one
sent yesterday (and probably some others too) does indeed
minimize the error along one coordinate axis (and in fact
works only if that axis does not lie in the plane that you
want to fit). If you want to minimize the distances of the
points from the plane, use the normal form for the plane:

    n x - d = 0,

where n is the normalized normal vector, x is the coordinate
vector, and d the distance of the origin from the plane.  The
distance of any point y from this plane is simply given by
n y - d, so you must minimize

  __
  \    |           | 2
  /    | n y_i - d |
  --
  i

with the constraint that n is normalized, so you need a
least-squaresminimizer that can handle constraints (e.g. via
Lagrange multipliers). Or use a "dirty hack": if you know (or
assume) that d will never be zero (i.e. the plane will not
contain the coordinate origin), set d = 1 and use an
unnormalized normal vector, whose length then is the inverse
of the distance from the origin.

Konrad Hinsen                     | E-Mail:
[email-masked]
Departement de chimie             | Tel.: +1-514-343-6111
ext. 3953
Universite de Montreal            | Fax:  +1-514-343-7586
C.P. 6128, succ. A                |
Deutsch/Esperanto/English/Nederlands/
Montreal (QC) H3C 3J7             | Francais (phase
experimentale)

********************************************

        What you say you want sounds a lot like the
transformation of a molecule to the coordinate system in
which the moment-of-inertia tensor is diagonal.  If your (y1,
y2, y3) components are *orthogonal*, you ought to be able a
similar approach.  This kind of code is buried in all sorts
of programs.
It should be a simple matter to adapt a suitable code
fragment to you use.

                           FREDERIC A. VAN-CATLEDGE

Scientific Computing Division         ||   Office: (302)
695-1187 or 529-2076
Central Research & Development Dept.  ||          
The DuPont Company                    ||      FAX: (302)
695-9658
P. O. Box 80320                       ||
Wilmington DE 19880-0320              || Internet:
[email-masked] 

********************************************

Tamas,

        I think I have solved a problem fairly similar to the
one you are interested in.  My problem was less general and
involved finding thebest-fit plane to a group of four points
centred around a fifth point and arose because I am looking
at compression strain in tetracoordinate carbon systems.  For
example the molecule [4.4.4.4]fenestrane C9H12      
                               H
                    H2C__C__CH2
                      |  |  |
                     HC-----CH
                      |  |  |
                    H2C--C--CH2
                         H
contains a central "flattened" tetracoordinate carbon atom. 
One way to determine the flattening at any C(C)4 moiety is
to find the best-fit plane for the four alpha-C atoms passing
through the central C atom.  This can be done simply by
finding the eigenvectors of what we have termed the Geometry
Tensor.
        The Geometry Tensor is constructed by multiplying the
matrix D (the vectors to your points -- in this case the 4
alpha-C atoms) by its transpose to give a 3x3 Real Symmetric
Matrix (routines to find the eigenvectors/values of any RSM
can be found in eispak and similar sets of routines).  The
eigenvectors of this matrix will correspond to the minimum
(Best-Fit) -- this will have the smallest eigenvalue, maximum
and one other extremum (i.e. they will form a set of
cartesian axes -- x,y and z -- one of which will be the
normal to your best-fit plane. 
     The mathematics behind this is actually quite simple and
involes setting up the equations to minimise d' = (sum of
squared distances of your points to an orbitrary plane). 
After which the need to form and find the eigenvectors of the
"Geometry Tensor" becomes obvious. 
     Naturally this can be readily extended to any number of
atoms and any origin.  I have used this method to re-orient
molecules that were optimised in C1 symmetry (in cartesians)
which exhibit higher symmetry but where the axes/planes of
symmetry do not lie on the x,y or z axes of my final
cartesians.
        I might be able to send you my fortran code (it is
very simple) if you think this would help but I will need to
check the origin of the eigenvector solving routines that I
use (to make sure we do not break anyone's intellectual
property rights as I did not write these routines myself --
I think they came from eispak).

        Cheers,

                Danne

 Danne R Rasmussen, PhD Student                     phone: 
 +61 6 249-3771
 Research School of Chemistry                         fax: 
 +61 6 249-0750
 Australian National University
 Canberra ACT 0200                             e-mail:
[email-masked]

********************************************

Hi,

     Your assumption about the fitting of the least squares
of perpendicular distances is correct .. this is what is done
by linear regression.  A side note, linear regression is used
in 
teaching software, but programs are always built around the
matrix formulation.

     The way that you can fit perpendicular distances, as
well as do non-linear fits (constants other than linear
coeficients) is using optimization methods, such as steepest
descent, Newton-Rhapson, etc. This is just like doing
geometry optimization, only you are minimizing on your
perpendicular distances instead of energy.

     The other possibility is spline methods (cube splines
are most commonly used).  A spline method is not a least
square fit, it is an exact fit.  However, if there is noise
in the data you fit exactly to the noise.

     Hope this helps.

                    Dave Young
                    [email-masked]

********************************************

From: [email-masked]
 Hi Dr. Tamas Gunda,

    I too have been interested in the question of fitting a
plane to a set of points(3D). As you correctly pointed out
this is a problem in total least squares, not simple least
squares. This problem as been approached many times using a
variety of algorithims(its' fundamental to crystallography).

I have written a simple FORTRAN  program based on the method
of D. M. Blow, Acta. Cryst.(1960),13,168. One fundamental
paper is V. Schomaker, J. Waser,  R. E. Marsh and ?. Bergman,
Acta Cryst(1959),12,60. There are many other papers on this
subject. I hope this helps and good luck.

                                   Brian Schmitz

                                   Univ. of Rhode Island

***************************************************
summary end here

*****************************************************************************
   Tamas E. Gunda, Ph.D.               phone: (+36-52) 316666 ext 2479
   Research Group of Antibiotics       fax  : (+36-52) 310936
   L. Kossuth University               e-mail: [email-masked]
   POBox 36                                   
   H-4010 Debrecen
   Hungary
*****************************************************************************