2

There are several ways in Matlab to calculate "LU decomposition". Here is one: 

function [L,A]=LU_factor(A,n)    
L=eye(n);
for k=1:n
    if (A(k,k) == 0) Error('Pivoting is needed!'); end
    L(k+1:n,k)=A(k+1:n,k)/A(k,k);
    for j=k+1:n
        A(j,:)=A(j,:)-L(j,k)*A(k,:);
    end
end

But my teacher told us that using for in MATLAB may decrease the efficiency of the program. He told us to calculate the LU decomposition using fewer fors. He said that you can find the required indexes without using for and then with some tricks you won't need to use for at all.

My first question is: will using for really decrease the speed of the program? My second question is: how can I store the required indexes in an array and use them instead of a for loop?

Divakar
  • 218,885
  • 19
  • 262
  • 358
bido
  • 47
  • 6
  • 1
    There is a function [lu](http://uk.mathworks.com/help/matlab/ref/lu.html) which is probably the most efficient way to do the decomposition in matlab :) – nivag Nov 11 '15 at 15:01
  • `for` loops used to be very inefficient in MATLAB and there is a lot of legacy dislike for them because of that. These days are they are actually far more efficient (thanks to the JIT compiler) provided you pre-allocate your matrices. Sometimes loops can even outperform vectorized solutions. But often they are still slower. The only real way to know is to write both out and time them yourself using something like the `timeit` function. – Dan Nov 11 '15 at 15:26
  • I'd argue the for loops here are perfectly fine. The FAST way to do an lu decomposition is the `lu` command, and there's no way to write something faster in MATLAB code. The point of writing your own `lu` function can't be for speed. The point has got to be to build conceptual understanding and practice. In coding outside of MATLAB, for loops are perfectly fine. Obscure, complicated code to avoid for loops is just MATLAB specific weirdness, it's not some generalizable skill. – Matthew Gunn Nov 11 '15 at 18:21

1 Answers1

2

There are two ways you can remove the inner loop.

If your teacher happen to be a bsxfun-lover -

A(k+1:n,:) = A(k+1:n,:) - bsxfun(@times,L(k+1:n,k),A(k,:))

Otherwise, 

A(k+1:n,:) = A(k+1:n,:) - L(k+1:n,k)*A(k,:)
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • Its not really **L**U decomposition if you don't create and return `L`. You might want to just stick with `A(k+1:n,:) = A(k+1:n,:) - L(k+1:n,k)*A(k,:)` to remove the inner loop. Also `bsxfun` is a little unnecessary here and I'd be extremely surprised if that's what his teacher meant. – IKavanagh Nov 11 '15 at 15:57
  • @IKavanagh Ah yes, you are right. Added back the `L` part. – Divakar Nov 11 '15 at 16:30