Subject: Principal components summary

From: "[email-masked]" <[email-masked]>
Date: Tue, 29 Aug 1995 10:19:42 +1

Some days ago I sent the following summary to CCL, but somehow
only a small part has arrived. Here is it again.
---------------------------------------------------------------
Here is the summary of my question concerning the calculation of 
principal
components and plane fitting to points:

The 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 geo-
metry, 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.

Do not refer, please, to programs like Mathcad etc. I need it in
a numerical method form, as I'd like to use it as a procedure in a 
program.

***************************************************************
**
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)

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

From:             [email-masked] (Heather Gordon)
To:               [email-masked]
Subject:          pca

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
eigenvalues (the eigenvalues, eigenvectors 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]

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

Dear Tamas,
     I also am a chemist, a chemometrician in fact, and am interested 
in
chemistry
, statistics and computing all mixed togoether. Whenever I have a
mathematical problem like yours, I often find an answer in the 
excellent book
"Numerical Recipes in C: The Art of Scientific Computing" by W. 
H.Press,
S.A.Teukolsky, W.T.Vetterling and B.P.Flannery (Cambridge University
Press). There are also version for FORTRAN and, I think, BASIC. PCA 
is
also described in various multivariate statistics and chemometrics 
books, such
as "Multivariate Calibration" by H.Martens and T.Naes (Wiley) and
"Multivariate pattern recognition in chemometrics" edited by 
R.G.Brereton
(Elsevier).

     The source code is actually written in Visual Basic, but 
conversion to
Basic, Fortran or C should be quite straightforward

PCA: ijX = ikT * kjP 
where X is an i by j datamatrix, T is an i by k scores matrix and P 
is a k by
j loadings matrix. In your case, the number of components to extract 
("kmax"
in the source code) will be 2. For this listing, "evals" are the 
eigenvalues, one
for each PC, which I have defined as the sum of squares of the 
scores.
NIPALS calculates the PC's in order of importance (i.e. PC's with the 
biggest
eigenvalues come first). The loadings will be 2 vectors of unit 
length, which
represent the two new PC axis, and the scores matrix gives the 
coordinates of
the points on these new PC axis.

     Think about whether you wish to mean-centre the data prior to 
the
PCA or not - this will give different results!

     Maybe you have the answer to your problem anyway, by the time I
write this, but best of luck anyway!

******************************************************
Sub NIPALS (X!(), kmax%, mc%, T!(), P!(), evals!())

    Dim i%, imax%, j%, jmax%, k%, maxcol%
     imax = UBound(Xmat, 1)
     jmax = UBound(Xmat, 2)
    Dim sm#, ss#, cmaxss#, rmaxss#, smaxss#
    Dim cmax#(), rmax#(), smax#()

    ReDim cmax(1 To imax)
    ReDim rmax(1 To jmax)
    ReDim smax(1 To imax)

    ReDim T(1 To imax, 1 To kmax)
    ReDim P(1 To kmax, 1 To jmax)
    ReDim evals(1 To kmax)

    k = 1

    ' Mean-centre matrix if required
    If mc = 1 Then
     For j = 1 To jmax
         sm = 0

         For i = 1 To imax
          sm = sm + X(i, j)
         Next i

         sm = sm / imax
     Next j

     For j = 1 To jmax
         For i = 1 To imax
          X(i, j) = X(i, j) - sm
         Next i
     Next j
    End If

START1:

    ' Find column with greatest sum of squares
    cmaxss = 0#

    For j = 1 To jmax
     ss = 0#

     For i = 1 To imax
         ss = ss + X(i, j) ^ 2
     Next i
     
     If ss > cmaxss Then
         cmaxss = ss
         maxcol = j
     End If
    Next j

    For i = 1 To imax
     cmax(i) = X(i, maxcol)
    Next i

START2:

    ' Calculate row vector and its sum of squares
    rmaxss = 0#

    For j = 1 To jmax
     rmax(j) = 0#

     For i = 1 To imax
         rmax(j) = rmax(j) + (cmax(i) * X(i, j))
     Next i

     rmaxss = rmaxss + rmax(j) ^ 2
    Next j

    ' Scale this row vector to unit length, after checking that
it is not a
    ' zero vector
    If rmaxss = 0# Then
     MsgBox ("Nipals( ): Exited before all PC's found")
     Exit Sub
    End If

    For j = 1 To jmax
     rmax(j) = rmax(j) / Sqr(rmaxss)
    Next j

    ' Calculate a new estimate of the PC scores
    smaxss = 0#

    For i = 1 To imax
     smax(i) = 0#

     For j = 1 To jmax
         smax(i) = smax(i) + (X(i, j) * rmax(j))
     Next j

     smaxss = smaxss + smax(i) ^ 2
    Next i

    ' Test to see if PC has converged. If so, store scores and
loadings,
    ' adjust matrix, and look for next PC. If not, use last
scores estimate
    ' to put back into the algorithm. The convergence criterion
can be
    ' adjusted (i.e. 1e-6 used for single precision). 
    If Abs(smaxss - cmaxss) / cmaxss < .000000000001 Then
     evals(k) = smaxss

     For i = 1 To imax
         T(i, k) = smax(i)
     Next i

     For j = 1 To jmax
         P(k, j) = rmax(j)
     Next j

     For i = 1 To imax
         For j = 1 To jmax
          X(i, j) = X(i, j) - (smax(i) * rmax(j))
         Next j
     Next i

     If k = kmax Then
         Exit Sub
     End If

     k = k + 1

     GoTo START1
    Else
     For i = 1 To imax
         cmax(i) = smax(i)
     Next i

     cmaxss = smaxss

     GoTo START2
    End If

End Sub

================

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]

***************************************************************
**
 
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
***************************************************************
**
end summary

*****************************************************************************
   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
*****************************************************************************