How to speed up function evaluation by vectorization

Is there a way to speed up this code significantly? Maybe by more vectorization?

It is the evaluation of a Hesse matrix. The first loop is a "dot product" between basis elements in the form of 16x16 matrices (Lambda) and corresponding coefficients (T). The second (double-) loop is the evaluation of the Hesse matrix itself.

This is a more generic version of the code:

 Hesse=eye(256,256);
 B=rand(1296,256);
 T=rand(256,1);
 Lambda=rand(16,16,256);
 f=rand(1296,1);
 P_int=B*T;

 rho_int=zeros(16);
 %Loop 1
 for j=1:256
    rho_int=rho_int+T(j,1)*Lambda(:,:,j)/2^4;
 end

 rho_int_inv=inv(rho_int);

 %Loop 2
 for k=2:256
    for l=2:256
        Hesse(k,l)=dot(f, (B(:,k).*B(:,l))./(P_int.^2))+trace(rho_int_inv*Lambda(:,:,k)*rho_int_inv*Lambda(:,:,l)/4^4);
    end
 end

[Andrei Bobrov]

P_int - what is it?

[Frederik Lohof]

In my setting P_int is a vector of probabilities that is used in the calculation of the Hesse matrix. The original function to the Hesse matrix is a likelihood function which is obtained in the context of quantum state tomography (quantum state estimation). The likelihood function is a measure for the likelihood of a certain quantum state (parametrized by T) given the frequencies f of certain measurement outcomes. These are obtained (in my case simulated) by measuring n identical copies of the same (unknown) quantum state in different measurement bases. P_int is a vector of probabilities of measuring a certain outcome given a certain state (B is a matrix yielding the connection between the parametrization T and the expected probabilities P_int).

Edit: And I made a typo in the initialization: It is P_int=B*T :)

Answer 1

[Kirby Fears]

The second loop seems fine. Vectorization of loop 1 below.

 % rho_int=zeros(16);
 % Loop 1 replaced by vector operation
 rho_int=sum(repmat(reshape(T,1,1,numel(T)),...
     size(Lambda,1),size(Lambda,2)).*Lambda,3)/2^4;
 rho_int_inv=inv(rho_int);

You could speed up the second loop with the parallel computing toolbox.

Hope this helps.

[Frederik Lohof]

Thanks, I will look at that.