NAG Library Routine Document
g03bcf (rot_procrustes)
1
Purpose
g03bcf computes Procrustes rotations in which an orthogonal rotation is found so that a transformed matrix best matches a target matrix.
2
Specification
Fortran Interface
Subroutine g03bcf ( 
stand, pscale, n, m, x, ldx, y, ldy, yhat, r, ldr, alpha, rss, res, wk, ifail) 
Integer, Intent (In)  ::  n, m, ldx, ldy, ldr  Integer, Intent (Inout)  ::  ifail  Real (Kind=nag_wp), Intent (Inout)  ::  x(ldx,m), y(ldy,m), yhat(ldy,m), r(ldr,m)  Real (Kind=nag_wp), Intent (Out)  ::  alpha, rss, res(n), wk(m*m+7*m)  Character (1), Intent (In)  ::  stand, pscale 

C Header Interface
#include nagmk26.h
void 
g03bcf_ (const char *stand, const char *pscale, const Integer *n, const Integer *m, double x[], const Integer *ldx, double y[], const Integer *ldy, double yhat[], double r[], const Integer *ldr, double *alpha, double *rss, double res[], double wk[], Integer *ifail, const Charlen length_stand, const Charlen length_pscale) 

3
Description
Let $X$ and $Y$ be $n$ by $m$ matrices. They can be considered as representing sets of $n$ points in an $m$dimensional space. The $X$ matrix may be a matrix of loadings from say factor or canonical variate analysis, and the $Y$ matrix may be a postulated pattern matrix or the loadings from a different sample. The problem is to relate the two sets of points without disturbing the relationships between the points in each set. This can be achieved by translating, rotating and scaling the sets of points. The $Y$ matrix is considered as the target matrix and the $X$ matrix is rotated to match that matrix.
First the two sets of points are translated so that their centroids are at the origin to give
${X}_{c}$ and
${Y}_{c}$, i.e., the matrices will have zero column means. Then the rotation of the translated
${X}_{c}$ matrix which minimizes the sum of squared distances between corresponding points in the two sets is found. This is computed from the singular value decomposition of the matrix:
where
$U$ and
$V$ are orthogonal matrices and
$D$ is a diagonal matrix. The matrix of rotations,
$R$, is computed as:
After rotation, a scaling or dilation factor,
$\alpha $, may be estimated by least squares. Thus, the final set of points that best match
${Y}_{c}$ is given by:
Before rotation, both sets of points may be normalized to have unit sums of squares or the
$X$ matrix may be normalized to have the same sum of squares as the
$Y$ matrix. After rotation, the results may be translated to the original
$Y$ centroid.
The $i$th residual, ${r}_{i}$, is given by the distance between the point given in the $i$th row of $Y$ and the point given in the $i$th row of $\hat{Y}$. The residual sum of squares is also computed.
4
References
Krzanowski W J (1990) Principles of Multivariate Analysis Oxford University Press
Lawley D N and Maxwell A E (1971) Factor Analysis as a Statistical Method (2nd Edition) Butterworths
5
Arguments
 1: $\mathbf{stand}$ – Character(1)Input

On entry: indicates if translation/normalization is required.
 ${\mathbf{stand}}=\text{'N'}$
 There is no translation or normalization.
 ${\mathbf{stand}}=\text{'Z'}$
 There is translation to the origin (i.e., to zero).
 ${\mathbf{stand}}=\text{'C'}$
 There is translation to origin and then to the $Y$ centroid after rotation.
 ${\mathbf{stand}}=\text{'U'}$
 There is unit normalization.
 ${\mathbf{stand}}=\text{'S'}$
 There is translation and normalization (i.e., there is standardization).
 ${\mathbf{stand}}=\text{'M'}$
 There is translation and normalization to $Y$ scale, then translation to the $Y$ centroid after rotation (i.e., they are matched).
Constraint:
${\mathbf{stand}}=\text{'N'}$, $\text{'Z'}$, $\text{'C'}$, $\text{'U'}$, $\text{'S'}$ or $\text{'M'}$.
 2: $\mathbf{pscale}$ – Character(1)Input

On entry: indicates if least squares scaling is to be applied after rotation.
 ${\mathbf{pscale}}=\text{'S'}$
 Scaling is applied.
 ${\mathbf{pscale}}=\text{'U'}$
 No scaling is applied.
Constraint:
${\mathbf{pscale}}=\text{'S'}$ or $\text{'U'}$.
 3: $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of points.
Constraint:
${\mathbf{n}}\ge {\mathbf{m}}$.
 4: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of dimensions.
Constraint:
${\mathbf{m}}\ge 1$.
 5: $\mathbf{x}\left({\mathbf{ldx}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: $X$, the matrix to be rotated.
On exit: if
${\mathbf{stand}}=\text{'N'}$,
x will be unchanged.
If
${\mathbf{stand}}=\text{'Z'}$,
$\text{'C'}$,
$\text{'S'}$ or
$\text{'M'}$,
x will be translated to have zero column means.
If
${\mathbf{stand}}=\text{'U'}$ or
$\text{'S'}$,
x will be scaled to have unit sum of squares.
If
${\mathbf{stand}}=\text{'M'}$,
x will be scaled to have the same sum of squares as
y.
 6: $\mathbf{ldx}$ – IntegerInput

On entry: the first dimension of the array
x as declared in the (sub)program from which
g03bcf is called.
Constraint:
${\mathbf{ldx}}\ge {\mathbf{n}}$.
 7: $\mathbf{y}\left({\mathbf{ldy}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: the target matrix, ${\mathbf{y}}$.
On exit: if
${\mathbf{stand}}=\text{'N'}$,
y will be unchanged.
If
${\mathbf{stand}}=\text{'Z'}$ or
$\text{'S'}$,
y will be translated to have zero column means.
If
${\mathbf{stand}}=\text{'U'}$ or
$\text{'S'}$,
y will be scaled to have unit sum of squares.
If
${\mathbf{stand}}=\text{'C'}$ or
$\text{'M'}$,
y will be translated and then after rotation translated back. The output
y should be the same as the input
y except for rounding errors.
 8: $\mathbf{ldy}$ – IntegerInput

On entry: the first dimension of the arrays
y and
yhat as declared in the (sub)program from which
g03bcf is called.
Constraint:
${\mathbf{ldy}}\ge {\mathbf{n}}$.
 9: $\mathbf{yhat}\left({\mathbf{ldy}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the fitted matrix, $\hat{Y}$.
 10: $\mathbf{r}\left({\mathbf{ldr}},{\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the matrix of rotations,
$R$, see
Section 9.
 11: $\mathbf{ldr}$ – IntegerInput

On entry: the first dimension of the array
r as declared in the (sub)program from which
g03bcf is called.
Constraint:
${\mathbf{ldr}}\ge {\mathbf{m}}$.
 12: $\mathbf{alpha}$ – Real (Kind=nag_wp)Output

On exit: if
${\mathbf{pscale}}=\text{'S'}$ the scaling factor,
$\alpha $; otherwise
alpha is not set.

On exit: the residual sum of squares.
 14: $\mathbf{res}\left({\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the residuals,
${r}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n$.
 15: $\mathbf{wk}\left({\mathbf{m}}\times {\mathbf{m}}+7\times {\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayWorkspace

 16: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{n}}<{\mathbf{m}}$, 
or  ${\mathbf{m}}<1$, 
or  ${\mathbf{ldx}}<{\mathbf{n}}$, 
or  ${\mathbf{ldy}}<{\mathbf{n}}$, 
or  ${\mathbf{ldr}}<{\mathbf{m}}$, 
or  ${\mathbf{stand}}\ne \text{'N'}$, $\text{'Z'}$, $\text{'C'}$, $\text{'U'}$, $\text{'S'}$ or $\text{'M'}$, 
or  ${\mathbf{pscale}}\ne \text{'S'}$ or $\text{'U'}$. 
 ${\mathbf{ifail}}=2$

On entry, either
x or
y contain only zeropoints (possibly after translation) when normalization is to be applied.
 ${\mathbf{ifail}}=3$

The $\hat{Y}$ matrix contains only zeropoints when least squares scaling is applied.
 ${\mathbf{ifail}}=4$

The singular value decomposition has failed to converge. This is an unlikely error exit.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
The accuracy of the calculation of the rotation matrix largely depends upon the singular value decomposition. See the
F08 Chapter Introduction for further details.
8
Parallelism and Performance
g03bcf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the
X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
Note that if the matrix
${X}_{c}^{\mathrm{T}}Y$ is not of full rank, then the matrix of rotations,
$R$, may not be unique even if there is a unique solution in terms of the rotated matrix,
${\hat{Y}}_{c}$. The matrix
$R$ may also include reflections as well as pure rotations, see
Krzanowski (1990).
If the column dimensions of the $X$ and $Y$ matrices are not equal, the smaller of the two should be supplemented by columns of zeros. Adding a column of zeros to both $X$ and $Y$ will have the effect of allowing reflections as well as rotations.
10
Example
Three points representing the vertices of a triangle in two dimensions are input. The points are translated and rotated to match the triangle given by $\left(0,0\right)$, $\left(1,0\right)$, $\left(0,2\right)$ and scaling is applied after rotation. The target matrix and fitted matrix are printed along with additional information.
10.1
Program Text
Program Text (g03bcfe.f90)
10.2
Program Data
Program Data (g03bcfe.d)
10.3
Program Results
Program Results (g03bcfe.r)