optimize sequence of matlab operations for solving a system of linear equations

I need to apply some sequence of operations with big matrices in matlab in order to solve a system of linear equations. In the beginning I though that this would be quite fast considering that matlab applies some internal optimization (at least that was my understanding and what I was reading in the net). However, it seems that I was wrong. The initial sequence of operations that I am applying are the following:

nbf = size(F,1);

if issparse(F)
   F = sparse(diag(r))*F;
else
   F = diag(r)*F;
end
if issparse(F)
   K = sparse(eye(nbf)) - F;
else
   K = eye(nbf) - F;
end

R = K\E;

where F is a square matrix MxM, where M = 50000 or more and almost all the times is not sparse. r and E are Mx1 vectors, and E is sparse. These operations on a powerful machine take around 25 minutes which is quite slow for my needs(in the ideal case I would like to have some real time computation but I do not think I can achieve that, I would be happy though if someone proves me wrong). Thus, I am trying to find some ways to reduce this computation time. A first idea was to see what I can get if I replace the F = diag(r)*F; and K = eye(nbf) - F; with F = bsxfun(@mtimes,diag(r),F); and K = bsxfun(@minus,eye(nbf),F); respectively. Regarding solving the inverse matrix operation R = K\E I came upon this https://it.mathworks.com/help/matlab/ref/decomposition.html or if you are not able to use the decomposition() function you could replace it with the corresponding factorize() from here , thus it could be replaced with R = decomposition(K)\E. However, this improved the computation time just by little 22 mins instead of 25. Therefore, I would like to ask if someone knows any suggestions or other ways how I could optimize these operations for computation speed.

Answer 1

[Christine Tobler]

I assume that most of the time is spent in the line R = K\E (as this has cubic complexity in the size of A, while the other operations have square complexity). So any changes to the previous lines will not have much influence. You could use the MATLAB Profiler to verify this.

I'm afraid there's not much to be done about the speed of backslash for a general dense matrix of this size, unless the matrix K has some special structure, or if you are reusing the same matrix K multiple times.

Is E a vector, or a matrix with many columns? How is R used after computing it?

[ThT]

Hi Christine, thank you for the response. Actually from what I've noticed (and measured) the most of the time is consumed in this operation F = diag(rho)*F; This really takes long time, thus I was trying to see if I can do something about it. The issue is that my F matrix is gonna be always more than 50000x50000 elements and in some cases it even reach up to 150000x150000.

E is a vector 1x50000 or 1xM, where M is the one of the sizes of matrix F. R is the result I need.

The actual question is can such a multiplication with so many elements go down to real time calculation (by using any mex libraries or something) or I am just chasing something impossible.

[Christine Tobler]

Hi ThT,

The line F = diag(rho)*F is quite inefficient, since it constructs the diagonal matrix diag(rho), only to discard it again right away. As you mentioned in the initial post, this should be replaced by F = bsxfun(@times, rho(:), F) (or, in newer MATLAB versions, simply by F = rho(:).*F). Similarly, it's better to always construct large identity matrices as sparse matrices, even if they are subtracted from a dense matrix.

Here's how I would reformulate the code above:

K = speye(nbf) - r(:).*F;
R = K\E;

For an example matrix of size 10000-by-10000, the first line took 0.75 seconds on my machine, while the second line took 3.660881 seconds.