A problem \(f(x)\) can typically only be approximated in finite precision (e.g., floating-point representation). A complete set of instructions for transforming data into a result is called an algorithm. In most cases, we can represent the output from the algorithm by another function \(\tilde{f}(x)\), and for the next few lectures, we will consider the algorithm to be executed in finite precision.
Suppose our problem is to evaluate the value of a polynomial given a real input \(x\), \[f(x) = 5x^3 + 4x^2 + 3x + 2.\]
Here are two possible algorithms for evaluating this polynomial:
This savings may seem small, but saving 25% of the total operations can be huge when the numbe of operations is in the millions or billions!
The technique we just saw in the previous example is known as Horner’s rule or algorithm We have that \(p(x) = c_1 + c_2x + c_3x^2 + \cdots + c_nx^{n-1} = c_1 + x(c_2 + \cdots + x(c_{n-2} + x(c_{n-1} + c_nx))).\)
"""
horner(c,x)
Evaluate a polynomial whose coefficients are given in ascending (according to
associated monomial degree) order in 'c', at the point 'x' using Horner's
algorithm.
"""
function horner(c,x)
n = length(c)
y = c[n]
for k in n-1:-1:1
y = x*y + c[k]
end
return y
endhorner (generic function with 1 method)In the code above, the length function returns the number of elements in vector c. We use c[n] to access the \(n\)th element of vector c. The polynomial value is evaluated recursively in the for loop. Note the format for the range for k – it ranges from \(n-1\) to 1 with steps of size \(-1\).
Let’s use this function to evaluate the value a given polynomial!
c = [-1,3,-3,1]
horner(c,1.6)0.21600000000000041If we solve a problem \(f\) using a computer algorithm and we see a large error in the result, we might suspect that the problem has poor conditioning. However, it could also be that the algorithm introduced additional error. When error in the output of an algorithm exceeds what the problem conditioning explains, we say that the algorithm is unstable.
We’re returning to the problem of calculating roots of a quadratic polynomial; that is finding values \(t\) such that \(at^2 + bt + c = 0\). In Lecture 1, we showed that this problem is ill-conditioned if and only if the roots are close together relative to their size.
Thus, find the roots of the polynomial \(p(x) = (x- 10^6)(x - 10^{-6}) = x^2 - (10^6 + 10^{-6})x + 1\) is a well-conditioned problem. As we saw previously, the quadratic formula is an algorithm for this problem.
a = 1; b = -(1e6+1e-6); c = 1;
@show x₁ = (-b + sqrt(b^2 - 4a*c)) / 2a;
@show x₂ = (-b - sqrt(b^2 - 4a*c)) / 2a;x₁ = (-b + sqrt(b ^ 2 - (4a) * c)) / (2a) = 1.0e6
x₂ = (-b - sqrt(b ^ 2 - (4a) * c)) / (2a) = 1.00000761449337e-6The larger root has no error, but now we measure the relative error in the smaller root!
error = abs(1e-6 - x₂) / 1e-6
@show accurate_digits = -log10(error);accurate_digits = -(log10(error)) = 5.118358987126217@show u₁ = b^2;
@show u₂ = u₁ - 4;
@show u₃ = sqrt(u₂);
@show u₄ = -u₃ - b;
@show u₅ = u₄/2;u₁ = b ^ 2 = 1.000000000002e12
u₂ = u₁ - 4 = 9.99999999998e11
u₃ = sqrt(u₂) = 999999.999999
u₄ = -u₃ - b = 2.00001522898674e-6
u₅ = u₄ / 2 = 1.00000761449337e-6Since \(a = c = 1\), we’ll assume these are represented exactly, and consider only the condition number with respect to \(b\). We can see where the instability in the smaller root calculation came from by considering the condition number of each of the sub-problems encountered.
| Problem \(f\) | \(\kappa_f\) |
|---|---|
| \(u_1 = f_1(b) = b^2\) | \(\kappa_{f_1}(b) = 2\) |
| \(u_2 = f_2(u_1) = u_1 - 4\) | \(\kappa_{f_2}(u_1) \approx 1\) |
| \(u_3 = f_3(u_2) = \sqrt{u_2}\) | \(\kappa_{f_3}(u_2) = 1/2\) |
| \(u_4 = f_4(u_3) = -(u_3 + b)\) | \(\kappa_{f_4}(u_3) \approx 5 \times 10^{11}\) |
| \(u_5 = f_5(u_4) = u_4/2\) | \(\kappa_{f_5}(u_4) = 1\) |
We expect to lose 11 digits of accuracy in the fourth step of this algorithm. Here the issue is the subtractive cancellation between \(\sqrt{b^2 - 4ac}\) and \(b\)!
We can compute this root with no error using a different algorithm!
@show x₂ = c / (a * x₁);
abs(x₂ - 1e-6) / 1e-6x₂ = c / (a * x₁) = 1.0e-60.0These two algorithms are equivalent when using real numbers and exact arithmetic, but the outputs they calculate in practice are perturbed by finite precision representation in each step and depend upon the specific order of operation.
This may seem scary and complicated, but most simple operations are well-conditioned most of the time! Exceptions are usually due to \(|f(x)|\) being much smaller than \(|x|\), and the most common culprit (by far!) is subtractive cancellation.
If a problem \(f(x)\) has poor conditioning, event just the act of rounding the data to floating-point representation may introduce large errors in the result. It’s not reasonable to expect that algorithms \(\tilde{f}\) will have small error in the sense that \(\tilde{f}(x) \approx f(x)\).
Numerical analysts instead prefer to characterize the error in a different way – instead of asking “Did you get nearly the right answer?”, we ask “Did you answer nearly the right question?”
Backward error analysis causes us to ask “What is the problem our algorithm actually solved?” and to measure the distance between the ideal data and this alternative input to \(f\).
using Pkg
Pkg.add("Polynomials")
using Polynomials Resolving package versions...
Updating `~/.julia/environments/v1.11/Project.toml`
[f27b6e38] + Polynomials v4.0.13
No Changes to `~/.julia/environments/v1.11/Manifest.toml`We need access to some code from the Julia package Polynomials. We first must use the pre-installed Julia package Pkg which we allow by using the command using Pkg. The command Pkg.add("Polynomials") installs the package Polynomials for the Julia installation to access. Note: we only have to do this once! Then, forevermore, we can use the command using Polynomials to allow access to these functions.
We can build the polynomial \(p\) from its roots, \(-2, -1, 1, 1, 3, 6\). We do this so we know the exact roots of the polynomial.
r = [-2.0, -1, 1, 1, 3, 6]
p = fromroots(r)Now, we compute the roots of \(p\) using an algorithm and see how much error has been introduced into each root.
r̃ = sort(roots(p))6-element Vector{Float64}:
-1.9999999999999993
-1.0000000000000002
0.9999999876576552
1.0000000123423434
3.0000000000000036
5.999999999999993println("Root errors:")
@. abs(r - r̃) / rRoot errors:6-element Vector{Float64}:
-3.3306690738754696e-16
-2.220446049250313e-16
1.2342344812843464e-8
1.2342343369553532e-8
1.1842378929335002e-15
1.1842378929335002e-15Next, we build a polynomial \(\tilde{p}\) from these approximate roots and check how far these output coefficients have deviated from the true coefficients.
p̃ = fromroots(r̃)c, c̃ = coeffs(p), coeffs(p̃)
println("Coefficients errors:")
@. abs(c - c̃) / cCoefficients errors:7-element Vector{Float64}:
1.7763568394002505e-15
-1.9737298215558337e-16
-1.6524249668839539e-15
6.459479416000911e-16
2.9605947323337506e-16
-5.551115123125783e-16
0.0Even though some computed roots were relatively far from the exact values, they are roots of a polynomial that is nearby to the ideal polynomial! In other words, we solved a problem nearby to the original problem, even if the results were quite far apart.
As an example, the algorithm \(f(x) = x + 1\) is stable, but not backward stable. If \(|x| < \epsilon_{\text{mach}}/2,\) then the computed result is \(\tilde{f}(x) = 1\) since there are no floating points between \(1\) and \(1 + \epsilon_{\text{mach}}\).
Hence, the only choice for a real \(\tilde{x}\) so that \(f(\tilde{x}) = \tilde{f}(x) = 1\) is \(\tilde{x} = 0\). Then \(|\tilde{x} - x|/|x| = 1\) – 100% backward error!
There are two important (and equivalent) views of matrix multiplication. Let \(\mathbf{A} \in \mathbb{R}^{m \times n}, \mathbf{B} \in \mathbb{R}^{n \times p}\) and \(\mathbf{C} = \mathbf{A}\mathbf{B} \in \mathbb{R}^{m \times p}\).
The entries of \(\mathbf{C} = \mathbf{A}\mathbf{B}\) may be calculated by the inner product, \[\mathbf{C}_{ij} = \sum_{k = 1}^n \mathbf{A}_{ik}\mathbf{B}_{kj}.\]
The alternative view of matrix multiplication is as a sum of rank-one matrices formed as outer products of corresponding columns of \(\mathbf{A}\) and rows of \(\mathbf{B}\), \[\mathbf{C} = \sum_{k = 1}^n \mathbf{A}_{:k} \mathbf{B}_{k:}.\]
Using the previous interpretations of matrix multiplications, we can better understand the special case of matrix-vector multiplication. Consider computing \(\mathbf{A}\mathbf{v}\) for matrix \(\mathbf{A} \in \mathbb{R}^{m times n}\) and vector \(\mathbf{v} \in \mathbb{R}^n\).
We have that \[\mathbf{A}\mathbf{v} = v_1 \mathbf{a}_1 + \cdots + v_n \mathbf{a}_n\] where \(\mathbf{a}_i\) is the \(i\)th columns of \(\mathbf{A}\).
We may transpose the matrix-vector product to get:
Recall the following important theorem:
Recall that this means \(\mathbf{A} = \mathbf{Q}^\top \mathbf{D} \mathbf{Q}\).
