Skip to content

Instantly share code, notes, and snippets.

Created Apr 22, 2015
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