Speed up sparse matrix calculations

Is it possible to speed up large sparse matrix calculations by e.g. placing parantheses optimally?

What I'm asking is: Can I speed up the following code by forcing Matlab to do the operations in a specified order (for instance "from right to left" or something similar)?

I have a sparse square symmetric matrix H (4000x4000), that previously has been factorized, and a sparse vector M with length equal to the dimension of H. What I want to do is the following:

k = 0;
cList = 1:3500     % List of contingecies to study
nc = length(cList);
nb = size(bus,1); % bus is a matrix with information about buses in the system
resVa = zeros(nb,nc);
[L U P] = lu(H);      % H is sparse (thus also L, U and P). Dimension = nbxnb

while k < nc
  k = k + 1;
  cont = cList(k);

    bf = branch(cont,1);  % branch is a matrix with information about branches
    bt = branch(cont,2);  % Col 1 = from-bus, col 2 = to-bus

    dy = -H(bf,bt);      % 
    M = sparse([bf,bt],1,[1,-1],nb,1);

    z = -M'*(U \ (L \ (P * M)));
    c = (1/dy + z)^-1;

    % V = Vm*exp(j*Va); - A complex vector of dimension nbx1
    % Sbus is a complex vector of dimension nbx1
    % Ybus is a complec matrix of dimension nbxnb
    mis = (V .*conj((Ybus)*V) - Sbus) ./ Vm;
    P_mis = real(mis);  

    converged = 0;
    i = 0;

    while (~converged && i < max_it)
       i = i + 1;
       %%The lines I hope to optimize:
       dVa = - (U \ (L \ (P * P_mis))); 
       dVaComp = (U \ (L \ (P * M * c * M' * dVa)));
       Va = Va + dVa + dVaComp;

       V = Vm .* exp(1j * Va);
       mis = (V .*conj((Ybus)*V) - Sbus) ./ Vm;    
       P_mis = real(mis);

       normP = norm(P, inf);
       if normP < tol
         converged = 1;
         break;
       end
    end
  resVa(:,k) = Va * 180 / pi;
end

Some additional information regarding the matrices:

All diagonal elements of H are non-zeros (it's still square, sparse and symmetrical).

Ybus and H have the same structure, but Ybus is complex and H is real.

Vm is updated using the imaginary part of mis, and some other matrices, but I've left it out for simplicity.

EDIT:

I achieved almost 75% reduction in computation time by using the extended syntax for lu when factorizing H.

[L U P Q] = lu(H);

Please let me know if any additional savings are possible to achieve.

Thanks!

[Cedric]

What size is H typically?

[John Doe]

Typically around 4000x4000. Each column contains about 4 non-zero elements.

The code above (except the lu-factorization) is in a loop that runs 4000 times. The two lines in the end ( dVa, and dVaComp) run ten times for each of those 4000 rounds, thus a total of 40000 times (P is updated within that loop). I added some information in the question now.

[Jan]

@Robert P: It would be much more convenient and useful to post the Matlab code instead of a mathematical notation. Instead of optimizing your code, we have to write it from scratch. On one hand this is more work for us, on the other hand the created code will not necessarily match your one sufficiently. So please post the code, which you want to get optimized.

You mathematical notation is even not clear:

start loop # 2 (10 iterations)
  dVa = - (U \ (L \ (P * P))); 
  dVaComp = (U \ (L \ (P * M * c * M' * dVa)));
end loop # 2

It is not specified, what is changed in the iterations. Then it is not meaningful to run the loop 10 times at all.

Calculating the inverse explicitly is not useful in general: it is slower and less accurate compared to solving the corresponding linear equation system.

[John Doe]

Thanks for the hints Jan! I have updated the question, please let me know if I should do further changes, or provide additional information.