因子旋转的理论基础

写本博客的目的是：在使用varimax方法时不知道原理到底是什么？查了好多教材和网上的资料，都没有具体的理论推导，无法真正理解此方法的内核计算。最后终于在factor analysis as a statistical method书中的chapter6找到相关内容，并进行整理，最后给出相关matlab代码。

1. varimax 方法

1.1 varimax的理论基础

在varimax旋转方法出现之前，通常使用旋转的图形方法，目的是消除负的因子载荷(loading)并且尽可能少的非0因子载荷来描述数据。但这种方法具有一定的主观性。

varimax方法与图形方法在某种程度上具有相似的特点。 与原始的因子载荷相比，用这种方法对因子进行旋转得到的新的因子载荷要么相对大要么相对小。正如我们所看到的，这种旋转的方法是通过最大化loading的平方的一个特定函数并且迭代完成的。在Kaiser’s的原始旋转方法中，因子是成对地旋转（两个两个的旋转）直到loadings收敛到最终值。下面我们描述的是"The simultaneously factor varimax solution"（同时因子varimax解），显然，这个方法的最主要特点是在每一次迭代时，所有因子同时被旋转，并非Kaiser’s的成对成对地旋转。它比原始的varimax方法更快，并且没有累积的舍入误差(因为原始方法是成对旋转的)。

将 矩阵${\mathbf{\Lambda }}_0(p\times k)$ 看做是未旋转的因子载荷，则${\mathbf{\Lambda }}_0^{\prime }$的第$i$列（即${\mathbf{\Lambda }}_0$ 的第$i$行）是一个$k$维的列向量，用${\mathbf{I}}_i$表示。 定义$\mathbf{M}$是一个$k\times k$的正交旋转矩阵，它的第$r$th列用${\mathbf{m}}_r(r=1,\cdots ,k)$。新的旋转的因子载荷矩阵$\mathbf{\Lambda }=\left[{\lambda }_{ir}\right]$，可以通过如下计算：
$$\mathbf{\Lambda }={\mathbf{\Lambda }}_0\mathbf{M}$$
因此，
$${\lambda }_{ir}={\mathbf{l}}_i^{\prime }{\mathbf{m}}_r.$$

定义一个标量$d_r(r=1,\cdots ,k)$：

$$d_r=\sum _{i=1}^p{\lambda }_{ir}^2=\sum _i{\left({\mathbf{I}}_i^{\prime }{\mathbf{m}}_r\right)}^2$$

$d_r$表示的是$\mathbf{\Lambda }$的第$r$th列的因子载荷的平方和。

在同时varimax方法中的标准$x$被最大化，形式如下：

$$\begin{array}{rl}x& =\sum _{r=1}^k\left[\sum _{i=1}^p{\left({\lambda }_{ir}^2-d_r/p\right)}^2\right]\\ & =\sum _r\left[\sum _i{\lambda }_{ir}^4-d_r^2/p\right]\\ & =\sum _r\left[\sum _i{\left({\mathbf{I}}_i^{\prime }{\mathbf{m}}_r\right)}^4-d_r^2/p\right]\end{array}$$
从这个定义中，我们可以看出：$x$代表的是${\lambda }_{ir}^2$与其对应的列均值的偏差平方和。这个$x$的最大化与$\mathbf{M}$有关。因为$\mathbf{M}$是正交阵，所以它的列满足如下条件：

$$\begin{array}{rlr}{\mathbf{m}}_r^{\prime }{\mathbf{m}}_r& =1,& \\ {\mathbf{m}}_r^{\prime }{\mathbf{m}}_s& =0& (r\ne s)\end{array}$$

使用拉格朗日乘子的方法，表达式如下：

$$y=x-2\sum _x\sum _s{a}_{rs}{\mathbf{m}}_r^{\prime }{\mathbf{m}}_s$$

$a$系数代表的是不确定的乘子，且$a_{rs}=a_{sr}$.

对$y$关于${\mathbf{m}}_s$求导并令其等于0，求解如下：

$$\begin{array}{rl}\partial y/\partial {\mathbf{m}}_s& =4\sum _i\left[{\left({\mathbf{I}}_i^{\prime }{\mathbf{m}}_s\right)}^3{\mathbf{I}}_i\right]-\left(4d_s/p\right)\sum _i\left[\left({\mathrm{I}}_i^{\prime }{\mathbf{m}}_s\right){\mathrm{I}}_i\right]-4\sum _r\left(a_{rs}{\mathbf{m}}_r\right)\\ & =4\sum _i\left(c_{is}-d_s{\lambda }_{is}/p\right){\mathbf{I}}_i-4\sum _r{a}_{rs}{\mathbf{m}}_r\end{array}$$
其中，

$$c_{ir}={\left({\mathbf{I}}_i^{\prime }{\mathbf{m}}_r\right)}^3={\lambda }_{ir}^3.$$

经检查发现：$\partial y/\partial {\mathbf{m}}_s$是下面矩阵的第$s$th列：

$$4\left[{\mathbf{\Lambda }}_0^{\prime }\mathbf{C}-(1/p){\mathbf{\Lambda }}_0^{\prime }\mathbf{\Lambda }\mathbf{D}-\mathbf{M}\mathbf{A}\right]$$
其中，$\mathbf{C}$是一个$p\times k$矩阵，其元素为$c_{ir}$，$\mathbf{A}$是$k\times k$的对称矩阵，且元素为$a_{rs}$，
$\mathbf{D}$是对角矩阵，对角元为$d_1,d_2,,\cdots ,d_k$.因此将$s$的所有值都考虑进去，有：

$$\begin{array}{l}\partial y/\partial \mathbf{M}=4(\mathbf{B}-\mathbf{M}\mathbf{A})\end{array}$$
其中，
$$\begin{array}{l}\mathbf{B}={\mathbf{\Lambda }}_0^{\prime }[\mathbf{C}-(1/p)\mathbf{\Lambda }\mathbf{D}]\end{array}$$

因此，最大化$x$的正交矩阵$\mathbf{M}$满足如下等式：

$$\mathbf{M}\mathbf{A}=\mathbf{B}$$

可以证明：如果$\mathbf{A}$是正定的，由于$\mathbf{A}$也是对称的，那么上述的等式对应于$x$的一个最大值。对上式左乘一个${\mathbf{M}}^{\prime }$，有：

$$\begin{array}{rl}\mathbf{A}& ={\mathbf{M}}^{\prime }\mathbf{B}\\ & ={\mathbf{\Lambda }}^{\prime }\mathbf{C}-(1/p){\mathbf{\Lambda }}^{\prime }\mathbf{\Lambda }\mathbf{D}\end{array}$$
因此，当$x$被最大化时，等式右边应该是对称矩阵。$\mathbf{A}$的第$r$th个对角元为：

$$\sum _i{\lambda }_{ir}^4-(1/p){\left(\sum _i{\lambda }_{ir}^2\right)}^2$$

基于$x$的定义，我们发现：$\mathrm{tr}(\mathbf{A})$就是$x$的最大值。

满足$\mathbf{M}\mathbf{A}=\mathbf{B}$的$\mathbf{M}$,$\mathbf{A}$和$\mathbf{B}$可以通过迭代获得。我们首先用一个近似的初始矩阵${\mathbf{M}}_1$,这将产生一个$\mathbf{\Lambda }$初始的近似${\mathbf{\Lambda }}_1={\mathbf{\Lambda }}_0{\mathbf{M}}_1$。并且，从${\mathbf{\Lambda }}_1$中，我们获得$\mathbf{C}$,$\mathbf{D}$和$\mathbf{B}$的近似${\mathbf{C}}_1$,${\mathbf{D}}_1$和${\mathbf{B}}_1$,x现在我们来求对称和正定矩阵${\mathbf{A}}_1$和第二个正交旋转矩阵${\mathbf{M}}_2$:

$${\mathbf{M}}_2{\mathbf{A}}_1={\mathbf{B}}_1$$

注意对上述方程的等式两边分别左乘自己的转置矩阵，有如下等式：

$${\mathbf{A}}_1^2={\mathbf{B}}_1^{\prime }{\mathbf{B}}_1$$
现在，${\mathbf{B}}_1^{\prime }{\mathbf{B}}_1$是$k\times k$的对称和正定矩阵，因此我们可以用如下形式来表示：$\mathbf{U}{\Delta }^2{\mathbf{U}}^{\prime }$，其中$\mathbf{U}$是正交矩阵，且$\Delta$是具有正的对角元素的对角矩阵。因此我们可以用以下形式来表示${\mathbf{A}}_1$:

$${\mathbf{A}}_1=\mathbf{U}\mathbf{\Delta }{\mathbf{U}}^{\prime }$$

则${\mathbf{M}}_2$为：

$${\mathbf{M}}_2={\mathbf{B}}_1{\mathbf{A}}_1^{-1}$$

然后用${\mathbf{M}}_2$替代${\mathbf{M}}_1$，在对上述程序进行重复，我们会得到一个${\mathbf{M}}_s$的矩阵序列。可以证明：当${\mathbf{M}}_1$与$\mathbf{M}$充分接近时，序列收敛到$\mathbf{M}$。且对应的矩阵序列${\mathbf{A}}_s$,它的$\mathrm{tr}({\mathbf{A}}_s)$是不断上升并且最终达到一个稳定值，这个稳定值便是标准$x$的最大值。

在实际中，选择一个好的$\mathbf{M}$的初始近似并不是必须的，通常情况下，${\mathbf{M}}_1=\mathbf{I}$ 并且${\mathbf{\Lambda }}_1={\mathbf{\Lambda }}_0$，通常只需要很少的几次迭代就可以了。

Horst建议，在最大化过程开始时，最好对原始的因子载荷矩阵${\mathbf{\Lambda }}_0$按行标准化，这样每一个行向量${\mathbf{I}}_i^{\prime }$就被尺度化为单位长度了。在最大化结束时，当获得$\mathbf{M}$后，最终的$\mathbf{\Lambda }$也需要重新按照行进行标准化，这样communalities的值才是正确的。

1.2 varimax的缺点

varimax有两个主要的缺点：

• 第一个缺点：如果在旋转过程中加入其他额外的因子，则通过varimax方法获得的因子载荷的pattern可能会发生很大的变化。
• 第二个缺点：当在因子中出现一个非常重要的一般因子（general factor）情况下，varimax方法可能并不work；

2. promax方法

在许多情况下，因子载荷的pattern可以通过进一步将因子转化为倾斜或者是相关的因子来进一步简化。可是这种简化需要付出一定的代价因为我们必须估计因子之间的相关系数,在对结果进行解释时也必须考虑这些因素。最终，人们的注意力转向了分析标准的发现，这些标准可以用来产生唯一的、不那么主观的解。有很多这样的方法，我们集中注意力介绍promax方法，因为它在实际中似乎效果更好。

2.1 promax方法的理论基础

promax方法的起始以varimax旋转的结果为起始。假设varimax方法产生的loadings是$\mathbf{\Lambda }=[{\lambda }_{ir}]$，并且构造另外一个$p\times k$的矩阵$\mathbf{Q}=\left[q_{ir}\right]$，它的元素按如下公式定义：

$$q_{ir}=\left|{\lambda }_{ir}^{m-1}\right|{\lambda }_{ir}$$

其中，$m$是满足$m>1$的某个整数。因此：如果${\lambda }_{ir}=0$则$q_{ir}=0$ ，否则，$q_{ir}$和${\lambda }_{ir}$有相同的符号和相同的绝对值。

现在，我们向寻找一个变换矩阵$\mathbf{U}$（此矩阵不必是正交的），对于$r=1,\dots ,k$，$\mathbf{\Lambda }\mathbf{U}$的第$r$列与$\mathbf{Q}$的第$r$列在最小二乘的意义上保持最大的一致性。一般来说，是当loadings相对较大时，$\mathbf{\Lambda }\mathbf{U}$的目的使其更大；当loadings相对较小时，$\mathbf{\Lambda }\mathbf{U}$使其更小。

$m$值的选择是一个反复试错的过程，但经验证明当简化loadings的pattern时，m大于4经常生成高度相关的因子，这在实际中是不被期望的。

将$\mathbf{Q}$ 和$\mathbf{U}$的第$r$列分别表示为${\mathbf{q}}_r$ 和${\mathbf{u}}_r$。然后，对于每一个$r$,我们选择的${\mathbf{u}}_r$使其最小化以下的表达式：
$${\left({\mathbf{q}}_r-\mathbf{\Lambda }{\mathbf{u}}_r\right)}^{\prime }\left({\mathbf{q}}_r-\mathbf{\Lambda }{\mathbf{u}}_r\right)$$
上式关于${\mathbf{u}}_{\mathbf{r}}$求导得到：

$$-2{\mathbf{\Lambda }}^{\prime }\left({\mathbf{q}}_r-\mathbf{\Lambda }{\mathbf{u}}_r\right)$$

令其等于0得到：
$$\left({\mathbf{\Lambda }}^{\prime }\mathbf{\Lambda }\right){\mathbf{u}}_r={\mathbf{\Lambda }}^{\prime }{\mathbf{q}}_r$$

因此，考虑所有的$r$，有：

$$\begin{array}{l}\left({\mathbf{\Lambda }}^{\prime }\mathbf{\Lambda }\right)\mathbf{U}={\mathbf{\Lambda }}^{\prime }\mathbf{Q}\\ \text{ or }\\ \mathbf{U}={\left({\mathbf{\Lambda }}^{\prime }\mathbf{\Lambda }\right)}^{-1}{\mathbf{\Lambda }}^{\prime }\mathbf{Q}\end{array}$$

事实上，重新尺度化$\mathbf{U}$的列是非常方便的，以至于转换后的因子具有单位方差。这可以通过找对角矩阵$\mathbf{D}$来实现，且它的对角元素是正的，且满足：
$${\mathbf{D}}^2=\mathrm{diag}\left[{\left({\mathbf{U}}^{\prime }\mathbf{U}\right)}^{-1}\right]$$

代替$\mathbf{U}$的变换矩阵为：
$$\mathbf{M}=\mathbf{U}\mathbf{D}$$
转换的因子载荷矩阵为：

$${\mathbf{\Lambda }}^{\ast }=\mathbf{\Lambda }\mathbf{M}$$

所以基于上述讨论有以下公式：

$$\begin{array}{rl}\mathbf{\Lambda }{\mathbf{\Lambda }}^{\prime }& ={\mathbf{\Lambda }}^{\ast }{\mathbf{M}}^{-1}{\mathbf{M}}^{\prime -1}{\mathbf{\Lambda }}^{\ast \prime }\\ & ={\mathbf{\Lambda }}^{\ast }{\left({\mathbf{M}}^{\prime }\mathbf{M}\right)}^{-1}{\mathbf{\Lambda }}^{\ast \prime }\\ & ={\mathbf{\Lambda }}^{\ast }\mathbf{\Phi }{\mathbf{\Lambda }}^{\ast \prime }\end{array}$$
其中，
$$\mathbf{\Phi }={\left({\mathbf{M}}^{\prime }\mathbf{M}\right)}^{-1}={\mathbf{D}}^{-1}{\left({\mathbf{U}}^{\prime }\mathbf{U}\right)}^{-1}{\mathbf{D}}^{-1}$$

$\mathbf{\Phi }$是变换矩阵的协方差矩阵。由$\mathbf{D}$的定义，很明显，$\mathbf{\Phi }$具有单位对角元素，因此$\mathbf{\Phi }$是新因子(被标准化后)的相关矩阵。

3. Transformation of correlated factors using a pattern matrix

4. MATLAB 实现

function [B,T] = rotateFactors(A,varargin)
if nargin > 1
    [varargin{:}] = convertStringsToChars(varargin{:});
end

names = {'method' 'normalize' 'reltol' 'maxit'};
dflts = {'varimax' [] [] []};
[method,normalize,reltol,maxit] = internal.stats.parseArgs(names, dflts, varargin{:});

methodNames = {'varimax'};
[method,~] = internal.stats.getParamVal(method,methodNames,'Method');

switch method
    case 'varimax'
        [B, T] = varimax(A, normalize, reltol, maxit);
end

end

%-----------------------------------------------------------------------

function [B, T] = varimax(A, normalize, reltol, maxit)
% VARIMAX method of orthogonal rotation of FA or PCA loadings.

[d, m] = size(A);

% Defaults to normalized varimax rotation
if nargin < 2 || isempty(normalize), normalize = 'on'; end
if nargin < 3 || isempty(reltol), reltol = sqrt(eps(class(A))); end
if nargin < 4 || isempty(maxit), maxit = 250; end

% Normalize the factor loadings
switch normalize
case {'on',1}
    h = sqrt(sum(A.^2, 2));
    A = bsxfun(@rdivide, A, h);
case {'off',0}
%     A = A;
otherwise
    error(message('stats:rotatefactors:BadNormalize'));
end

T = eye(m); % the intial rotation matrix is identity.
R = zeros(m); % the coefficient matrix

converged = false;
for n = 1:maxit
    Aold = A;
    A = Aold * T;
    D = sum(A.^2,1); % the sum of square of the column of B
    C = A.^3;
    P = Aold'*(C-(1/d).*A.*D);
    [U,L] = svd(P'*P);
    Rold = R;
    R = U.*sqrt(diag(L))'*U'; 
    T = P/R;
    
    if (trace(R)-trace(Rold)< reltol)
        B = A;
        converged = true;
        break;
    end
end  

if ~converged
    error(message('stats:rotatefactors:IterationLimit'));
end

% Unnormalize the rotated loadings
switch normalize
    case {'on',1}
        B = bsxfun(@times, B, h);
end
end