GPA: Core Algorithms and Random-Start Wrappers

In GPArotation: Gradient Projection Factor Rotation

Core Algorithms and Random-Start Wrappers

Description

Gradient projection rotation optimization routines for orthogonal and oblique factor rotation. These functions can be used directly to rotate a loadings matrix, or indirectly through a rotation objective passed to a factor estimation routine such as factanal.

Usage

    GPForth(A, Tmat=diag(ncol(A)), normalize=FALSE, eps=1e-5, maxit=1000, 
       method="varimax", methodArgs=NULL)
    GPFoblq(A, Tmat=diag(ncol(A)), normalize=FALSE, eps=1e-5, maxit=1000, 
       method="quartimin", methodArgs=NULL)

    GPFRSorth(A, Tmat=diag(ncol(A)), normalize=FALSE, eps=1e-5, maxit=1000, 
       method="varimax", methodArgs=NULL, randomStarts=0, ...)
    GPFRSoblq(A, Tmat=diag(ncol(A)), normalize=FALSE, eps=1e-5, maxit=1000, 
       method="quartimin", methodArgs=NULL, randomStarts=0, ...)
    

Arguments

A	initial factor loadings matrix for which the rotation criterion is to be optimized.
Tmat	initial rotation matrix.
normalize	see details.
eps	convergence is assumed when the norm of the gradient is smaller than eps.
maxit	maximum number of iterations allowed in the main loop.
method	rotation objective criterion.
methodArgs	a list of additional arguments passed to the rotation objective.
randomStarts	number of random starts (GPFRSorth and GPFRSoblq only).
...	additional arguments passed to GPForth or GPFoblq, such as eps and maxit when calling via GPFRSorth or GPFRSoblq.

Details

The GPFRSorth and GPFRSoblq functions serve as the primary user interfaces for orthogonal and oblique rotations, respectively. They act as wrappers for the core GP algorithms (GPForth for orthogonal rotation and GPFoblq for oblique rotation), extending them with the ability to perform multiple random starts. Any additional arguments provided to these wrappers are passed directly down to the underlying GP algorithms. While the wrappers are generally recommended, the core functions GPForth and GPFoblq can also be invoked directly.

All of these functions require an initial loadings matrix, A, which fixes the equivalence class over which the optimization is performed. This matrix must be the solution to an orthogonal factor analysis problem, such as one obtained from factanal or another factor estimation routine.

Mathematically, a general rotation of a matrix A is defined as A %*% solve(t(Th)). In the case of orthogonal rotation, the initial rotation matrix Tmat is orthonormal, which simplifies the rotation formula to A %*% Th. In all scenarios, the final rotation matrix Th is computed by the GP rotation algorithm.

An accessible introduction to gradient projection algorithms for factor rotation is provided in Mansolf and Reise (2016).

The normalize argument

The normalize argument specifies whether and how the loadings matrix should be normalized prior to rotation, and subsequently denormalized after rotation.

• If FALSE (the default), no normalization is performed.

• If TRUE, Kaiser normalization is applied so that the squared row entries of the normalized matrix A sum to 1.0. This procedure is sometimes referred to as Horst normalization.

• If provided as a vector (which must have a length equal to the number of indicators, i.e., the number of rows in A), the columns of A are divided by this vector before rotation and multiplied by it afterward.

• If provided as a function, it can be used to apply a custom normalization scheme. The function must take A as an argument and return a vector, which is then applied in the same manner as the vector input described above. See NormalizingWeight for an example implementing Cureton-Mulaik normalization.

For a detailed investigation into how normalization affects factor rotations, including its potential impact on the qualitative interpretation of loadings, see Nguyen and Waller (2022).

The method argument

The method argument takes a string specifying the rotation objective function. By default, oblique rotations use "quartimin", while orthogonal rotations default to "varimax". The package supports a comprehensive suite of rotation objectives: "oblimin", "quartimin", "target", "pst", "oblimax", "entropy", "quartimax", "Varimax", "simplimax", "bentler", "tandemI", "tandemII", "geomin", "cf", "infomax", "mccammon", "bifactor", "lp", and "varimin".

Internally, this string is prefixed with "vgQ." to invoke the actual calculation function (see vgQ for underlying mathematical details). It is important to note that several rotation criteria—specifically "oblimin", "target", "pst", "simplimax", "geomin", "cf", and "lp"—require one or more supplementary arguments. These additional arguments can be seamlessly passed via the methodArgs list in the wrapper functions. Default values and direct usage examples for these arguments can be found in the rotations documentation.

The randomStarts argument

Because factor rotation criteria frequently suffer from local minima, trying multiple starting configurations can help identify a superior solution. The randomStarts argument, available exclusively in the GPFRSorth and GPFRSoblq wrappers, facilitates this robust search approach.

• By default, randomStarts = 0, which defaults to using the identity matrix as the initial rotation matrix Tmat. The initial rotation matrix Tmat can also be set by the user.

• Setting randomStarts = 1 initializes Tmat with a single random matrix.

• Setting randomStarts > 1 attempts multiple random starts and returns the rotated loadings matrix that achieved the lowest criterion value f across all attempts. Note that this returned solution is technically still a local minimum, and is not guaranteed to be the global minimum. Users are encouraged to review the random start diagnostics detailed in the package examples.

Under the hood, an internal, unexported engine named .GPA_RS_engine safely manages the random start loop, tracks convergence diagnostics, and handles factor correlation matrix naming.

While the core algorithms GPForth and GPFoblq do not support the randomStarts argument directly, users can manually supply a single random initial rotation matrix to them using Tmat = Random.Start(ncol(A)).

Legacy functions

The original implementations authored by Bernaards and Jennrich (2005) have been retained as GPForth.legacy and GPFoblq.legacy. These functions are kept purely for historical reference and backward compatibility for reproducibility. They are not exported into the package namespace, meaning they must be explicitly invoked using the triple-colon operator:

GPArotation:::GPForth.legacy(A, method = "varimax")
GPArotation:::GPFoblq.legacy(A, method = "quartimin")

The results generated by these legacy functions should be numerically identical to those produced by the current implementations. You can see a direct comparison of this in the examples section.

Value

A GPArotation object which is a list with elements:

loadings	The rotated loadings matrix, one column per factor. If random starts were requested, this is the solution with the lowest criterion value.
Th	The rotation matrix, satisfying loadings %*% t(Th) = A for orthogonal rotation and loadings = A %*% solve(t(Th)) for oblique rotation.
Table	A matrix recording the iteration history: iteration number, criterion value, log10 of the gradient norm, and step size (alpha).
method	A string indicating the rotation criterion.
orthogonal	A logical indicating if the rotation is orthogonal.
convergence	A logical indicating if convergence was obtained.
Phi	t(Th) %*% Th, the covariance matrix of the rotated factors. Omitted (NULL) for orthogonal rotations.
Gq	The gradient of the criterion at the rotated loadings.
randStartChar	A named vector summarising random start results: randomStarts, Converged, atMinimum, localMins. Only present when randomStarts > 1.

Author(s)

Coen A. Bernaards and Robert I. Jennrich with some R modifications by Paul Gilbert

References

Bernaards, C.A. and Jennrich, R.I. (2005) Gradient Projection Algorithms and Software for Arbitrary Rotation Criteria in Factor Analysis. Educational and Psychological Measurement, 65, 676–696. doi: 10.1177/0013164404272507

Jennrich, R.I. (2001). A simple general procedure for orthogonal rotation. Psychometrika, 66, 289–306. doi: 10.1007/BF02294840

Jennrich, R.I. (2002). A simple general method for oblique rotation. Psychometrika, 67, 7–19. doi: 10.1007/BF02294706

Mansolf, M., & Reise, S. P. (2016). Exploratory Bifactor Analysis: The Schmid-Leiman Orthogonalization and Jennrich-Bentler Analytic Rotations. Multivariate Behavioral Research, 51(5), 698–717. doi: 10.1080/00273171.2016.1215898

Nguyen, H.V. and Waller, N.G. (2023). Local minima and factor rotations in exploratory factor analysis. Psychological Methods. 28(5), 1122–1141. doi: 10.1037/met0000467

See Also

rotations, Random.Start, factanal, Harman8, box26, CCAI, NetherlandsTV

Examples

  # --- Basic rotation calls ---
  data(Harman, package = "GPArotation")           # 8 physical variables
  quartimax(Harman8)                              # direct rotation call
  GPFRSorth(Harman8, method = "quartimax")        # equivalent via wrapper
  GPFRSoblq(Harman8, method = "quartimin", normalize = TRUE)
  loadings(quartimin(Harman8, normalize = TRUE))  # extract loadings directly

  # --- Passing criterion arguments via methodArgs ---
  # Crawford-Ferguson family: kappa selects the criterion.
  # For box26: p = 26 variables, m = 3 factors.
  # Equamax:  kappa = m / (2 * p) = 3 / 52
  # Parsimax: kappa = (m - 1) / (p + m - 2) = 2 / 27
  data(Thurstone, package = "GPArotation") # 26 variable box problem 
  GPFRSoblq(box26, method = "cf", methodArgs = list(kappa = 3/52))  # Equamax
  GPFRSoblq(box26, method = "cf", methodArgs = list(kappa = 2/27))  # Parsimax
  
  # --- Two-step vs single-step factanal for oblique rotation ---
  #
  # The recommended approach for oblique rotation is the two-step procedure:
  # (1) obtain unrotated loadings from factanal, then
  # (2) rotate separately using GPArotation.
  # This gives full control over the rotation, including random starts.
  #
  # Prior to R 4.5.1, the single-step approach (rotation inside factanal)
  # had a bug in factor reordering after oblique rotation. This was reported
  # by Bernaards and others and fixed by the R core team in R 4.5.1.

  data("WansbeekMeijer", package = "GPArotation")

  # Step 1: unrotated 3-factor solution
  fa.unrotated <- factanal(factors = 3, covmat = NetherlandsTV,
                           normalize = TRUE, rotation = "none")

  # Step 2: oblique Crawford-Ferguson rotation with kappa = 0.3
  # (non-standard kappa, not corresponding to any named special case)
  set.seed(44)
  fa.cf <- cfQ(loadings(fa.unrotated), kappa = 0.3, normalize = TRUE,
               randomStarts = 100)
  fa.cf

  # Single-step via factanal - correct in R >= 4.5.1
  if (getRversion() >= "4.5.1") {
    set.seed(44)
    fa.factanal <- factanal(factors = 3, covmat = NetherlandsTV, rotation = "cfQ",
      control = list(rotate = list(normalize = TRUE, kappa = 0.3, randomStarts = 100)))
    # The two approaches should agree after sorting
    fa.sorted <- print(fa.cf, sortLoadings = TRUE)
    cat("Maximum difference in loadings between two-step and single-step:\n")
    print(max(abs(abs(fa.sorted$loadings) - abs(fa.factanal$loadings))))
  } else {
    cat("Single-step factanal oblique rotation requires R >= 4.5.1.\n")
    cat("Use the two-step procedure above for correct results.\n")
  }

  # --- Displaying rotation output ---
  origdigits <- options("digits")
  data("CCAI", package = "GPArotation")
  fa.unrotated <- factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "none")
  res <- oblimin(loadings(fa.unrotated), gam = -0.5, randomStarts = 20)
                                    # gam = -0.5: more orthogonal than quartimin
  res                               # default print
  print(res)                        # equivalent to above
  print(res, Table = TRUE)          # include iteration table
  print(res, rotateMat = TRUE)      # include rotating matrix
  print(res, digits = 2)            # rounded to 2 decimal places
  summary(res)                      # pattern and structure matrices for oblique rotation
  summary(res, Structure = FALSE)   # pattern matrix only
  options(digits = origdigits$digits)
    
  # --- Random start diagnostics ---
  # When randomStarts > 1, the output includes randStartChar which summarizes
  # the random start results:
  #   randomStarts : number of random starts attempted
  #   Converged    : number of starts that converged
  #   atMinimum    : number of starts at the same lowest minimum
  #   localMins    : number of distinct local minima found
  data(Thurstone, package = "GPArotation")
  res <- GPFRSoblq(box26, method = "geomin", normalize = TRUE, randomStarts = 50)
  res$randStartChar

  # --- Factor ordering ---
  # Raw GPArotation output is unsorted — factors may appear in any order
  # depending on the starting matrix. Use print() to obtain sorted loadings.
  # Once sorted, repeated calls to print() are stable.
  set.seed(334)
  xusl <- quartimin(Harman8, normalize = TRUE, randomStarts = 100)

  loadings(xusl)                                        # unsorted raw output
  max(abs(print(xusl)$loadings - xusl$loadings)) == 0  # FALSE: print() reorders
  xsl <- print(xusl)                                    # capture sorted result
  max(abs(print(xsl)$loadings - xsl$loadings)) == 0    # TRUE: already sorted
  
  # --- Normalization ---
  # Kaiser normalization
  data("CCAI", package = "GPArotation")
  fa.unrotated <- factanal(factors = 3, covmat = CCAI_R, n.obs = 461, rotation = "none")
  oblimin(loadings(fa.unrotated), normalize = TRUE, randomStarts = 100)
  
  # --- Legacy engine --- 
  # Demonstrates numerical equivalence of current and legacy implementations.
  # The maximum absolute difference should be zero or at machine epsilon.
  data("Harman", package = "GPArotation")

  # Rotate using current engine
  res_new    <- oblimin(Harman8)

  # Rotate using legacy engine
  res_legacy <- GPArotation:::GPFoblq.legacy(Harman8, method = "oblimin")

  # Numerically identical
  max(abs(loadings(res_new) - loadings(res_legacy)))

GPArotation documentation built on April 29, 2026, 9:08 a.m.