Skip to content

Instantly share code, notes, and snippets.

Created April 22, 2015 21:59
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
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
#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)
#Scale the matrix A
function scale_matrix(A)
return scaled_ones(size(A,1)) * A
function rsq(A,B)
tbatab = transpose(B) * A * transpose(A) * B
return sum(sqrt(abs(eigvals(tbatab))))
function proc_reflect(A,B)
Abar = scale_matrix(A)
Bbar = scale_matrix(B)
return rsq(Abar,Bbar) / sum(diag(transpose(Bbar) * Bbar))
function proc_rotate_reflect(A,B)
X = transpose(scale_matrix(A)) * scale_matrix(B)
Xsvd = svdfact!(X)
return transpose(Xsvd.Vt) * transpose(Xsvd.U)
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))
#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