Rotation to get equal loadings in the first principal component or factor

Question

I have the observations of $n$ variables $x_i(t)$ where $t$ is the time, and $i=1,2,\dots,n$ is the number of the variable. They're very correlated, so I wanted to use PCA. The first component explains a ton of variance, and its loadings are sort of equal $w_i\approx1/n$, but not exactly equal. I wish they were exactly equal!

So, I'd like to rotate the factors such that the first component loadings would be exactly equal, i.e. $w_i=1/n$. Obviously, this is going to mess up other components (factors), but it's probably Ok. My question is how to rotate PCA factors or set up factor analysis so that my first factor would have equal loadings, and the other components were similar to PC2 and PC3 from PCA?

I'm using MATLAB, but any other language would be Ok with me.

UPDATE: I want my factors still be uncorrelated as in PCA. I'm assuming that since PC1 is already quite close to equal weights, then it must be possible to distort PC2 and PC3 as less as possible. PC2 and PC3 capture the slope and the curvature of variables, so it would be good to preserve this quality, but it's not essential.

Answer 1

The standard "factor rotation" take on this problem would be to perform PCA, choose some number of components (e.g. three), normalize them to unit variance, and then rotate the loadings according to some specific criterium. Usually rotation criteria optimize loading sparseness, but I guess one can try instead to bring the first loading vector to a specific direction -- in your case $$\mathbf v_1 = [1/n, \ldots, 1/n]^\top.$$ However, note that the subspace spanned by the first three PC directions might not include $\mathbf v_1$ at all! Most likely it does not. If so, then there is no "factor rotation" sensu strictu that can solve your problem.

For that reason I think you need to use a more specific approach.

---

My suggestion is to construct your first component (C1, as opposed to PC1) by projecting the data $\mathbf X$ onto the main diagonal vector $\mathbf v_1$ by computing $\mathbf X \mathbf v_1$. Then you can take the "lefover" data and perform a standard PCA on it, taking resulting principal components as your components C2, C3, etc.

There are two natural ways to define "leftover" data, depending on what property of PCA you want to be preserved.

Preserving the orthogonality of the PC directions

Project the data on the hyperplane orthogonal to $\mathbf v_1$, i.e. by defining $$\mathbf X_\text{leftover}=\mathbf X-\mathbf X\mathbf v_1 \mathbf v_1^\top = \mathbf X(\mathbf I - \mathbf v_1 \mathbf v_1^\top).$$ This will produce orthogonal component axes (because PCA will be done in the subspace orthogonal to $\mathbf v_1$), but not necessarily uncorrelated components: C2 and C3 can (and will) be somehow correlated with C1. If the original PCA produces PC1 that is very close to C1, then these correlations will be weak.

Preserving the uncorrelatedness of the PCs

The set of vectors $\mathbf w$ that are orthogonal to the unit vector $\mathbf v_1$ is defined by the condition $\mathbf w^\top \mathbf v_1$; to project the data onto the subspace spanned by this set we used the projector $\mathbf P = \mathbf I - \mathbf v_1 \mathbf v^\top$.

The set of vectors $\mathbf w$ that yield projections uncorrelated with the projection onto $\mathbf v_1$ is defined by the condition $\mathbf w^\top \boldsymbol \Sigma \mathbf v_1$, where I defined the covariance matrix of the data $\boldsymbol \Sigma = \mathbf X^\top \mathbf X / (n-1)$. It follows by analogy that to project the data onto the subspace spanned by this set we need to use the projector $\mathbf P = \mathbf I - \tilde{\mathbf v}_1 \tilde{\mathbf v}_1^\top$, where $\tilde{\mathbf v}_1=\boldsymbol \Sigma \mathbf v_1/\|\boldsymbol \Sigma \mathbf v_1\|$ is simply $\boldsymbol \Sigma \mathbf v_1$ normalized to unit length.

This insight allows us to compute the leftover data as follows: $$\mathbf X_\text{leftover}=\mathbf X\left(\mathbf I - \frac{\boldsymbol \Sigma \mathbf v_1 \mathbf v_1^\top \boldsymbol \Sigma}{\|\boldsymbol \Sigma \mathbf v_1\|^2}\right).$$

Again, this will yield uncorrelated components, but not necessarily orthogonal projection directions.

---

Note, as a consistency check, that if $\mathbf v_1$ is a unit eigenvector of $\boldsymbol \Sigma$, then $\tilde{\mathbf v}_1=\boldsymbol \Sigma \mathbf v_1/\|\boldsymbol \Sigma \mathbf v_1\| = \mathbf v_1$, meaning that for standard PCA the two approaches coincide.

Answer 2

In my matrix-calculator "MatMate" I've implemented some "experimental" rotations where I play with the usual orthogonal rotation-criteria "pca", "quartimax", "varimax", giving them a minus-sign to minimize instead of maximizing: "-pca", "-quartimax", "-varimax" .
The "-varimax" comes near to what you want: it finds that coordinate-system where the variance of the squared loadings is minimized. Unfortunately in the example dataset this does not arrive at a factor with exactly equal loadings.

But I've also one, which maximizes/minimizes the variance of the unsquared loadings - then the first factor has maximal variance in its loadings and the last factor minimal variance - which means, the loadings are in fact equal. I've called this provisorically "medimax" - rotation (the rotation-criterion is not more difficult that that of "varimax" and should be easily implementable in something like Matlab, I think) .

If you like you can experiment with MatMate after downloading, see here; your data/correlation-/loadingsmatrices can (possibly easiest) be imported via clipboard as a tab-or blank-separated number-array.

Here are some sample rotation of a 5x5 loadingsmatrix.

load = rot(load,"drei")   // "drei"="dreieck" means here "triangular"
  // this is also my initial loadingsmatrix taken by cholesky-cecomposition
  // of the correlation matrix, say "cor":
  // load = cholesky(cor)
            f1         f2           f3          f4         f5 
item1   |   1,0000    0,0000      0,0000      0,0000      0,0000 |
item2   |   0,5413    0,8408      0,0000      0,0000      0,0000 |
item3   |  -0,1765   -0,3278      0,9281      0,0000      0,0000 |
item4   |   0,4574    0,2695     -0,1543      0,8333      0,0000 |
item5   |   0,4155    0,1662     -0,1271      0,3301      0,8214 |

load = rot(load,"pca") 
        // in MatMate the usual pc-solution is taken by rotation! ( with "pc-criterion")
        |   0,7419    0,3537     -0,3776      0,2203      0,3651 |
        |   0,7796   -0,0712     -0,4365     -0,0794     -0,4362 |
        |  -0,5322    0,8188     -0,0462     -0,1317     -0,1637 |
        |   0,7918    0,0927      0,2645     -0,5233      0,1434 |
        |   0,7200    0,2159      0,5367      0,3371     -0,1827 |

load = rot(load,"varimax")
  // the rotation with the varimax-criterion shows, that the items
  // are nearly uncorrelated: each one has its own factor/axis to load high on  
        |   0,0535   -0,9316     -0,2481      0,1878     -0,1803 |
        |   0,1839   -0,2595     -0,9166      0,1967     -0,1409 |
        |  -0,9748    0,0498      0,1550     -0,1203      0,0939 |
        |   0,1414   -0,1965     -0,1979      0,9152     -0,2543 |
        |   0,1024   -0,1753     -0,1325      0,2350     -0,9413 |

load = rot(load,"-varimax")
  // this is one of the "experimental" rotations: minimizing the variance
  // of the squared loadings. We do not find an axis/factor with perfect
  // equal loadings though.
        |   0,6802   -0,3249     -0,3772     -0,3569      0,4025 |
        |   0,7219    0,4043     -0,3121      0,2589      0,3883 |
        |  -0,7161   -0,4121     -0,3744      0,2675      0,3251 |
        |   0,7163   -0,3474      0,3742      0,3578      0,3133 |
        |   0,6839   -0,4007     -0,2887      0,3670     -0,3921 |

Now we introduce the obscure "medimax"-rotation:

load = rot(load,"medimax")
  // here we have another experimental rotation, which finds indeed
  // a factor of equal loadings. Here it is the last factor (while
  // the first factor has maximal possible variance on it
        |   0,5054   -0,3506     -0,4818     -0,3938     -0,4843 |
        |   0,6789   -0,4363      0,1422      0,3065     -0,4843 |
        |  -0,8733   -0,0104     -0,0146     -0,0498     -0,4843 |
        |   0,6403    0,2734      0,3413     -0,4052     -0,4843 |
        |   0,5388    0,5529     -0,3663      0,1879     -0,4843 |

load = rot(load,"-medimax")
  // the minus-sign /inversion of the rotation-criterion leads here only to 
  // inversion of the order of the factors.
        |   0,4843   -0,3938     -0,4818      0,3506      0,5054 |
        |   0,4843    0,3065      0,1422      0,4363      0,6789 |
        |   0,4843   -0,0498     -0,0146      0,0104     -0,8733 |
        |   0,4843   -0,4052      0,3413     -0,2734      0,6403 |
        |   0,4843    0,1879     -0,3663     -0,5529      0,5388 |

Conceptually, such a factor of equal loadings could be interpreted as a "bias-factor" when an itemscale is evaluated. It says something like in which extension people tend to answer all questions/items with the same value

[update]

If we proceed as indicated in @amoeba's answer to append a pca-rotation on the remaining factors/axes (2..5) then this gives

load = rot(load,"-medimax")
load = rot(load,"pca",1..5,2..5)  
    \\ using items 1..5 for the computation of the rotation criterion and
    \\ using only factors 2..5 for the actual rotations
        |   0,4843   -0,5579      0,3475     -0,5423      0,1983 |
        |   0,4843   -0,6450      0,4441      0,2849     -0,2664 |
        |   0,4843    0,8657      0,0039     -0,1031      0,0739 |
        |   0,4843   -0,6411     -0,2687      0,2481      0,4698 |
        |   0,4843   -0,5530     -0,5510     -0,2341     -0,3183 |

The covariances of the factors are then in the following matrix:

chk = load' * load
              f1       f2            f3          f4          f5 
 f1     |   1,1727   -0,7416     -0,0117     -0,1678      0,0763 |
 f2     |  -0,7416    2,1935      0,0000      0,0000     -0,0000 |
 f3     |  -0,0117    0,0000      0,6938     -0,0000     -0,0000 |
 f4     |  -0,1678    0,0000     -0,0000      0,5023      0,0000 |
 f5     |   0,0763   -0,0000     -0,0000      0,0000      0,4378 |

---

Here is pseudocode how to compute the rotation, this should easily be transferable to MatLab etc.

The final rotation is done by iterating the following scheme for one pass until convergence.

One pass is: for all pairs of factors $(F_i,F_j) , j>i$ copy $F_i \to X $ and $F_j \to Y$ and compute a rotation-criterion from the (selected) rows of invidual terms $x_r,y_r$ of the vectors $X,Y$ which are the current loadings in this factors. From this the rotation-angle $\varphi$ and from this a pair of $cp=\cos(\varphi),sp=\sin(\varphi)$ is computed.

Then use that $cp,sp$ to actually rotate $X,Y \to F_i',F_j'$.

rot(F, "medimax")
    // rotate the complete matrix F ,take COLUMNS as axes
rot(F, "medimax", itemlist , factorlist ) 
    // only for selected items in itemlist compute the criterion
      // examples itemlist=1..4  or itemlist =2..5´9..11
    // rotate only in the subspace of factorlist 
      // examples factorlist=1..4  or factorlist =2..5´9..11

Pseudoalgorithm (1 pass of Iterations):

for all pairs of (the selected) factors F_i und F_j with j>i 
  copy columnvector F_i into columnvector X
  copy columnvector F_j into columnvector Y

  sx,sy,sx2,sy2 :=0
  for all rows/items { of the itemlist}
        x2 := x² - y²                      y2 := 2 * x * y 
        sx := sx + x                       sy := sy + y
        sx2 := sx2 + x2                    sy2 := sy2 + y2 
  end for       

  sx2 := sx2 - (sx² - sy²)            sy2:= sy2 - (2 * sx * sy)

  len := sqrt ( sx2² + sy2²)
  c2p := s2x / len                 s2p := s2y / len
  cp  := sqrt ( 0.5(1 + c2p))      sp  := sqrt ( 0.5(1 - c2p))
    if sy<0 then sp = - sp 

  // --------------------------------------------------------
  for all items/rows { rotate the whole system in that axes F_i,F_j}
     let r denote the row-index:
     F_i,r := cp * X_r - sp * Y_r         
     F_j,r := cp * Y_r + sp * X_r 
  end for

  end pair F_i,F_j
  end pass