[julia-users] Matrix crossproduct by groups of rows

[Giuseppe Ragusa]

I have been struggling to find a fast and elegant way to implement the 
following algorithm.

I have an Array{Float64, 2}, say X = randn(100,2) and an Array{Int64, 1}, 
call it cl which denotes row-blocks of X. More concretely:

  0.863377   0.867817   1.0
 -0.310559  -0.863393   1.0
  1.74963    0.0464976  1.0
 -1.29083   -0.496978   1.0
  ⋮                        
  1.04181    0.682401   5.0
  0.427546   0.237533   5.0
  0.650279  -0.727299   5.0
  0.121416   1.61274    5.0


What I have to calculate is the crossproducts X_1'X_1,  X_2'*X_2, ..., 
X_N'X_N, where X_i = X[find(cl.==i),:]. 

Something like this work, but I feel is way from optimal:

X = randn(100,2)
rep(v, n) = [v[div(i,n)+1] for i=0:n*length(v)-1]
cl=rep(1:5, 20)
out = zeros(2,2)


function block_crossprod!(out, X, cl)
    for j in distinct(cl)
       out += X[find(cl.==j),:]'*X[find(cl.==j),:]
    end
end

Any suggestions will be greatly appreciated.

Thank you.