Skip to content

Instantly share code, notes, and snippets.

@gvdr
Created Apr 22, 2015
Embed
What would you like to do?
Compute the ordinary Procrustes sum of squares of two matrices in Julialang. All credits to procOPA {shapes}
#This is just a quick translation of an ordinary Procrustes sum of squares method, as done in the R function procOPA.
#All credits to Dryden, I. L. (2014). shapes package. R Foundation for Statistical Computing, Vienna, Austria. Contributed package. Version 1-1.10. URL http://www.R-project.org
#The function are, for now, as they are. I'll be working on them in the future...
#Compute a matrix with ones on the main diagonal
#and -1/n elsewhere
function scaled_ones(n)
eye(n) - (1/n) * ones(n,n)
end
#Scale the matrix A
function scale_matrix(A)
return scaled_ones(size(A,1)) * A
end
function rsq(A,B)
tbatab = transpose(B) * A * transpose(A) * B
return sum(sqrt(abs(eigvals(tbatab))))
end
function proc_reflect(A,B)
Abar = scale_matrix(A)
Bbar = scale_matrix(B)
return rsq(Abar,Bbar) / sum(diag(transpose(Bbar) * Bbar))
end
function proc_rotate_reflect(A,B)
X = transpose(scale_matrix(A)) * scale_matrix(B)
Xsvd = svdfact!(X)
return transpose(Xsvd.Vt) * transpose(Xsvd.U)
end
function proc_OSS(A,B)
R = proc_rotate_reflect(A,B)
s = proc_reflect(A,B)
Ahat = scale_matrix(A)
Bhat = scale_matrix(B) * R * s
resids = Ahat - Bhat
return sum(diag(transpose(resids) * resids))
end
#P.S. this is a clone of an anonymous gist I mistakenly published. I'll clean up stuff as soon as github's fellows allow me to do it.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment