This file contains the Matlab code that appears in the paper:

Gradient Projection Algorithms and Software for Arbitrary
Rotation Criteria in Factor Analysis.

by

Coen A. Bernaards and Robert I. Jennrich

Website: http://www.stat.ucla.edu/research

--------------------------------------------------------------
What follows is the general GPF subroutine for orthogonal
rotation. This is use without change.
--------------------------------------------------------------

function [Lh,Th,table]=GPForth(A,T)

al=1;
L=A*T;
[f,Gq]=vgQ(L);
G=A'*Gq;
table=[]
for iter=0:500
    M=T'*G;
    S=(M+M')/2;
    Gp=G-T*S;
    s=norm(Gp,'fro');
    table=[table;iter f log10(s) al];
    if s<10^(-5),break,end

    al=2*al;
    for i=0:10
       X=T-al*Gp;
       [U,D,V]=svd(X,0);
       Tt=U*V';
       L=A*Tt;
       [ft,Gq]=vgQ(L);
       if ft<f-.5*s^2*al,break,end
       al=al/2;
       end
    T=Tt;
    f=ft;
    G=A'*Gq;
    end
    Th=T;
    Lh=A*T;

--------------------------------------------------------------
What follows is  a vgQ subroutine defining the quartimax
criterion and gradient.
--------------------------------------------------------------

function [q,Gq]=vgQ(L)

q=-norm(L.^2,'fro')^2/4;
Gq=-L.^3;

--------------------------------------------------------------
What follows is the the main program for quartimax rotation
of Harman's 8 physical variables.
--------------------------------------------------------------

A=[ .830 -.396
     .818 -.469
     .777 -.470
     .798 -.401
     .786  .500
     .672  .458
     .594  .444
     .647  .333];

T=eye(2);
[L,T,table]=GPForth(A,T);

table,L

--------------------------------------------------------------
Turning to the oblique case what follows is the general GPF
subrutine for oblique rotation. This is used without change.
--------------------------------------------------------------

function [Lh,Phi,Th,table]=GPFoblq(A,T)

al=1;
table=[];
Ti=inv(T);
L=A*Ti';
[f,Gq]=vgQ(L);
G=-(L'*Gq*Ti)';
for iter=0:500
    Gp=G-T*diag(sum(T.*G));
    s=norm(Gp,'fro');
    table=[table;iter f log10(s) al];
    if s<10^(-5),break,end

    al=2*al;
    for i=0:10
       X=T-al*Gp;
       v=1./sqrt(sum(X.^2));
       Tt=X*diag(v);
       Ti=inv(Tt);
       L=A*Ti';
       [ft,Gq]=vgQ(L);
       if ft<f-.5*s^2*al,break,end
       al=al/2;
       end
    T=Tt;
    f=ft;
    G=-(L'*Gq*Ti)';
    end
Th=T;
Lh=L;
Phi=T'*T;

--------------------------------------------------------------
What follows is the vgQ subroutine defining the quartimin
criterion and its gradient.
--------------------------------------------------------------

function [q,Gq]=vgQ(L)

L2=L.^2;
[p,k]=size(L);
N=ones(k,k)-eye(k);
X=L2*N;
q=sum(sum(L2.*X))/4;
Gq=L.*X;

--------------------------------------------------------------
What follows is the main program for quartimin rotation of
Harman's 8 physical variables.
--------------------------------------------------------------

A=[ .830 -.396
     .818 -.469
     .777 -.470
     .798 -.401
     .786  .500
     .672  .458
     .594  .444
     .647  .333];
T=eye(2);

[L,phi,T,table]=GPFoblq(A,T);

table,L,phi

--------------------------------------------------------------
The second half of the Bernaards and Jennrich paper describes
how to produce vgQ subroutines foe a variety of rotation
methods. Matlab code for these are given here. We begin with
a main program for testing these using Thurstone's box
problem. Included is code for making sorted absolute loading
plots. After the vgQ subroutines there is a program for
promax rotation.
--------------------------------------------------------------

A=[
.659 -.736  .138
.725  .180 -.656
.665  .537  .500
.869 -.209 -.443
.834  .182  .508
.836  .519  .152
.856 -.452 -.269
.848 -.426  .320
.861  .416 -.299
.880 -.341 -.354
.889 -.417  .436
.875  .485 -.093
.667 -.725  .109
.717  .246 -.619
.634  .501  .522
.936  .257  .165
.966 -.239 -.083
.625 -.720  .166
.702  .112 -.650
.664  .536  .488
];

T=eye(3);

[L,phi,T,table]=GPFoblq(A,T);

table,L

[p,k]=size(L);
sal=sort(abs(reshape(L,p*k,1)));
plot(sal,'o')

--------------------------------------------------------------
oblimin family
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

gm=0;

L2=L.^2;
[p,k]=size(L);
N=ones(k,k)-eye(k);
C=ones(p,p)/p;
X=(eye(p)-gm*C)*L2*N;
q=sum(sum(L2.*X))/4;
Gq=L.*X;

--------------------------------------------------------------
Crawford-Ferguson family.
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

kp=0;

L2=L.^2;
[p,k]=size(L);
M=ones(p,p)-eye(p);
N=ones(k,k)-eye(k);

X=(1-kp)*L2*N+kp*M*L2;
q=sum(sum(L2.*X))/4;
Gq=L.*X;

--------------------------------------------------------------
minimum entropy
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

L2=L.^2;
q=-sum(sum(L2.*log(L2)))/2;
Gq=-L.*log(L2)-L;

--------------------------------------------------------------
Bentler's invariant pattern simplicity criterion.
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

L2=L.^2;
X=L2'*L2;
D=diag(diag(X));
q=-(log(det(X))-log(det(D)));
Gq=-L.*(L2*(inv(X)-inv(D)));

--------------------------------------------------------------
Comrey's tandem II criterion
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

L2=L.^2;
[p,k]=size(L);
X0=ones(p,p);
X1=L*L';
X2=L2*L2';
q=sum(sum(L2.*((X0-X1.^2)*L2)));
Gq=4*L.*((X0-X1.^2)*L2)-4*(X1.*X2)*L;

--------------------------------------------------------------
geomin
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

ep=.01;

L2=L.^2;
[p,k]=size(L);
u=ones(k,1);
e=exp(log(L2+ep)*u/k);
q=sum(e);
Gq=(2/k)*L.*(e*u')./(L2+ep);

--------------------------------------------------------------
target rotation
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

global H

q=norm(L-H,'fro')^2;
Gq=2*(L-H);

--------------------------------------------------------------
partially specified target rotation
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

global H W

q=norm(W.*(L-H),'fro')^2;
Gq=2*W.*(L-H);

--------------------------------------------------------------
simplimax
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

global m;

L2=L.^2;
[p,k]=size(L);
ls=reshape(L2,p*k,1);
ls=sort(ls);
lsm=ls(m);
W=(L2<=lsm);
q=sum(sum(L2.*W))/2;
Gq=L.*W;

--------------------------------------------------------------
infomax
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

S=L.^2;

[p,k]=size(L);
u=ones(p,1);
v=ones(k,1);
s1=S*v;
s2=u'*S;
s=s2*v;

P=S/s;
p1=s1/s;
p2=s2/s;

Q0=-sum(sum(P.*log(P)));
Q1=-sum(p1.*log(p1));
Q2=-sum(p2.*log(p2));

q=log(k)+Q0-Q1-Q2;

H=-(log(P)+1);
al=u'*(S.*H)*v/s^2;
G0=H/s-al*u*v';

h1=-(log(p1)+1);
al1=s1'*h1/s^2;
G1=h1*v'/s-al1*u*v';

h2=-(log(p2)+1);
al2=h2*s2'/s^2;
G2=u*h2/s-al2*u*v';

Gq=(2*L).*(G0-G1-G2);

--------------------------------------------------------------
McCammon minimum entropy ratio
--------------------------------------------------------------
function [q,Gq]=vgQ(L)

[p,k]=size(L);
S=L.^2;

u=ones(p,1);
v=ones(k,1);
M=u*u';
R=M*S;
P=S./R;
H=-(log(P)+1);

Q1=-sum(sum(P.*log(P)));
G1=H./R-M*(S.*H./R.^2);

s2=sum(S);
s=sum(s2);
p2=s2/s;
h=-(log(p2)+1);
al=h*p2';

Q2=-sum(p2.*log(p2));
G2=u*h/s-al*u*v';

q=log(Q1)-log(Q2);
Gq=2*L.*(G1/Q1-G2/Q2);

--------------------------------------------------------------
Promax rotation. What follows is not a vgQ subroutine, but
rather a subroutine for procrustes rotation. This form of
rotation is not based on optimizing a criterion and hence
requires a special subroutine. The procrustes subroutine may
be used for promax rotation by choosing an appropriate target
H, for example the element-wise cube of a varimax rotation of
A.
--------------------------------------------------------------
function [L,T]=procrustes(H,A)

S=inv(A'*A)*A'*H;
d=sqrt(diag(inv(S'*S)));
D=diag(d);
L=A*S*D;
T=inv(S*D)';