Chapter 2

This is a GNU Octave version of the Fundamentals of Numerical Computation Jupyter Notebooks found at https://tobydriscoll.net/project/fnc/.

Example 2.1.1 Vandermonde Interpolation

We create two vectors for data about the population of China. The first has the years of census data, the other has the numbers of millions of people.

It's convenient to measure time in years since 1980.

Now we have four data points , so and we seek an interpolating cubic polynomial. We construct the associated Vandermonde matrix:

To solve for the vector of polynomial coefficients, we use a backslash:

The algorithms used by the backslash operator are the main topic of this chapter. For now, observe that the coefficients of the cubic polynomial vary over several orders of magnitude, which is typical in this context. By our definitions, these coefficients are given in ascending order of power in . MATLAB always expects the decreasing-degree order, so we convert ours to this convention here.

We can use the resulting polynomial to estimate the population of China in 2005:

The official figure is 1297.8, so our result is not bad.

We can visualize the interpolation process. First, we plot the data as points. We'll shift the variable back to actual years.

We want to superimpose a plot of the polynomial. In order to add to a plot, we must use the hold command:

To plot the interpolating polynomial, we create a vector with many points in the time interval using linspace.

Let's clear the figure (clf) and redo it, this time continuing the curve outside of the original date range. We'll also annotate the graph (using title, xlabel, ylabel and legend) to make its purpose clear.

While the interpolation is plausible, the extrapolation to the future is highly questionable! As a rule, extrapolation more than a short distance beyond the original interval is not reliable.

Example 2.2.1 MATLAB Demo

Square brackets are used to enclose elements of a matrix or vector. Use spaces or commas for horizontal concatenation, and semicolons or new lines to indicate vertical concatenation.

A vector is considered to be a matrix with one singleton dimension.

Concatenated elements within brackets may be matrices for a block representation, as long as all the block sizes are compatible.

The dot-quote .' transposes a matrix. A single quote ' on its own performs the hermitian (transpose and complex conjugation). For a real matrix, the two operations are the same.

There are some convenient shorthand ways of building vectors and matrices other than entering all of their entries directly or in a loop. To get a row vector with evenly spaced entries between two endpoints, you have two options.

Accessing an element is done by giving one (for a vector) or two index values in parentheses. The keyword end as an index refers to the last position in the corresponding dimension.

The indices can be vectors, in which case a block of the matrix is accessed.

If a dimension has only the index : (a colon), then it refers to all the entries in that dimension of the matrix.

The matrix and vector senses of addition, subtraction, scalar multiplication, multiplication, and power are all handled by the usual symbols. If matrix sizes are such that the operation is not defined, an error message will result.

A*B causes an error. (We trap it here using a special syntax.)

A square matrix raised to an integer power is the same as repeated matrix multiplication.

In many cases, one instead wants to treat a matrix or vector as a mere array and simply apply a single operation to each element of it. For multiplication, division, and power, the corresponding operators start with a dot.

A*C would be an error.

The two operands of a dot operator have to have the same size---unless one is a scalar, in which case it is expanded or ``broadcast'' to be the same size as the other operand.

Most of the mathematical functions, such as cos, sin, log, exp and sqrt, also operate elementwise on a matrix.

Example 2.3.2 Backlash Solve

For a square matrix , the command A\b is mathematically equivalent to .

One way to check the answer is to compute a quantity known as the residual. It is (hopefully) close to machine precision, scaled by the size of the entries of the data.

If the matrix is singular, a warning is produced, but you get an answer anyway.

When you get a warning, it's important to check the result rather than blindly accepting it as correct.

Example 2.3.3 Triangular Solve

It's easy to get just the lower triangular part of any matrix using the tril command.

We'll set up and solve a linear system with this matrix.

It's not clear what the error in this answer is. However, the residual, while not zero, is virtually in size.

Next we'll engineer a problem to which we know the exact answer. You should be able to convince yourself that for any and , (missing equation)

Note: Here we need to insert some functions.

Everything seems OK here. But another example, with a different value for , is more troubling.

It's not so good to get four digits of accuracy after starting with sixteen! But the source of the error is not hard to track down. Solving for performs in the first row. Since is so much smaller than , this a recipe for losing digits to subtractive cancellation.

Example 2.4.1 Gaussian Elimination

We create a 4-by-4 linear system with the matrix

and with the right-hand side

We define an augmented matrix by tacking on the end as a new column.

The goal is to introduce zeros into the lower triangle of this matrix. By using only elementary row operations, we ensure that the matrix always represents a linear system that is equivalent to the original. We proceed from left to right and top to bottom. The first step is to put a zero in the (2,1) location using a multiple of row 1:

We repeat the process for the (3,1) and (4,1) entries.

The first column has the zero structure we want. To avoid interfering with that, we no longer add multiples of row 1 to anything. Instead, to handle column 2, we use multiples of row 2. We'll also exploit the highly repetitive nature of the operations to write them as a loop.

We finish out the triangularization with a zero in the (4,3) place. It's a little silly to use a loop for just one iteration, but the point is to establish a pattern.

Recall that is an augmented matrix: it represents the system , where

The solutions to this system are identical to those of the original system, but this one can be solved by backward substitution.

Example 2.4.2

We revisit the previous example using algebra to express the row operations on A.

We use the identity and its columns heavily.

The first step is to put a zero in the (2,1) location using a multiple of row 1:

We repeat the process for the (3,1) and (4,1) entries.

And so on, following the pattern as before.

Example 2.4.3 LU Factorization

We add a new function.

Because MATLAB doesn't show all the digits by default, it's best to compare two quantities by taking their difference.

(Usually we can expect zero only up to machine precision. However, all the exact numbers in this example are also floating-point numbers.)

To solve a linear system, we no longer need the matrix A.

Example 2.5.3 Flops 1

Here is a simple algorithm implementing the multiplication of an matrix and an vector. We use only scalar operations here for maximum transparency.

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

Since the matrix has 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 .

Let's run an experiment with the built-in matrix-vector multiplication. We use tic and toc to start and end timing of the computation.

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

Example 2.5.4 Flops 2

Let's repeat the experiment of the previous figure for more, and larger, values of .

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

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

The full story of the execution times is complicated, but the asymptotic approach to is unmistakable.

Example 2.5.5 Flops 3

We'll test the conclusion of flops experimentally, using the built-in lu function instead of the purely instructive lufact.

We plot the performance on a log-log graph and compare it to . The result could vary significantly from machine to machine.

Example 2.6.1 Pivot Fail

Here is the previously solved system.

It has a perfectly good solution, obtainable through LU factorization.

first we need a function:

If we swap the second and fourth equations, nothing essential is changed, and MATLAB still finds the solution.

However, LU factorization fails.

Example 2.6.2 Pivot Fix

Here is the system that ``broke" LU factorization for us.

When we use the built-in lu function with three outputs, we get the elements of the PLU factorization.

We can solve this as before by incorporating the permutation.

However, if we use just two outputs with lu, we get as the first result.

MATLAB has engineered the backslash so that systems with triangular or permuted triangular structure are solved with the appropriate style of triangular substitution.

The pivoted factorization and triangular substitutions are done silently and automatically when backslash is called on the original matrix.

Example 2.7.1 Vector Norms

In MATLAB one uses the norm command.

Example 2.7.2 Matrix Norms

The default norm returned by the norm command is the 2-norm.

You can get the 1-norm as well.

The 1-norm is equivalent to

Similarly, we can get the -norm and check our formula for it.

Here we illustrate the geometric interpretation of the 2-norm. First, we will sample a lot of vectors on the unit circle in .

We can apply to every column of simply by using

We superimpose the image of the unit circle with the circle whose radius is , and display multiple plots with the subplot command.

Example 2.8.1 Conditioning

MATLAB has a function cond to compute . The family of Hilbert matrices is famously badly conditioned. Here is the case.

Next we engineer a linear system problem to which we know the exact answer.

Now we perturb the data randomly but with norm .

We solve the perturbed problem using built-in pivoted LU and see how the solution was changed.

Here is the relative error in the solution.

And here are upper bounds predicted using the condition number of the original matrix.

Even if we don't make any manual perturbations to the data, machine epsilon does when we solve the linear system numerically.

Because , it's possible to lose 8 digits of accuracy in the process of passing from and to . That's independent of the algorithm; it's inevitable once the data are expressed in double precision.

Now we choose an even more poorly conditioned matrix from this family.

Before we compute the solution, note that exceeds 1/eps. In principle we might end up with an answer that is completely wrong.

MATLAB will notice the large condition number and warn us not to expect much from the result.

In fact the error does exceed 100%.

Example 2.9.1 Banded

Here is a matrix with both lower and upper bandwidth equal to one. Such a matrix is called tridiagonal.

We can extract the elements on any diagonal using the diag command. The "main'' or central diagonal is numbered zero, above and to the right of that is positive, and below and to the left is negative.

We can also put whatever numbers we like onto any diagonal with the diag command.

Here is what happens when we factor this matrix without pivoting.

Observe that the factors preserve the lower and upper bandwidth of the original matrix. However, if we introduce row pivoting, this structure may be destroyed.

Example 2.9.2 Banded Timing

We'll use a large banded matrix to observe the speedup possible in LU factorization.

If we factor the matrix as is, MATLAB has no idea that it could be exploiting the fact that it is tridiagonal.

But if we convert the matrix to a sparse one, the time required gets a lot smaller.

Example 2.9.3 Symmetric LU

We begin with a symmetric .

Carrying out our usual elimination in the first column leads us to

But now let's note that if we transpose this result, we have the same first column as before! So we could apply again and then transpose back.

Using transpose identities, this is just

Now you can see how we proceed down and to the right, eliminating in a column and then symmetrically in the corresponding row.

Finally, we arrive at a diagonal matrix.

Example 2.9.4 Cholesky

Here is a simple trick for turning any square matrix into a symmetric one.

Picking a symmetric matrix at random, there is little chance that it will be positive definite. Fortunately, the built-in Cholesky factorization chol always detects this property. The following would cause an error if run:

There is a different trick for making an SPD matrix from (almost) any other matrix.

A word of caution: chol does not check symmetry; in fact, it doesn't even look at the lower triangle of the input matrix.