Efficiency of matrix computations - Fundamentals of Numerical Computation

Predicting how long an algorithm will take to solve a particular problem, on a particular computer, as written in a particular way in a particular programming language, is an enormously difficult undertaking. It’s more practical to predict how the required time will scale as a function of the size of the problem. In the case of a linear system of equations, the problem size is $n$ , the number of equations and variables in the system. Because expressions of computational time are necessarily approximate, it’s customary to suppress all but the term that is dominant as $n\to\infty$ . We first need to build some terminology for these expressions.

2.5.1Asymptotic analysis¶

One immediate consequence is that $f\sim g$ implies $f=O(g)$ .^[1]

Example 2.5.2

Consider $f(n) = \sin (1/n)$ , $g(n)=1/n$ and $h(n) = 1/n^2$ . For large $n$ , Taylor’s theorem with remainder implies that

f(n) = \frac{1}{n} - \cos(1/\xi)\frac{1}{6 n^3},

(2.5.2)

where $n<\xi<\infty$ . But

\lim_{n\to \infty} \frac{f}{g} = \lim_{n\to \infty} 1-\cos(1/\xi)\frac{1}{6 n^2} = 1,

(2.5.3)

and so $f \sim g$ . On the other hand, comparing $f$ and $h$ , we find

\lim_{n\to \infty} \frac{f}{h} = \lim_{n\to \infty} n-\cos(1/\xi)\frac{1}{6 n} = \infty,

(2.5.4)

so we cannot say that $f = O(h)$ . A consideration of $h/f$ will show that $h = O(f)$ , however.

It’s conventional to use asymptotic notation that is as specific as possible. For instance, while it is true that $n^2+n=O(n^{10})$ , it’s more informative, and usually expected, to say $n^2+n=O(n^2)$ . There are additional notations that enforce this requirement strictly, but we will just stick to the informal understanding.

There is a memorable way to use asymptotic notation to simplify sums:

\begin{split} \sum_{k=1}^n k&\sim \frac{n^2}{2} = O(n^2), \text{ as $n\to\infty$}, \\ \sum_{k=1}^n k^2 &\sim \frac{n^3}{3} = O(n^3), \text{ as $n\to\infty$}, \\ &\vdots \\ \sum_{k=1}^n k^p &\sim \frac{n^{p+1}}{p+1} = O(n^{p+1}), \text{ as $n\to\infty$}. \end{split}

(2.5.5)

These formulas greatly resemble the definite integral of $x^p$ .

2.5.2Flop counting¶

Traditionally, in numerical linear algebra, we count arithmetic operations to measure an algorithm’s running time.

Example 2.5.4 (Floating-point operations in matrix-vector multiplication)

Julia

MATLAB

Python

Example 2.5.4

Here is a straightforward implementation of matrix-vector multiplication.

n = 6
A = randn(n, n)
x = rand(n)
y = zeros(n)
for i in 1:n
    for j in 1:n
        y[i] += A[i, j] * x[j]    # 1 multiply, 1 add
    end
end

Each of the loops implies a summation of flops. The total flop count for this algorithm is

\sum_{i=1}^n \sum_{j=1}^n 2 = \sum_{i=1}^n 2n = 2n^2.

(1)

Since the matrix $\mathbf{A}$ has $n^2$ elements, all of which have to be involved in the product, it seems unlikely that we could get a flop count that is smaller than $O(n^2)$ in general.

Let’s run an experiment with the built-in matrix-vector multiplication. Note that Julia is unusual in that loops have a variable scope separate from its enclosing code. Thus, for n in n below means that inside the loop, the name n will take on each one of the values that were previously assigned to the vector n.

n = 1000:1000:5000
t = []
for n in n
    A = randn(n, n)  
    x = randn(n)
    time = @elapsed for j in 1:80; A * x; end
    push!(t, time)
end

The reason for doing multiple repetitions at each value of $n$ in the loop above is to avoid having times so short that the resolution of the timer is significant.

@pt :header = ["size", "time (sec.)"] [n t]

Looking at the timings just for $n=2000$ and $n=4000$ , they have ratio

@show t[n.==4000] ./ t[n.==2000];

t[n .== 4000] ./ t[n .== 2000] = [2.8181890443804565]

If the run time is dominated by flops, then we expect this ratio to be

\frac{2(4000)^2}{2(2000)^2}=4.

(2)

Example 2.5.4

Here is a straightforward implementation of matrix-vector multiplication.

n = 6;
A = magic(n);
x = ones(n,1);
y = zeros(n,1);
for i = 1:n
    for j = 1:n
        y(i) = y(i) + A(i,j)*x(j);   % 2 flops
    end
end

Each of the loops implies a summation of flops. The total flop count for this algorithm is

\sum_{i=1}^n \sum_{j=1}^n 2 = \sum_{i=1}^n 2n = 2n^2.

(1)

Since the matrix $\mathbf{A}$ has $n^2$ elements, all of which have to be involved in the product, it seems unlikely that we could get a flop count that is smaller than $O(n^2)$ in general.

Let’s run an experiment with the built-in matrix-vector multiplication, using tic and toc to time the operation.

n = (500:500:5000)';
t = zeros(size(n));
for k = 1:length(n)
    A = randn(n(k), n(k));  x = randn(n(k), 1);
    tic    % start a timer
    for j = 1:200      % repeat 100 times
        A*x;
    end
    time = toc;           % read the timer
    t(k) = time / 200;   % seconds per instance
end

The reason for doing multiple repetitions at each value of $n$ in the loop above is to avoid having times so short that the resolution of the timer is significant.

table(n, t, 'variablenames', ["size", "time"])

Looking at the timings just for $n=2500$ and $n=5000$ , they have ratio

t(n==5000) / t(n==2500)

If the run time is dominated by flops, then we expect this ratio to be

\frac{2(5000)^2}{2(2500)^2}=4.

(2)

Example 2.5.4

Here is a straightforward implementation of matrix-vector multiplication.

n = 6
A = random.rand(n, n)
x = ones(n)
y = zeros(n)
for i in range(n):
    for j in range(n):
        y[i] += A[i, j] * x[j]   # 2 flops

Each of the loops implies a summation of flops. The total flop count for this algorithm is

\sum_{i=1}^n \sum_{j=1}^n 2 = \sum_{i=1}^n 2n = 2n^2.

(1)

Since the matrix $\mathbf{A}$ has $n^2$ elements, all of which have to be involved in the product, it seems unlikely that we could get a flop count that is smaller than $O(n^2)$ in general.

Let’s run an experiment with the built-in matrix-vector multiplication. We assume that flops dominate the computation time and thus measure elapsed time.

N = 400 * arange(1, 11)
t = []
print("  n           t")
for i, n in enumerate(N):
    A = random.randn(n, n)  
    x = random.randn(n)
    start = timer()
    for j in range(50): A @ x
    t.append(timer() - start)
    print(f"{n:5}   {t[-1]:10.3e}")

  n           t
  400    1.054e-03
  800    2.666e-02
 1200    3.578e-02

 1600    7.288e-02

 2000    8.508e-02

 2400    1.993e-01

 2800    6.765e-02

 3200    2.149e-01

 3600    1.095e-01

 4000    1.280e-01

The reason for doing multiple repetitions at each value of $n$ above is to avoid having times so short that the resolution of the timer is a factor.

Looking at the timings just for $n=2000$ and $n=4000$ , they have ratio:

print(t[9] / t[4])

1.50431067089939

If the run time is dominated by flops, then we expect this ratio to be

\frac{2(4000)^2}{2(2000)^2}=4.

(2)

Suppose that the running time $t$ of an algorithm obeys a function that is $O(n^p)$ . For sufficiently large $n$ , $t\approx Cn^p$ for a constant $C$ should be a good approximation. Hence

t \approx Cn^p \qquad \Longleftrightarrow \qquad \log t \approx p(\log n) + \log C.

(2.5.6)

So, we expect that a graph of $\log t$ as a function of $\log n$ will be a straight line of slope $p$ .

Example 2.5.5 (Asymptotics in log-log plots)

Julia

MATLAB

Python

Example 2.5.5

Let’s repeat the experiment of the previous figure for more, and larger, values of $n$ .

randn(5,5)*randn(5);  # throwaway to force compilation

n = 400:200:6000
t = []
for n in n
    A = randn(n, n)  
    x = randn(n)
    time = @elapsed for j in 1:50; A * x; end
    push!(t, time)
end

Plotting the time as a function of $n$ on log-log scales is equivalent to plotting the logs of the variables.

scatter(n, t, label="data", legend=false,
    xaxis=(:log10, L"n"), yaxis=(:log10, "elapsed time (sec)"),
    title="Timing of matrix-vector multiplications")

You can see that while the full story is complicated, the graph is trending to a straight line of positive slope. For comparison, we can plot a line that represents $O(n^2)$ growth exactly. (All such lines have slope equal to 2.)

plot!(n, t[end] * (n/n[end]).^2, l=:dash,
    label=L"O(n^2)", legend=:topleft)

Example 2.5.5

Let’s repeat the previous experiment for more, and larger, values of $n$ .

n = (400:400:6000)';
t = zeros(size(n));
for k = 1:length(n)
    A = randn(n(k), n(k));  x = randn(n(k), 1);
    tic    % start a timer
    for j = 1:100      % repeat ten times
        A*x;
    end
    time = toc;          % read the timer
    t(k) = time / 100;   % seconds per instance
end

Plotting the time as a function of $n$ on log-log scales is equivalent to plotting the logs of the variables.

clf    % clear any existing figure
loglog(n, t, 'o-')
xlabel('size of matrix')
ylabel('time (sec)')
title('Timing of matrix-vector multiplications')

hold on
loglog(n, 0.5 * t(end) * (n / n(end)).^2, 'k--')
axis tight
legend('data', 'O(n^2)', 'location', 'southeast')

Example 2.5.5

Let’s repeat the experiment of the previous example for more, and larger, values of $n$ .

N = arange(400, 6200, 200)
t = zeros(len(N))
for i, n in enumerate(N):
    A = random.randn(n,n)  
    x = random.randn(n)
    start = timer()
    for j in range(20): A@x
    t[i] = timer() - start

Plotting the time as a function of $n$ on log-log scales is equivalent to plotting the logs of the variables, but is formatted more neatly.

fig, ax = subplots()
ax.loglog(N, t, "-o", label="observed")
ylabel("elapsed time (sec)");
xlabel("$n$");
title("Timing of matrix-vector multiplications");

ax.loglog(N, t[-1] * (N/N[-1])**2, "--", label="$O(n^2)$")
ax.legend();  fig

2.5.3Solution of linear systems¶

Recall the steps of Algorithm 2.4.2 for the system $\mathbf{A}\mathbf{x}=\mathbf{b}$ :

Factor $\mathbf{L}\mathbf{U}=\mathbf{A}$ using Gaussian elimination.
Solve $\mathbf{L}\mathbf{z}=\mathbf{b}$ for $\mathbf{z}$ using forward substitution.
Solve $\mathbf{U}\mathbf{x}=\mathbf{z}$ for $\mathbf{x}$ using backward substitution.

The second and third steps are solved by Function 2.3.1 and Function 2.3.2.

Julia

MATLAB

Python

Take forwardsub, for instance. It has a single flop in line 10. Line 12 computes

sum(L[i, j] * x[j] for j in 1:i-1)

This line requires $i-1$ multiplications and $(i-2)$ additions, for a total of $2i-3$ flops. Line 13 adds two more flops, to make that $2n-1$ . These lines are performed within a loop as $i$ ranges from 2 to $n$ , so the total count is

1 + \sum_{i=2}^n (2i-1) = \sum_{i=1}^n (2i-1) = -n + 2 \sum_{i=1}^n i.

Take forwardsub, for instance. All of the flops are in line 12:

(b(i) - L(i, 1:i-1) * x(1:i-1)) / L(i, i);

The first operation to be performed in this line is the inner product L(i, 1:i-1) * x(1:i-1), which requires $i-1$ multiplications and $(i-2)$ additions, for a total of $2i-3$ flops. Then there is one subtraction and one division, which adds two more flops. These lines are performed within a loop as $i$ ranges from 1 to $n$ , so the total count is

\sum_{i=1}^n (2i-1) = -n + 2 \sum_{i=1}^n i.

Take forwardsub, for instance. Line 11 computes

s = L[i, :i] @ x[:i]

This is the inner product between the first $i$ elements of the $i$ th row of L and the first $i$ elements of x, requiring $i$ multiplications and $(i-1)$ additions. Line 212 adds two more flops, for a total of $2i+1$ in each pass through the loop. In this loop, $i$ ranges from 0 to $n-1,$ so the total count is

\sum_{i=0}^{n-1} (2i+1) = \sum_{i=1}^n (2i-1) = -n + 2 \sum_{i=1}^n i.

It is not hard to find an exact formula for the sum above, but we use the asymptotic expression (2.5.5) to simplify it to $\sim n^2$ . After all, since flop counting is only an approximation of true running time, why bother with the more complicated exact expression? An analysis of backward substitution yields the same result.

Before counting flops for the LU factorization, we have to admit that Function 2.4.1 is not written as economically as it could be. Recall from our motivating example in Example 2.4.3 that we zero out the first row and column of $\mathbf{A}$ with the first outer product, the second row and column with the second outer product, and so on. There is no good reason to do multiplications and additions with values known to be zero.

Julia

MATLAB

Python

Suppose we replace lines 14–18 of lufact in Algorithm 2.4.1 with

14
15
16
17
18
for k in 1:n-1
    U[k, k:n] = Aₖ[k, k:n]
    L[k:n, k] = Aₖ[k:n, k] / U[k, k]
    Aₖ[k:n, k:n] -= L[k:n, k] * U[k, k:n]'
end

We will use the following handy fact.

Line 16 above divides each element of the vector Aₖ[k:n, k] by a scalar. Hence, the number of flops equals the length of the vector, which is $n-k+1$ .

Line 17 has an outer product followed by a matrix subtraction. The definition (A.5) of the outer product makes it clear that that computation takes one flop (multiplication) per element of the result, which here results in $(n-k+1)^2$ flops. The number of subtractions is identical.

Altogether, the factorization takes

\sum_{k=1}^{n-1} n-k + 1 + 2(n-k+1)^2

flops.

Suppose we replace lines 14–18 of lufact in Algorithm 2.4.1 with

14
15
16
17
18
for k = 1:n-1
    U(k, k:n) = A_k(k, k:n);
    L(k:n, k) = A_k(k:n, k) / U(k, k);
    A_k(k:n, k:n) = A_k(k:n, k:n) - L(k:n, k) * U(k, k:n);
end

We will use the following handy fact.

Line 16 above divides each element of the vector A_k(k:n, k) by a scalar. Hence, the number of flops equals the length of the vector, which is $n-k+1$ .

Altogether, the factorization takes

\sum_{k=1}^{n-1} n-k + 1 + 2(n-k+1)^2

flops.

Suppose we replace lines 14–17 of lufact in Algorithm 2.4.1 with

14
15
16
17
for k in range(n-1):
    U[k, k:n] = A_k[k, k:n]
    L[k:n, k] = A_k[k:n, k] / U[k, k]
    A_k[k:n, k:n] -= np.outer(L[k:n, k], U[k, k:n])

We will use the following handy fact.

Line 16 above divides each element of the vector A_k[k:n, k] / U[k, k] by a scalar. Hence, the number of flops equals the length of the vector, which is $n-k$ .

Line 17 has an outer product followed by a matrix subtraction. The definition (A.5) of the outer product makes it clear that that computation takes one flop (multiplication) per element of the result, which here results in $(n-k)^2$ flops. The number of subtractions is identical.

Altogether, the factorization takes

\sum_{k=0}^{n-2} n-k + 2(n-k)^2 = \sum_{k=1}^{n-1} n-k + 1 + 2(n-k+1)^2

flops.

There are different ways to simplify the total count,

\sum_{k=1}^{n-1} n-k + 1 + 2(n-k+1)^2.

(2.5.7)

We will make a change of summation index using $j=n-k$ . The endpoints of the sum are $j=n-1$ when $k=1$ and $j=1$ when $k=n-1$ . Since the order of terms in a sum doesn’t matter, we get

\begin{align*} \sum_{j=1}^{n-1} 1+j+2(j+1)^2 &= \sum_{j=1}^{n-1} 3 + 5j + 2j^2 \\ & \sim 3(n-1) + \frac{5}{2}(n-1)^2 + \frac{2}{3}(n-1)^3 \\ & \sim \frac{2}{3}n^3. \end{align*}

(2.5.8)

We have proved the following.

Example 2.5.6 (Floating-point operations in LU factorization)

Julia

MATLAB

Python

Example 2.5.6

We’ll test the conclusion of $O(n^3)$ flops experimentally, using the built-in lu function instead of the purely instructive lufact.

lu(randn(3, 3));   # throwaway to force compilation

n = 400:400:4000
t = []
for n in n
    A = randn(n, n)  
    time = @elapsed for j in 1:12; lu(A); end
    push!(t, time)
end

We plot the timings on a log-log graph and compare it to $O(n^3)$ . The result could vary significantly from machine to machine, but in theory the data should start to parallel the line as $n\to\infty$ .

scatter(n, t, label="data", legend=:topleft,
    xaxis=(:log10, L"n"), yaxis=(:log10, "elapsed time"))
plot!(n, t[end ]* (n/n[end]).^3, l=:dash, label=L"O(n^3)")

Example 2.5.6

We’ll test the conclusion of $O(n^3)$ flops experimentally, using the built-in lu function instead of the purely instructive lufact.

n = (200:100:2400)';
t = zeros(size(n));
for k = 1:length(n)
    A = randn(n(k), n(k));  
    tic    % start a timer
    for j = 1:6,  [L, U] = lu(A);  end
    time = toc;
    t(k) = time / 6;  
end

clf
loglog(n, t, 'o-')
hold on
loglog(n, 0.5 * t(end) * (n/n(end)).^3, 'k--')
axis tight
xlabel('size of matrix'), ylabel('time (sec)')
title('Timing of LU factorization')
legend('lu','O(n^3)','location','southeast');

Example 2.5.6

We’ll test the conclusion of $O(n^3)$ flops experimentally using the lu function imported from scipi.linalg.

from scipy.linalg import lu
N = arange(200, 2600, 200)
t = zeros(len(N))
for i, n in enumerate(N):
    A = random.randn(n,n)  
    start = timer()
    for j in range(5): lu(A)
    t[i] = timer() - start

loglog(N, t, "-o", label="obseved")
loglog(N, t[-1] * (N / N[-1])**3, "--", label="$O(n^3)$")
legend();
xlabel("$n$");
ylabel("elapsed time (sec)");
title("Timing of LU factorizations");

In practice, flops are not the only aspect of an implementation that occupies significant time. Our position is that counting flops as a measure of performance is a useful oversimplification. We will assume that LU factorization (and as a result, the solution of a linear system of $n$ equations) requires a real-world time that is roughly $O(n^3)$ . This growth rate is a great deal more tolerable than, say, $O(2^n)$ , but it does mean that for (at this writing) $n$ greater than 10,000 or so, something other than general LU factorization will have to be used.

2.5.4Exercises¶

Exercise 2.5.5

✍ This problem is about evaluation of a polynomial $c_1 + c_2 x + \cdots + c_{n}x^{n-1}$ .

(a) Here is a little code to do the evaluation.

Julia

MATLAB

Python

y = c[1]
xpow = 1
for i in 2:n
    xpow *= x
    y += c[i]*xpow
end

Assuming that x is a scalar, how many flops does this function take, as a function of $n$ ?

(b) Compare the count from (a) to the flop count for Horner’s algorithm, Function 1.3.1.

Footnotes¶

More precisely, $O(g)$ and $\sim g$ are sets of functions, and $\sim g$ is a subset of $O(g)$ . That we write $f=O(g)$ rather than $f\in O(g)$ is a quirk of convention.
↩