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=2000, 
       method="varimax", methodArgs=NULL, algorithm = "bb", fwindow = 10)
    GPFoblq(A, Tmat=diag(ncol(A)), normalize=FALSE, eps=1e-5, maxit=2000, 
       method="quartimin", methodArgs=NULL, algorithm = "bb", fwindow = 10)

    GPFRSorth(A, Tmat=NULL, normalize=FALSE, eps=1e-5, maxit=2000, method="varimax", 
       methodArgs=NULL, randomStarts=0, algorithm = "bb", fwindow = 10, ...)
    GPFRSoblq(A, Tmat=NULL, normalize=FALSE, eps=1e-5, maxit=2000, method="quartimin", 
       methodArgs=NULL, randomStarts=0, algorithm = "bb", fwindow = 10, ...)
    

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).
algorithm	Line-search strategy used to calculate step size alpha.
fwindow	integer. Number of previous criterion values used in the non-monotone line search.
...	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)).

The algorithm argument

The algorithm argument controls how the step size \alpha is initialized at each iteration of the gradient projection loop. Three options are available:

"legacy"

The original Bernaards and Jennrich (2005) heuristic: \alpha is doubled at the start of each iteration and halved until the Armijo sufficient decrease condition is satisfied. Simple and reliable, particularly for smooth orthogonal criteria. This is the default for all orthogonal rotations and for target and PST rotations.

"bb"

Barzilai-Borwein (BB) step size initialization. Uses the ratio of successive changes in the rotation matrix and projected gradient to estimate a good starting step size automatically:

\alpha_{BB} = \frac{\|\Delta T\|^2}{|\langle \Delta T, \Delta G_p \rangle|}

where \Delta T = T_k - T_{k-1} and \Delta G_p = G_{p,k} - G_{p,k-1} are the changes in the rotation matrix and projected gradient between consecutive iterations. The BB estimate is then used as the starting point for the Armijo line search. Benchmarking shows meaningful speedups for oblique rotation criteria on larger matrices (50 variables or more) with many random starts, particularly for non-smooth criteria such as simplimax. This is the default for all oblique rotations except target and PST.

"cayley"

The Cayley transform retraction for orthogonal rotations (Wen and Yin, 2013). Instead of projecting back onto the orthogonal manifold via SVD, the Cayley transform produces an exactly orthogonal update in closed form:

T_{k+1} = \left(I + \frac{\alpha}{2} W\right)^{-1} \left(I - \frac{\alpha}{2} W\right) T_k

where W = G_p T_k^T - T_k G_p^T is a skew-symmetric matrix derived from the projected gradient. The Cayley transform avoids the SVD decomposition at each iteration and can be faster for larger orthogonal problems. Available for orthogonal rotations only.

• algorithm = "bb" with fwindow = 10 is the default for all rotation criteria.

• algorithm = "cayley" is available for orthogonal rotations and can be specified explicitly by the user.

• algorithm = "legacy" is available for exact reproducibility of pre-2027 results.

Users can override the default by passing algorithm explicitly to any rotation function or directly to GPForth and GPFoblq.

The fwindow argument

The fwindow argument controls the width of the non-monotone line search window (Zhang and Hager, 2004). At each iteration, the Armijo sufficient decrease condition is evaluated against the maximum criterion value over the last fwindow iterations rather than just the immediately preceding value:

f(T_{k+1}) < \max_{j \in W_k} f(T_j) - 0.5\, s^2\, \alpha

where W_k is the window of the last fwindow iterates and s is the projected gradient norm.

• fwindow = 1 (default for "legacy") reduces to the standard monotone Armijo condition — the criterion must decrease at every iteration.

• fwindow > 1 allows temporary increases in the criterion value, which can help the algorithm escape flat regions and shallow local minima. This is particularly beneficial in combination with algorithm = "bb".

• fwindow = 10 is the default for algorithm = "bb" and algorithm = "cayley".

Benchmarking across matrix sizes and rotation criteria suggests that fwindow = 10 with algorithm = "bb" provides the best balance of speed and solution quality for oblique rotations on larger models. For orthogonal rotations with algorithm = "legacy", fwindow = 1 is recommended as larger windows can slow convergence.

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 when algorithm = "legacy" and fwindow = 1 are used.

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

Barzilai, J. and Borwein, J.M. (1988) Two-point step size gradient methods. IMA Journal of Numerical Analysis, 8, 141–148. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1093/imanum/8.1.141")}

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/0013164404272507")}

Jennrich, R.I. (2001) A simple general procedure for orthogonal rotation. Psychometrika, 66, 289–306. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BF02294840")}

Jennrich, R.I. (2002) A simple general method for oblique rotation. Psychometrika, 67, 7–19. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/BF02294706")}

Mansolf, M. and Reise, S.P. (2016) Exploratory Bifactor Analysis: The Schmid-Leiman Orthogonalization and Jennrich-Bentler Analytic Rotations. Multivariate Behavioral Research, 51(5), 698–717. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1037/met0000467")}

Wen, Z. and Yin, W. (2013) A feasible method for optimization with orthogonality constraints. Mathematical Programming, 142, 397–434. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s10107-012-0584-1")}

Zhang, H. and Hager, W.W. (2004) A nonmonotone line search technique and its application to unconstrained optimization. SIAM Journal on Optimization, 14(4), 1043–1056. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1137/S1052623403428208")}

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 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(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(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(fa.unrotated, normalize = TRUE, randomStarts = 100)
    data("CCAI", package = "GPArotation")
  res.u <- factanal(covmat = CCAI_R, factors = 4, rotation = "none",
                    n.obs = 461)

  # TandemI --- non-smooth orthogonal, BB and Cayley both faster
  res.t1l <- tandemI(res.u, algorithm = "legacy", fwindow =  1, maxit = 10000)
  res.t1b <- tandemI(res.u, algorithm = "bb",     fwindow = 10)
  res.t1c <- tandemI(res.u, algorithm = "cayley", fwindow = 10)
  cat("TandemI   legacy:", nrow(res.t1l$Table) - 1, "iterations\n")
  cat("TandemI   bb:    ", nrow(res.t1b$Table) - 1, "iterations",
      " max diff:", format(max(abs(res.t1l$loadings - res.t1b$loadings)),
                           scientific = TRUE, digits = 2), "\n")
  cat("TandemI   cayley:", nrow(res.t1c$Table) - 1, "iterations",
      " max diff:", format(max(abs(res.t1l$loadings - res.t1c$loadings)),
                           scientific = TRUE, digits = 2), "\n")

  # Entropy --- orthogonal, BB and Cayley both faster
  res.t2l <- entropy(res.u, algorithm = "legacy", fwindow =  1, maxit = 10000)
  res.t2b <- entropy(res.u, algorithm = "bb",     fwindow = 10)
  res.t2c <- entropy(res.u, algorithm = "cayley", fwindow = 10)
  cat("Entropy   legacy:", nrow(res.t2l$Table) - 1, "iterations\n")
  cat("Entropy   bb:    ", nrow(res.t2b$Table) - 1, "iterations",
      " max diff:", format(max(abs(res.t2l$loadings - res.t2b$loadings)),
                           scientific = TRUE, digits = 2), "\n")
  cat("Entropy   cayley:", nrow(res.t2c$Table) - 1, "iterations",
      " max diff:", format(max(abs(res.t2l$loadings - res.t2c$loadings)),
                           scientific = TRUE, digits = 2), "\n")

  # Oblimin --- smooth oblique, BB faster, same solution
  res.ol <- oblimin(res.u, algorithm = "legacy", fwindow =  1, maxit = 10000)
  res.ob <- oblimin(res.u, algorithm = "bb",     fwindow = 10)
  cat("Oblimin   legacy:", nrow(res.ol$Table) - 1, "iterations\n")
  cat("Oblimin   bb:    ", nrow(res.ob$Table) - 1, "iterations",
      " max diff:", format(max(abs(res.ol$loadings - res.ob$loadings)),
                           scientific = TRUE, digits = 2), "\n")

  # Simplimax --- non-smooth oblique with many local minima.
  # BB finds a substantially better solution in fewer iterations.
  # Difference in loadings reflects different local minima, not noise.
  res.sl <- simplimax(res.u, algorithm = "legacy", fwindow =  1, maxit = 10000)
  res.sb <- simplimax(res.u, algorithm = "bb",     fwindow = 10)
  cat("Simplimax legacy:", nrow(res.sl$Table) - 1, "iterations",
      " f =", round(res.sl$Table[nrow(res.sl$Table), 2], 6), "\n")
  cat("Simplimax bb:    ", nrow(res.sb$Table) - 1, "iterations",
      " f =", round(res.sb$Table[nrow(res.sb$Table), 2], 6), "\n")
  cat("BB found a lower criterion value --- better solution quality.\n")

GPArotation documentation built on June 18, 2026, 9:06 a.m.