Complexity of Gaussian elimination (LU factorization)

Question

Can someone tell me how to calculate the order of a) $LU$ decomposition as well as b) Gaussian elimination of a square matrix $A$? I am at a loss ...

Given:: $A$ is a $n\times n$ matrix and ultimate aim is to find a solution set of $Ax=b$. I couldn't solve the first part but for the 2nd part, I think the number of basic operations will be (as per my calculations) $$\frac{n}{2}(n+3)(2n-1)$$

Am I right? It shall even be helpful if you could cite a link or reference to some text which contain the calculations of the order of the processes.

Answer 1

You are slightly off as the dominant term should be $\frac{2}{3}n^3$. Most text in numerical linear algebra only provide the dominant error term. I am aware of a single text, specifically:

Fraleigh and Beauregard: Linear Algebra 2nd Edition, Addison-Wesley, 1990.

which does the exact flop count. I imagine that the 3rd edition from 1995 is no different.

Counting flops is difficult to do correctly. Here are some tips.

1. Write the algorithm out by hand.
2. Count the number of flops in each arithmetic statement.
3. Count the length of each loop.
4. Determine the number of times each arithmetic statement is executed.

Personally, I can never get the flop count right, unless I execute those steps.

To settle the matter decisively you can write an implementation and force the computer to count flops, comparing the actual count with your current conjecture. I shall demonstrate using MATLAB.

The following routine does an LU factorization without pivoting overwriting the matrix A with the LU factorization.

function [Y, count, diff]=mylu(A)

n=size(A,2); count=0;
for j=1:n-1
    for i=j+1:n
        A(i,j)=A(i,j)/A(j,j);
        count=count+1;
    end
    for k=j+1:n
        for i=j+1:n
            A(i,k)=A(i,k)-A(i,j)*A(j,k);
            count=count+2;
        end
    end
end
Y=A;

% Total multipliers computed (n-1)*n/2;
% Total updates made 1 + 2^2 + .... + (n-1)^2 = (n-1)*n*(2n-1)/6
% Total flops (n-1)*n*(2n-1)/3 + (n-1)*n/2;

diff=(n-1)*n*(2*n-1)/3+(n-1)*n/2-count;

The following routines solves the lower unit triangular linear system $Ly=b$ where the strictly lower triangular part of $L$ is embedded in the strict lower triangular part of $A$.

function [y, count, diff]=myforward(A,b)

n=size(A,2); y=b; count=0;
% Loop length n
for j=1:n
    % Loop length n-j
    for i=j+1:n
        % Two flops here
        y(i)=y(i)-A(i,j)*y(j);
        % This statement is executed 1 + 2 + ... (n-1) = (n-1)n/2 times
        count=count+2;
    end
end
% Total flop count = n(n-1)
diff=(n-1)*n-count;

The following routine solves the upper triangular linear system $Ux = y$ where $U$ is the upper triangular part of $A$.

function [x, count, diff]=mybackward(A,y)

n=size(A,2); count=0;

x=y;
for j=n:-1:1
    % One flop here
    x(j)=x(j)/A(j,j);
    % Update counter
    count=count+1;
    % Loop length j-1;
    for i=1:j-1
        % Two flops here. This is done (n-1)n/2 times
        x(i)=x(i)-A(i,j)*x(j);
        % Update the counter
        count=count+2;
    end
end

% Flopcount (n-1)n + n = n^2;
diff=n^2-count;