Talk:Varimax rotation

R

FWIW, the R code for varimax appears to be the following. —Ben FrantzDale (talk) 02:40, 1 August 2008 (UTC)[reply]

> varimax
function (x, normalize = TRUE, eps = 1e-05) 
{
    nc <- ncol(x)
    if (nc < 2) 
        return(x)
    if (normalize) {
        sc <- sqrt(drop(apply(x, 1, function(x) sum(x^2))))
        x <- x/sc
    }
    p <- nrow(x)
    TT <- diag(nc)
    d <- 0
    for (i in 1:1000) {
        z <- x %*% TT  # Matrix multiplication
        B <- t(x) %*% (z^3 - z %*% diag(drop(rep(1, p) %*% z^2))/p)
        sB <- La.svd(B)
        TT <- sB$u %*% sB$vt 
        dpast <- d
        d <- sum(sB$d)
        if (d < dpast * (1 + eps)) 
            break
    }
    z <- x %*% TT
    if (normalize) 
        z <- z * sc
    dimnames(z) <- dimnames(x)
    class(z) <- "loadings"
    list(loadings = z, rotmat = TT)
}

Other algorithms

I'm looking at this... it looks like the goal is, given a subspace defined by a subset of the eigenvectors, find a rotation that maximizes the variance of the terms of the eigenvectors. This sounds a lot like rotation recovery using SVD...

Since all eigenvectors have the same length, the variance is maximized when n times the variance is maximized. If v is an n-by-d matrix with eigenvectors as its columns, and with a column for each of the d dimensions we are interested in, n times the variance of the terms of each is found on the diagonal of the matrix $v^{\top }v$ .

In rotation recovery, you find a rotation between point clouds by subtracting the mean from both and then taking the outer product of the two arrays. That is, $x^{\top }y$ is the outer product between them, and you want to find the rotation that registers them. It turns out that the rotation you want is the one that maximizes the trace of this outer product. You can find R in SO(3) to solve

\operatorname {argmax} _{R\in SO(3)}\operatorname {trace} (x^{\top }Ry)

If you take the svd to find

USV^{\top }=x^{\top }y

,

then the (possibly-improper) rotation is

R=UV^{\top }

.

This is related to what the R code is doing, but not quite... it's making TT be the rotation matrix and z be the new basis. That is, it finds R and z such that

z=xR

.

The iteration does

z:=xR

B_{ij}:=x_{ik}(z_{kj}^{3}-z_{k\ell }\delta _{\ell m}z_{mj}^{2}/ndof)=x_{ik}(z_{kj}^{3}-z_{k\ell }z_{\ell j}^{2}/ndof)

[U,S,V]:=\operatorname {SVD} (B)

R:=UV^{\top }

d:=\operatorname {trace} (B)

At this time of day, it's not clear to me why this works, and if so, why it needs to iterate. —Ben FrantzDale (talk) 04:51, 28 July 2009 (UTC)[reply]

So yea, the references are generally vague. I added a good one. VARIMAX tries to minimize the variance of the squared components, not the variance of the components. Here's a Python version:

def varimax(Phi, gamma = 1.0, q = 20, tol = 1e-6):
    from scipy import eye, asarray, dot, sum
    from scipy.linalg import svd
    p,k = Phi.shape
    R = eye(k)
    d=0
    for i in range(q):
        d_old = d
        Lambda = dot(Phi, R)
        u,s,vh = svd(dot(Phi.T,asarray(Lambda)**3 - (gamma/p) * dot(Lambda, diag(diag(dot(Lambda.T,Lambda))))))
        R = dot(u,vh)
        d = sum(s)
        if d_old!=0 and d/d_old < 1 + tol: break
    return dot(Phi, R)

Have fun. —Ben FrantzDale (talk) 05:41, 4 August 2009 (UTC)[reply]