how do I find the angles between an original and a rotated PCA loadings matrix?

Question

Suppose I have two matrices of PCA loadings loa.orig, and loa.rot, and I know that loa.rot is a rotation (by-hand or otherwise) of loa.orig.

(loa.orig might also have been already orthogonally rotated by varimax or something, but I don't think that matters).

I know want to know the angles by which loa.orig has been rotated to arrive at loa.rot.

I understand from this comment on another question that "rotations don't commute", and therefore, the order of pair-wise (plane-wise) rotations also matters.

So to reproduce loa.rot from loa.orig I would need to know a series of necessary rotations, ideally in the order given by rots in the below.

Here's an MWE:

library(psych) # this allows manual rotation

# suppose I have some ORIGINAL loadings matrix, from a principal components analysis, with three retained components
loa.orig <- cbind(c(0.6101496, 0.7114088, 0.3356003, 0.7318809, 0.5980133, 0.4102817, 0.7059148, 0.6080662, 0.5089014, 0.587025, 0.6166816, 0.6728603, 0.7482675, 0.5409658, 0.6415472, 0.3655053, 0.6313868), c(-0.205317, 0.3273207, 0.7551585, -0.1981179, -0.423377, -0.07281187, -0.04180098, 0.5003459, -0.504371, 0.1942334, -0.3285095, 0.5221494, 0.1850734, -0.2993066, -0.08715662, -0.02191772, -0.2002428), c(-0.4692407, 0.1581682, -0.04574932, -0.1189175, 0.2449018, -0.5283772, 0.02826476, 0.1703277, 0.2305158, 0.2135566, -0.2783354, -0.05187637, -0.104919, 0.5054129, -0.2403471, 0.5380329, -0.07999642))

# I then rotate 1-2 by 90°, and 1-3 by 45°
loa.rot <- factor.rotate(f = loa.orig, angle = 90, col1 = 1, col2 = 2)
loa.rot <- factor.rotate(f = loa.rot, angle = 45, col1 = 1, col2 = 3)

# predictably, loa.rot and loa.orig are now different
any(loa.rot == loa.orig) # are any of them the same?

Obviously, in this case, I know the angles and order, but let's assume I don't. Also, let's assume that in the real use case there may be many components retained and rotated, not just three.

I am a bit unsure what would be a conventional way to report the order of component-pair-wise (plane-wise) rotation angles, but I'd imagine the list of possible combinations (~~not permutations~~) should do.

combs <- combn(x = ncol(loa.orig), m = 2, simplify = TRUE)  # find all possible combinations of factors
rots <- data.frame(t(combs), stringsAsFactors = FALSE)  # transpose
rots  # these rows give the *order* in which the rotations should be done

rots gives these permutations.

Would be great to know how to arrive at loa.rot from loa.orig, rotating component pairs as given by the rows in rots.

---

UPDATE: Attempt based on below answer

Based on the below answer, I've tried to put together a function and test it with a varimax rotation and a real dataset. (No reason in particular for varimax – I just wanted some rotation where we actually do not know the angles.).

I then test whether I can actually re-create the varimax rotation from the vanilla loadings, using the extracted angles.

library(psych)
data("Harman74.cor")  # notice the correlation matrix is called "cov", but doc says its a cor matrix
vanilla <- principal(r = Harman74.cor$cov, nfactors = 4, rotate = "none", )$loadings  # this is unrotated
class(vanilla) <- NULL  # print methods only causes confusion
varimax <- principal(r = Harman74.cor$cov, nfactors = 4, rotate = "varimax")$loadings  # this is rotated
class(varimax) <- NULL  # print methods only causes confusion

find.rot.instr <- function(original, rotated) {
  # original <- vanilla$loadings  # testing
  # rotated <- varimax$loadings  # testing
  getAngle <- function(A, B) acos(sum(A*B) / (norm(A, "F") * norm(B, "F"))) * 180/pi

  rots <- combn(x = ncol(original), m = 2, simplify = FALSE)  # find all possible combinations of factor pairs

  tmp <- original
  angles <- sapply(rots, function(cols) {
    angle <- getAngle(tmp[, cols], rotated[, cols])
    tmp <<- factor.rotate(tmp, angle = angle, col1 = cols[1], col2 = cols[2])
    return(angle)
  })
  return(angles)
}

vanilla.to.varimax.instr <- find.rot.instr(original = vanilla, rotated = varimax)  # these are the angles we would need to transform in this order

rots <- combn(x = ncol(vanilla), m = 2, simplify = FALSE)  # find all possible combinations of factor pairs
# this is again, because above is in function

# now let's implement the extracted "recipe"
varimax.recreated <- vanilla # start with original loadings
varimax.recreated == vanilla  # confirm that it IS the same
for (i in 1:length(rots)) {  # loop over all combinations, starting from the top
  varimax.recreated <- factor.rotate(f = varimax.recreated, angle = vanilla.to.varimax.instr[i], col1 = rots[[i]][1], col2 = rots[[i]][2])
}
varimax == varimax.recreated  # test whether they are the same
varimax - varimax.recreated  # are the close?

Unfortunately, they are not the same, not even similar :(

> varimax == varimax.recreated  # test whether they are the same
PC1   PC3   PC2   PC4
VisualPerception       FALSE FALSE FALSE FALSE
Cubes                  FALSE FALSE FALSE FALSE
PaperFormBoard         FALSE FALSE FALSE FALSE
Flags                  FALSE FALSE FALSE FALSE
GeneralInformation     FALSE FALSE FALSE FALSE
PargraphComprehension  FALSE FALSE FALSE FALSE
SentenceCompletion     FALSE FALSE FALSE FALSE
WordClassification     FALSE FALSE FALSE FALSE
WordMeaning            FALSE FALSE FALSE FALSE
Addition               FALSE FALSE FALSE FALSE
Code                   FALSE FALSE FALSE FALSE
CountingDots           FALSE FALSE FALSE FALSE
StraightCurvedCapitals FALSE FALSE FALSE FALSE
WordRecognition        FALSE FALSE FALSE FALSE
NumberRecognition      FALSE FALSE FALSE FALSE
FigureRecognition      FALSE FALSE FALSE FALSE
ObjectNumber           FALSE FALSE FALSE FALSE
NumberFigure           FALSE FALSE FALSE FALSE
FigureWord             FALSE FALSE FALSE FALSE
Deduction              FALSE FALSE FALSE FALSE
NumericalPuzzles       FALSE FALSE FALSE FALSE
ProblemReasoning       FALSE FALSE FALSE FALSE
SeriesCompletion       FALSE FALSE FALSE FALSE
ArithmeticProblems     FALSE FALSE FALSE FALSE
> varimax - varimax.recreated  # are the close?
PC1         PC3          PC2       PC4
VisualPerception        0.2975463  1.06789735  0.467850675 0.7740766
Cubes                   0.2317711  0.91086618  0.361004861 0.4366521
PaperFormBoard          0.1840995  0.98694002  0.369663215 0.5496151
Flags                   0.4158185  0.82820078  0.439876777 0.5312143
GeneralInformation      0.8807097 -0.33385999  0.428455899 0.7537385
PargraphComprehension   0.7604679 -0.30162120  0.389727192 0.8329341
SentenceCompletion      0.9682664 -0.39302764  0.445263121 0.6673116
WordClassification      0.7714312  0.03747430  0.460461099 0.7643221
WordMeaning             0.8010876 -0.35125832  0.396077591 0.8201986
Addition                0.4236932 -0.32573100  0.204307400 0.6380764
Code                    0.1654224 -0.01757153  0.194533996 0.9777764
CountingDots            0.3585004  0.28032822  0.301148474 0.5929926
StraightCurvedCapitals  0.5313385  0.55251701  0.452293566 0.6859854
WordRecognition        -0.3157408 -0.13019630 -0.034647588 1.1235253
NumberRecognition      -0.4221889  0.10729098 -0.035324356 1.0963785
FigureRecognition      -0.3213392  0.76012989  0.158748259 1.1327322
ObjectNumber           -0.3234966 -0.02363732 -0.007830001 1.1804147
NumberFigure           -0.2033601  0.59238705  0.170467459 1.0831672
FigureWord             -0.0788080  0.35303097  0.154132395 0.9097971
Deduction               0.3423495  0.41210812  0.363022937 0.9181519
NumericalPuzzles        0.3573858  0.57718626  0.393958036 0.8206304
ProblemReasoning        0.3430690  0.39082641  0.358095577 0.9133117
SeriesCompletion        0.4933886  0.56821932  0.465602192 0.9062039
ArithmeticProblems      0.4835965 -0.03474482  0.332889805 0.9364874

So clearly, I'm making a mistake.

Catsup

Answer 1

I think the problem is easier if you do the matrix algebra (I am doing this for orthogonal rotations, not oblique transformations, but the logic is the same.)

First, notice that any rotation, t, produces a rotated components matrix c such that c = Ct where C is the unrotated PCA solution.

Then, since C'C is the original correlation matrix, R, we can solve for t by premultiplying both sides by C' and then by the inverse of R. This leads to

t = C' R^-1 c

For your example, let

R <- Harman74.cor$cov
P <- principal(Harman74.cor$cov,nfactors=4,rotate="none")
p <- principal(Harman74.cor$cov,nfactors=4)  #the default is to do varimax
rotation <- t(P$loadings) %*% solve(R) %*% p$loadings

Then, to see if this correct

P$loadings %*% rotation - p$loadings

In addition, the coming release of psych now reports the rotation matrix, so we can compare rotation to p$rot.mat

round(p$rot.mat,2)
round(rotation,2)

These differ in the order of the components, but that is because psych reorders the components after rotation to reflect the size of the rotated component sum of squares.

Maeve

thanks, this approach (via rotation matrices) really makes more sense to me, especially because rotation angles do commutate (order matter). I've also described how the newly returned rot.mat in psych solves this (and a related) problem: #32157531 – Catsup

Answer 2

EDIT: Method for arbitrary numbers of dimensions

I've got a method now which will find the analogue to Euler angles for an arbitrary number of dimensions for a rotation matrix (though it becomes increasingly computationally intensive as the dimensions increase). This method works for the sample dataset and varimax for between 2 and 6 factors. I haven't tested for more than 6. For 5 and 6 factors, there seems to be a reflection added for the 5th column - I'm not sure what determines which columns get reflected so this is hardcoded in at present into the examples.

The method works as before by generating a rotation matrix using lm.fit. It then uses yacas to calculate using symbolic matrix multiplication the composite rotation matrix for an arbitrary rotation in the appropriate number of dimensions. The angles are than solved iteratively. There will always be one element of the matrix which is based on the sin of one angle, and then this can be used to calculate the values of the other angles iteratively.

The output includes the subset of columns from the input that were actually different, the composite rotation/reflection matrix generated using a linear model, the individual rotations in terms of angles and columns, the calculated composite rotation matrix, a reflection/column swapping matrix needed for some of the examples, and the sum of the square of the differences between the calculated rotation and the input (for all of the examples, this is of the order of 1e-20).

Unlike my original solution, this only provides one possible combination of order of rotations. The actual number of possibilities (even just for the equivalent of Tait-Bryan angles is factorial(n * (n-1) / 2) for n dimensions, which for 6 dimensions is approx. 1.3e12.

library("psych")
library("magrittr")
library("stringr")
library("Ryacas")
rot_mat_nd <- function(dimensions, composite_var = NULL,
                       rot_order = combn(dimensions, 2, simplify = FALSE)) {
  d <- diag(dimensions)
  storage.mode(d) <- "character"
  mats <- lapply(seq(rot_order), function(i) {
    l <- paste0("a", i)
    cmb <- rot_order[[i]]
    d[cmb[1], cmb[1]] <- sprintf("cos(%s)", l)
    d[cmb[1], cmb[2]] <- sprintf("-sin(%s)", l)
    d[cmb[2], cmb[1]] <- sprintf("sin(%s)", l)
    d[cmb[2], cmb[2]] <- sprintf("cos(%s)", l)
    paste0("{{",
           paste(apply(d, 1, paste, collapse = ","), collapse = "},{"),
           "}}")
  })

  yac_statement <- paste0("Simplify(", paste(mats, collapse = "*"), ")")
  if (!is.null(composite_var)) {
    yac_statement <- paste0(composite_var, " := ", yac_statement)
  }
  output <- yacas(yac_statement)
  list(mats = mats, composite = output, rot_order = rot_order)
}


find_angles_nd <- function(input, rotated, reflect = NULL) {
  matched_cols <- sapply(1:ncol(input), function(i)
    isTRUE(all.equal(input[, i], rotated[, i])))
  dimensions <- sum(!matched_cols)
  theor_rot <- rot_mat_nd(dimensions, "r")
  rv <- yacas("rv := Concat @ r")
  swap_mat <- matrix(0, dimensions, dimensions)
  swap_mat[cbind(1:dimensions, match(colnames(input), colnames(rotated)))] <- 1
  if (!is.null(reflect)) {
    swap_mat[, reflect] <- -swap_mat[, reflect]
  }
  input_changed <- input[, !matched_cols]
  rotated_changed <- rotated[, !matched_cols]
  rotated_swapped <- rotated_changed %*% swap_mat
  rot_mat <- lm.fit(input_changed, rotated_swapped)$coef
  rot_mat_v <- c(t(rot_mat))

  known_angles <- numeric()
  angles_to_find <- nrow(rot_mat) * (nrow(rot_mat) - 1) / 2

  iterations <- 0L
  angles_found <- -1
  while(length(known_angles) < angles_to_find & angles_found != 0) {
    iterations <- iterations + 1L
    message(sprintf("Iteration %d; angles remaining at start %d",
                iterations, angles_to_find - length(known_angles)))
    yacas(sprintf("rvwv := WithValue({%s}, {%s}, rv)",
                  paste(names(known_angles), collapse = ", "),
                  paste(known_angles, collapse = ", ")
    ))
    var_num <- yacas("MapSingle(Length, MapSingle(VarList, rvwv))") %>%
      as.expression %>%
      eval %>%
      unlist
    angles_found <- 0L
    for (i in which(var_num == 1)) {
      var <- as.character(yacas(sprintf("VarList(rvwv[%d])[1]", i)))
      if (!(var %in% names(known_angles))) {
        to_solve <- as.character(yacas(sprintf("rvwv[%d]", i)))
        fun_var <- str_extract(to_solve, sprintf("(sin|cos)\\(%s\\)", var))
        fun_c <- substr(fun_var, 1, 3)
        if (fun_c == "sin") {
          to_solve_mod <- str_replace(to_solve, fixed(fun_var), "x")
          solved <- as.character(yacas(sprintf("Solve(%s == %0.15f, x)[1]", to_solve_mod, rot_mat_v[i])))
          answer <- asin(eval(parse(text = str_replace(solved, "x == ", ""))))
          known_angles <- c(known_angles, setNames(answer, var))
          angles_found <- angles_found + 1L
        }
      }
    }
    message(sprintf("- found %d", angles_found))
  }
  calc_rot_mat <-
    matrix(unlist(simplify2array(eval(
      as.expression(theor_rot$composite),
      as.list(known_angles)
    ))), dimensions, byrow = TRUE)

  ssd <- sum((input_changed %*% calc_rot_mat %*% swap_mat - rotated_changed) ^ 2)

  angles <- known_angles[paste0("a", 1:angles_to_find)] / pi * 180

  list(rot_mat = rot_mat, calc_rot_mat = calc_rot_mat, swap_mat = swap_mat, angles = angles,
       rot_order = theor_rot$rot_order, input_changed = input_changed,
       rotated_changed = rotated_changed, sum_square_diffs = ssd)
}

factor_rotate_multiple <- function(input_changed, angles, rot_order, swap_mat) {
  rotated <- input_changed
  for (i in 1:length(angles)) {
    rotated <- factor.rotate(rotated, angles[i], rot_order[[i]][1], rot_order[[i]][2])
  }
  rotated %*% swap_mat
}

Examples for 2-6 dimensions

data("Harman74.cor")  # notice the correlation matrix is called "cov", but doc says its a cor matrix

example_nd <- function(dimensions, reflect = NULL) {
  find_angles_nd(
    unclass(principal(r = Harman74.cor$cov, nfactors = dimensions, rotate = "none")$loadings),
    unclass(principal(r = Harman74.cor$cov, nfactors = dimensions, rotate = "varimax")$loadings),
    reflect = reflect
  )
}

angles_2d <- example_nd(2)

angles_2d[c("angles", "rot_order", "sum_square_diffs")]
# shows the resultant angle in degrees, rotation order and the
# sum of the squares of the differences between the provided and calculated
# rotations

#$angles
#       a1 
#-39.88672 
#
#$rot_order
#$rot_order[[1]]
#[1] 1 2
#
#
#$sum_square_diffs
#[1] 8.704914e-20

angles_3d <- example_nd(2)

angles_3d[c("angles", "rot_order", "sum_square_diffs")]
#$angles
#       a1        a2        a3 
#-45.19881 -29.77423 -17.07210 
#
#$rot_order
#$rot_order[[1]]
#[1] 1 2
#
#$rot_order[[2]]
#[1] 1 3
#
#$rot_order[[3]]
#[1] 2 3
#
#
#$sum_square_diffs
#[1] 7.498253e-20

angles_4d <- example_nd(2)
angles_5d <- example_nd(2, reflect = 5)
angles_6d <- example_nd(2, reflect = 5)

Original

This problem can be split into two. The first part is to calculate the composite rotation matrix that relates the input to the output. I.e. in the equation A %*% B == C, working out B from A and C. For square matrices, this can be done with solve. In this case, however, with rows > columns the easiest approach is to use linear models and the lm.fit function.

The second part of the problem is identifying the rotations that were performed to result in that composite rotation matrix. There are an infinite number of possible combinations, but in general it is reasonable to work (for 3 columns) in terms of a series of rotations about the three axes, i.e. using Euler angles. Even then there are six possible orders for the rotations. For two columns, the problem is trivial since there will only be one rotation required and a single asin or acos is enough to identify the angle. For more than 3 columns, it is possible to generalise the Euler angle method but gets more complicated.

Here's a complete method in R using linear models to find the composite rotation matrix and the RSpincalc package to find the Euler angles. It assumes that the rotations have affected 3 columns as in the example given.

library("RSpincalc")
library("combinat")
find_rotations <- function(input, rotated) {
  matched_cols <- sapply(1:ncol(input), function(i) isTRUE(all.equal(input[, i], rotated[, i])))
  if (sum(!matched_cols) != 3) {
    stop("This method only works for rotations affecting 3 columns.")
  }
  rot_mat <- lm.fit(input[, !matched_cols], rotated[, !matched_cols])$coef
  rot_poss <- as.data.frame(do.call("rbind", permn(c("z", "y", "x"))), stringsAsFactors = FALSE)
  rot_poss$axes_for_EA <- apply(rot_poss, 1, function(x) paste(rev(x), collapse = ""))
  combo_cols <- as.data.frame(matrix(which(!matched_cols)[combn(3, 2)], 2))
  rot_poss[5:10] <- do.call("rbind", permn(combo_cols, unlist))
  names(rot_poss)[c(1:3, 5:10)] <- c(paste0("axis", 1:3),
    paste0("rot", rep(1:3, each = 2), "_c", rep(1:2, 3)))
  rot_poss[paste0("angle", 1:3)] <- t(round(sapply(rot_poss$axes, DCM2EA, DCM = rot_mat), 14)) * 180 / pi
  rot_poss[paste0("angle", 1:3)] <-
    lapply(1:3, function(i) {
      ifelse(rot_poss[, paste0("axis", i)] == "y", 1, -1) *
        rot_poss[, paste0("angle", 4 - i)]
    })
  rot_poss
}

For the OP's data:

rot_poss <- find_rotations(loa.orig, loa.rot)

Another more comprehensive demo:

set.seed(123)
input <- matrix(rnorm(65), 13)

library("magrittr")

rotated <- input %>%
  factor.rotate(33.5, 1, 3) %>%
  factor.rotate(63, 2, 3) %>%
  factor.rotate(-3, 1, 2)

rot_poss <- find_rotations(input, rotated)

rot_poss

#  axis1 axis2 axis3 axes_for_EA rot1_c1 rot1_c2 rot2_c1 rot2_c2 rot3_c1 rot3_c2
#1     z     y     x         xyz       1       2       1       3       2       3
#2     z     x     y         yxz       1       2       2       3       1       3
#3     x     z     y         yzx       2       3       1       2       1       3
#4     x     y     z         zyx       2       3       1       3       1       2
#5     y     x     z         zxy       1       3       2       3       1       2
#6     y     z     x         xzy       1       3       1       2       2       3
#     angle1    angle2   angle3
#1 -1.585361 30.816825 63.84410
#2 44.624426 50.431683 53.53631
#3 59.538985 26.581047 16.27152
#4 66.980038 14.511490 27.52974
#5 33.500000 63.000000 -3.00000
#6 30.826477 -1.361477 63.03177

possible_calls <- do.call("sprintf", c(list(fmt = "input %%>%%
  factor.rotate(%1$g, %4$g, %5$g) %%>%%
  factor.rotate(%2$g, %6$g, %7$g) %%>%%
  factor.rotate(%3$g, %8$g, %9$g)"),
  rot_poss[c(paste0("angle", 1:3),
    grep("^rot", names(rot_poss), value = TRUE))]))
cat(possible_calls, sep = "\n\n")

#input %>%
  #factor.rotate(-1.58536, 1, 2) %>%
  #factor.rotate(30.8168, 1, 3) %>%
  #factor.rotate(63.8441, 2, 3)

#input %>%
  #factor.rotate(44.6244, 1, 2) %>%
  #factor.rotate(50.4317, 2, 3) %>%
  #factor.rotate(53.5363, 1, 3)

#input %>%
  #factor.rotate(59.539, 2, 3) %>%
  #factor.rotate(26.581, 1, 2) %>%
  #factor.rotate(16.2715, 1, 3)

#input %>%
  #factor.rotate(66.98, 2, 3) %>%
  #factor.rotate(14.5115, 1, 3) %>%
  #factor.rotate(27.5297, 1, 2)

#input %>%
  #factor.rotate(33.5, 1, 3) %>%
  #factor.rotate(63, 2, 3) %>%
  #factor.rotate(-3, 1, 2)

#input %>%
  #factor.rotate(30.8265, 1, 3) %>%
  #factor.rotate(-1.36148, 1, 2) %>%
  #factor.rotate(63.0318, 2, 3)

lapply(possible_calls, function(cl) all.equal(eval(parse(text = cl)), rotated, tolerance = 1e-6))
# tolerance reduced because above calls have rounding to 6 significant figures

Note the last bit will output the calls using magrittr pipes that can reproduce the rotation. Note also that there are 6 possible orders of rotation each with different angles.

In order to get the angles to match, I had to flip the sign on y rotation. I also had to flip the order of rotation.

For the varimax update, my method works for nfactor = 3 but would need adjustment for 4 or above.

Medication