Procrustes Analysis

Steps in (home-made) procrustes analysis

Compare shape of data (represented as two matrices A and B) while ignoring the effects of translation, scaling and rotation.

Center A around 0: \[A_{centered} = A - Centroid(A) \]
Normalize the “size” of A by its Frobenius norm: \[A_{normalized} = A/(||A||_F/n) \]
Apply the same transformation to matrix B
Determine the best rotation matrix T, which rotates B to minimize the squared difference between A and B: \[ argmin_T (||A - TB||_F)\]

Mathematics for step 4

\[ \begin{aligned} argmin_T (||A - TB||_F) &= argmin_T (||A - TB||_F^2) \\ &= argmin_T (tr((A - TB)^T(A - TB)) \\ &= argmin_T (tr(A^TA + B^TT^TTB - 2A^TTB) \\ &= argmin_T (tr(A^TA) + tr(B^TB) - 2tr(A^TTB)) \\ &= argmax_T (tr(A^TTB)) \\ &= argmax_T (tr((BA^T)T)) \quad \quad \text{ (trace is invariant under cyclic permutations)}) \\ &= argmax_T (tr((U\Sigma V^T)T)) \quad \text{ (SVD of $BA^T$)}) \\ &= argmax_T (tr(\Sigma (V^TTU))) \quad \text{ (Elements in $V^TTU$ all less or equal to 1. Thus $tr(\Sigma (V^TTU)) \le tr(\Sigma I)$)} \\ &= VU^T \quad \quad \quad \quad \quad \quad\quad\quad \text{ (when $V^TTU = I$}) \\ \end{aligned} \]

Wrap into a function

procrustes <- function(A, B){
  # center and normalize A 
  A.centered <- t(scale(t(A), center = TRUE, scale = FALSE))
  A.size <- norm(A.centered, type = "F") / (ncol(A) * nrow(A))
  A.normalized <- A.centered / A.size

  # center and normalize B
  B.centered <- t(scale(t(B), center = TRUE, scale = FALSE))
  B.size <- norm(B.centered, type = "F") / (ncol(B) * nrow(B))
  B.normalized <- B.centered / B.size

  # Rotation matrix T 
  svd.results <- svd(B.normalized %*% t(A.normalized))
  U <- svd.results$u
  V <- svd.results$v
  T <- V %*% t(U)

  # B transformed
  B.transformed <- T %*% B.normalized

  # Error after superimposition
  RSS <- norm(A.normalized - B.transformed,  type = "F")

  # Return
  return(list(A.normalized = A.normalized, B.normalized = B.normalized, rotation.mtx = T, B.transformed = B.transformed, RSS = RSS))
}

Simulate data

# number of data points
k = 1000

# x
set.seed(7)
x <- rnorm(n = k, mean = 0, sd = 1)

# random noise
set.seed(9)
e <- rnorm(n = k, mean = 0, sd = 1)

# y
y <- abs(x) + e

# Matrix A
# Store x and y into a 2 x k matrix
A <- matrix(c(x, y), byrow = TRUE, nrow = 2)

# matrix B
# Arbitrarily rotate, scale, and translate A 
theta <- pi / 4
rot.mtx <- matrix(c(sin(theta), -cos(theta), cos(theta), sin(theta)), ncol=2)
B <- 0.5*rot.mtx %*% A - 6

Perform procrustes analysis on simulated data

procrustes.results <- procrustes(A, B)

Visualization

A.df <- as_data_frame(t(A))
B.df <- as_data_frame(t(B))

A.normalized.df <- as_data_frame(t(procrustes.results$A.normalized))
B.transformed.df <- as_data_frame(t(procrustes.results$B.transformed))


data.df <- rbind(A.df, B.df, A.normalized.df, B.transformed.df)
colnames(data.df) <- c('x', 'y')
data.df$matrix <- rep(c('A', 'B', 'A', 'B'), each = k)
data.df$treatment <- rep(c('Original', 'Superimposed'), each = 2*k)


data.plot <- ggplot(data.df, aes(x = x, y = y, color = matrix)) +
  geom_point(aes(shape = matrix), size = 4, alpha = 0.5) +
  theme_bw() +
  theme(aspect.ratio = 1) +
  facet_wrap(~treatment, scales = 'free') +
  labs(title = 'Procrustes analysis')

data.plot

Conclusion

Despite differences in original position, size, and orientation, A and B overlap very well after analysis since they have the same “shape”.

Reference

https://en.wikipedia.org/wiki/Procrustes_analysis https://en.wikipedia.org/wiki/Orthogonal_Procrustes_problem

Appendix

Error (Frobenius norm of the residual matrix)

procrustes.results$RSS

## [1] 2.062914e-12

Rotation matrix (Estimated vs. original)

t(procrustes.results$rotation.mtx)

##            [,1]      [,2]
## [1,]  0.7071068 0.7071068
## [2,] -0.7071068 0.7071068

rot.mtx

##            [,1]      [,2]
## [1,]  0.7071068 0.7071068
## [2,] -0.7071068 0.7071068