procrustes: Solving the Procrustes Problem

In pracma: Practical Numerical Math Functions

Solving the Procrustes Problem

Description

procrustes solves for two matrices A and B the ‘Procrustes Problem’ of finding an orthogonal matrix Q such that A-B*Q has the minimal Frobenius norm.

kabsch determines a best rotation of a given vector set into a second vector set by minimizing the weighted sum of squared deviations. The order of vectors is assumed fixed.

Usage

procrustes(A, B)

kabsch(A, B, w = NULL)

Arguments

A, B	two numeric matrices of the same size.
w	weights , influence the distance of points

Details

The function procrustes(A,B) uses the svd decomposition to find an orthogonal matrix Q such that A-B*Q has a minimal Frobenius norm, where this norm for a matrix C is defined as sqrt(Trace(t(C)*C)), or norm(C,'F') in R.

Solving it with B=I means finding a nearest orthogonal matrix.

kabsch solves a similar problem and uses the Procrustes procedure for its purpose. Given two sets of points, represented as columns of the matrices A and B, it determines an orthogonal matrix U and a translation vector R such that U*A+R-B is minimal.

Value

procrustes returns a list with components P, which is B*Q, then Q, the orthogonal matrix, and d, the Frobenius norm of A-B*Q.

kabsch returns a list with U the orthogonal matrix applied, R the translation vector, and d the least root mean square between U*A+R and B.

Note

The kabsch function does not take into account scaling of the sets, but this could easily be integrated.

References

Golub, G. H., and Ch. F. van Loan (1996). Matrix Computations. 3rd Edition, The John Hopkins University Press, Baltimore London. [Sect. 12.4, p. 601]

Kabsch, W. (1976). A solution for the best rotation to relate two sets of vectors. Acta Cryst A, Vol. 32, p. 9223.

See Also

svd

Examples

##  Procrustes
U <- randortho(5)               # random orthogonal matrix
P <- procrustes(U, eye(5))

##  Kabsch
P <- matrix(c(0, 1, 0, 0, 1, 1, 0, 1,
              0, 0, 1, 0, 1, 0, 1, 1,
              0, 0, 0, 1, 0, 1, 1, 1), nrow = 3, ncol = 8, byrow = TRUE)
R <- c(1, 1, 1)
phi <- pi/4
U <- matrix(c(1, 0, 0,
              0, cos(phi), -sin(phi),
              0, sin(phi),  cos(phi)), nrow = 3, ncol = 3, byrow = TRUE)

Q <- U %*% P + R
K <- kabsch(P, Q)
# K$R == R  and  K$U %*% P + c(K$R) == Q

pracma documentation built on Nov. 10, 2023, 1:14 a.m.