Markdown transcription of my handwritten CS 370 notes. Spelling and wording are preserved where legible; diagrams are described in brackets, and editorial clarifications are also in brackets.

Floating Point

Normalized form of a real number $F=\lbrace\beta,t,L,U\rbrace$:

where $d_i$ are digits in base $\beta$ ($0\leq d_i<\beta$), $d_1\neq0$, and $p$ is an integer where $L\leq p\leq U$.

• IEEE single precision: $\lbrace2,24,-126,127\rbrace$
• IEEE double precision: $\lbrace2,53,-1022,1023\rbrace$
• If $p>U$: overflow.
• If $p<L$: underflow.

Rounding Modes for FP Conversion

Rounding modes for FP conversion:

1. Round to nearest number available in $F$ ($0.5$ rounds up).
2. Truncation: discard all digits after the $t$th.

$x_{\text{exact}}$ means the true solution. $x_{\text{approx}}$ means the approximate, representable solution.

A result is correct to approximately $s$ digits if

or equivalently

Therefore,

so about two digits are correct.

For an FP system $F$, $E_{\text{rel}}$ is bounded such that for any $x\in\mathbb R$ and its FP representation $\operatorname{fl}(x)$,

where $E$ is machine epsilon and is defined as the smallest number such that

Therefore,

For an FP system $F=\lbrace\beta,t,L,U\rbrace$:

1. Under rounding, $E=\frac12\beta^{1-t}$.
2. Under truncation, $E=\beta^{1-t}$.

IEEE requires that

Here $w+z$ is the exact real result represented as FP. But in sequences of operations, error compounds:

Catastrophic cancellation error: subtracting input numbers with error of about the same magnitude.

Benign cancellation: subtracting input numbers with no error of about the same magnitude.

Rules of thumb:

1. Sum numbers of the same sign and magnitude.
2. If input numbers contain error, reformulate to avoid cancellation.

Consider two formulas for the Taylor series for $e^x$, where $x<0$:

The second is much better because it avoids cancellation.

Compare

as $x\to0$.

but $f(x)$ approximates this much worse because $\cos^2x\approx1$.

Other types of error:

• round-off error (FP representation);
• truncation error (i.e. after $n$ terms of Taylor);
• uncertainty/error in input;
• error/approximation in the model of the situation.

Conditioning and Stability

For a problem $P$ with input $I$ and output $O$, if small $\Delta I$ results in small $\Delta O$, then $P$ is well-conditioned. Otherwise, $P$ is ill-conditioned.

If an algorithm is unstable, any error in the input is magnified.

Conditioning and stability are independent properties.

Intro to Interpolation

Given a set of data points from an unknown function $y=p(x)$, can we approximate $p$’s value at other points?

Univalence theorem: Given $n$ pairs $(x_i,y_i)$ where $1\leq i\leq n$, there exists a unique polynomial $p(x)$ of degree $\leq n-1$ that interpolates the data.

For $n$ points, find all coefficients $c_i$ of

Plugging in each $(x_i,y_i)$ gives one linear equation:

This gives an $n\times n$ linear system,

or $V\vec c=\vec y$, where $V$ is the Vandermonde matrix.

which means there must be a solution. This proves the univalence theorem.

Monomial form:

Monomial basis functions: $1,x,x^2,\ldots$

Lagrange Basis Function

$L_k(x)$:

There is no $x-x_k$ or $x_k-x_k$ entry.

Notice:

since numerator/denominator $=1$, and

since $x_j-x_j$ will occur in the numerator.

Lagrange Form

$p(x)$ interpolates each $(x_i,y_i)$ by definition. For example,

Runge’s Phenomenon

Fitting a high-degree polynomial (many points) often gives excessive oscillation.

We solve this with piecewise functions.

Hermite Interpolation

Fit a cubic to each interval given values and derivatives at each point.

We use the same strategies as in Vandermonde to determine the coefficients of each cubic.

Define $p_i(x)$ as the polynomial on the $i$th interval, between points $x_i$ and $x_{i+1}$:

where

$y_i’$ refers to the linear slope, not derivative $s_i$.

Knots: points where the interpolant transfers from $p_i(x)$ to $p_{i+1}(x)$.

Nodes: points where data for the interpolant is specified, i.e. slope.

Hermite interpolation removes discontinuities at each control point by definition, since it creates some polynomial $p_i(x)$ with derivatives $s_i,s_{i+1}$ at each control point. Therefore it creates a smooth curve.

Cubic Splines

Fit a cubic $S_i(x)$ on each interval; require matching first and second derivatives; require only points as input.

For each interval $[x_i,x_{i+1}]$ for $S_i(x)$ and $S_{i+1}(x)$:

$S_i(x)$ joins $S_{i+1}(x)$ at $x_{i+1}$; the interval from $x_i$ to $x_{i+1}$ is $\Delta x_i$.

Require:

Two equations are needed for boundary conditions $(x_1,x_n)$:

• Clamped: $S’(x_1)=c_1$ and $S’(x_n)=c_2$, where $c_1,c_2$ are given constants.
• Free: $S’‘(x_1)=S’‘(x_n)=0$.
• Periodic: $S_1’(x_1)=S_n’(x_n)$ and $S_1’‘(x_1)=S_n’‘(x_n)$ - wrap around.
• “Not a knot”: $S_1’’‘(x)=S_2’’‘(x)$ and $S_{n-1}’’‘(x)=S_{n-2}’’‘(x)$ - the last two segments on an end become the same polynomial.

To solve for a cubic spline over $n$ points, we need:

• two equations per interpolating interval, $(1)$ and $(2)$, for $n-1$ intervals;
• two equations per interior point, $(3)$ and $(4)$, for $n-2$ interior points;
• two boundary equations.

Therefore,

equations total, to solve in one system.

Compare with Hermite interpolation over $n$ points with $n$ slopes, where we solve $n-1$ systems of four equations, also $4(n-1)$ equations total.

Notice cubic splines also produce smooth curves.

Cost of Computing Cubic Splines

We have to solve $4n-4$ equations for $4(n-1)$ unknowns, since each of the $n-1$ intervals involves a cubic

with four unknowns.

$A$ has size $(4n-4)^2$, which costs $O(n^3)$ to compute naively and $O(n)$ to compute in a smart way.

$O(n)$ strategy: Compute cubic splines using Hermite interpolation.

1. Write each spline $S_i(x)$ in Hermite closed form:

where

1. Find $s_i$ such that $S_i’‘(x)=S_{i+1}’‘(x)$ at interior points, i.e. second derivative matches. We do not require first derivative matching because it is satisfied by the Hermite equation.

Plugging in definitions for $c_i,d_i$,

Now force $S_i’‘(x_{i+1})=S_{i+1}’‘(x_{i+1})$ for $i=1$ to $n-2$:

Equating and simplifying,

so

Reindexing with $i\to i-1$, from $i=2$ to $n-1$, gives an equation per interior node:

To find $s_i$, we need boundary conditions.

1. Clamped:

where $s_1^,s_n^$ are given.

1. Free:

Subbing the Hermite closed form into (1),

so $c_1=0$. Applying the definition for $c_i$,

therefore

Subbing the Hermite closed form into (2),

which simplifies to

Cubic Spline Equations

For each interior node $(i=2,\ldots,n-1)$,

Boundary conditions $(i=1,i=n)$:

Clamped:

Free:

In matrix form,

where $\vec r$ is the right-hand side of the equation.

For $T$, each $T_{i,j}$ describes the $i$th equation of (1), or a boundary-condition equation, and describes the coefficient of $s_j$ in that equation. For $i=2,\ldots,n-1$,

and $T_{i,k}=0$ for $k\neq i-1,i,i+1$.

For clamped conditions:

For free conditions:

Notice $T$ always has a solution and is a tridiagonal matrix, which can be solved/stored efficiently.

Using cubic splines with Hermite interpolation, we solve one equation in $s_1,\ldots,s_n$ variables for each of the $n$ nodes, giving an $n\times n$ linear system. Therefore, total size of system is $n^2$.

Recall that using cubic splines requires solving a linear system of size $(4n-4)^2$. Therefore cubic splines plus Hermite interpolation is approximately four times smaller.

Parametric Curves and ODEs

Parametric curves: $p(t)=(x(t),y(t))$. This allows for points that have the same $(x,y)$ coordinate.

A common parameterization is to use arc length as $t$. We can use interpolation on parameterized curves $x(t)$ and $y(t)$ as well.

Given $(x_i,y_i)$ data points, parameterize $t_i$ with $(x_0,y_0)$ and $(x_n,y_n)$ by:

1. Use node index as parameterization, e.g. $t_i=i$.
2. Approximate arc-length parameterization: set $t_0=0$ at the first node, then

We can apply curve concepts to surfaces by subdividing them.

An ordinary differential equation is a function $f$ relating a variable and its derivative:

Most ODEs are too complex to have simple closed forms, so we develop methods to find approximates $(t_i,y_i)$ where $y_i\approx y(t_i)$. Given $(t_i,y_i)$ pairs, use interpolation methods to derive a full solution for all times.

Initial-value problem (IVP):

Find/approximate $y(t)$.

Time-stepping: Given initial conditions, step to the next time endpoint using the dynamics function and time step $h$.

Set $n=0$, $t=t_0$, $y=y_0$:

1. Compute $y_{n+1}$.
2. Increment time, $t_{n+1}=t_n+h$.
3. Advance $n=n+1$.

Forward Euler

• Single-step: use information from the current time only.
• Multistep: use information from previous steps.
• Implicit: solve the given function to find $y_{n+1}$.
• Explicit: evaluate the given function for $y_{n+1}$.
• Time step $h$ can be constant or vary.

Explicit, single-step scheme; more accurate as $h\to0$.

ODEs - Forward Euler

We can think of time stepping as finding an approximate integral for $y’$, since we approximate $y$ and

Geometrically, FE approximates using the tangent height $f(t_n,y_n)$:

A smooth solution curve with a Forward Euler step of width $h$, using the slope $f(t_2,y(t_2))$ to approximate the next value.

FE applies to systems:

In vector notation,

Example:

Recall the Taylor series

It approximates a function by knowing $f(a)$ and derivatives at $x=a$.

Deriving FE, 1: Approximate $y’$ using a finite difference:

Since $y’(t_n)=f(t_n,y(t_n))$,

so

Deriving FE, 2: Approximate $y’$ with a simplified Taylor series. Center at $t_n$ and evaluate at $t_{n+1}$:

Dropping $h^2$ and above,

Absolute/global error at step $n$ is $\lvert y_n-y(t_n)\rvert$. We do not know $y(t_n)$, so approximate it using Taylor series.

Local truncation error is the error at each step.

Error of Forward Euler:

The Taylor series centered at $t_n$ and evaluated at $t_{n+1}$ is

Assuming exact data at $t_n$, local error is the difference between FE and the Taylor series.

Thus local error of FE is $O(h^2)$.

Trapezoidal: keep the $y’‘(t_n)$ term, approximating

Then

with $O(h^3)$ LTE. Intuition: evaluate slope at both ends and step according to average slope. Trapezoidal is implicit because it requires $f(t_{n+1},y_{n+1})$.

Explicit schemes are simpler/faster per step and less stable. Implicit schemes are more complex/expensive per step and more stable, so can use larger time steps.

Improved Euler: avoid implicitness by estimating the endpoint with FE:

with LTE $O(h^3)$.

ODEs - More Schemes

Recall multivariable Taylor expansion:

Beginning with

and using

gives

which is the trapezoidal rule.

For Improved Euler, use the FE estimate $y_{n+1}^*=y_n+hf(t_n,y_n)$ and expand $f$ in the $y$ variable:

Substitution in trapezoidal gives the explicit Improved Euler formula with $O(h^3)$ LTE.

Backward Euler: Forward Euler, but take the height at $t_{n+1}$ instead of $t_n$, making it implicit; LTE $O(h^2)$:

Remember sign of $h$ alternates [when deriving backward differences].

Runge-Kutta

RK2 / midpoint: Take an FE step to the halfway point, evaluate slope there as $k_2$, then use it to take a full step:

Equivalently,

LTE is $O(h^3)$.

Fourth-order Runge-Kutta: Evaluate $y’(t)=f(t,y)$ at intermediate positions $k_1,\ldots,k_4$ and combine to find $y_{n+1}$. Fits a quadratic; LTE $O(h^5)$.

For constant step $h$, number of steps is $O(h^{-1})$. Therefore global error is local error times all steps, so the global order is one degree less:

Backward Differentiation Formulas

BDF1, BDF2, etc.; the number indicates global-error order. Generalizes Backward Euler using previous time steps.

Derive by interpolation:

1. Fit an interpolant $p(t)$ with Lagrange polynomials to unknown $(t_{n+1},y_{n+1})$ and one or more earlier points.
2. Determine $p’(t)$ by differentiating.
3. Require endpoint slope to match: $p’(t_{n+1})=f(t_{n+1},y_{n+1})$.
4. Rearrange for BDF.

For BDF1, fit $(t_n,y_n)$ and $(t_{n+1},y_{n+1})$:

Then

so $y_{n+1}=y_n+hf(t_{n+1},y_{n+1})$.

For BDF2, fit $(t_{n-1},y_{n-1})$, $(t_n,y_n)$, $(t_{n+1},y_{n+1})$. At $t_{n+1}$,

Set equal to $f(t_{n+1},y_{n+1})$ and rearrange to obtain BDF2.

Second-order Adams-Bashforth: explicit, multistep, LTE $O(h^3)$:

Scheme summary:

ODEs - Higher-Order Systems and Stability

Order of an ODE: highest derivative appearing in it. Convert higher-order ODEs into first-order ODEs:

1. For each variable $y$ with more than a first derivative, introduce $y_i=y^{(i-1)}$ for $i=1$ to $n$.
2. Substitute new $y_i$’s into the original ODE.
3. Relate the $y_i$’s through new equations.

Unstable: error $\epsilon$ in initial conditions leads to error growing exponentially as steps $n\to\infty$.

Use the test equation $y’=-\lambda y$. Exact $y(t)=y_0e^{-\lambda t}\to0$ as $t\to\infty$. A numerical approximation should also have $y_n\to0$.

To test stability: apply a time-stepping scheme, find a closed form for solution/error behaviour, and find conditions on $h$ ensuring error approaches zero.

ODEs - Stability and Truncation Error

For Forward Euler on $y’=-\lambda y$, $\lambda>0$:

It goes to zero only when

If $y(t_0)=y_0+E_0$, then

FE is conditionally stable.

For Backward Euler,

For $\lambda>0,h>0$, both solution and error go to zero. BE is stable.

For Improved Euler,

Stability requires

which gives $0<h\lambda<2$, as in FE. It has the same stability range for this problem even though LTE is $O(h^3)$.

For linear test equations, the bound relates to $h\lvert f_y\rvert$, e.g. $\lvert h(-a)\rvert<1$ for $y’=-ay$. For nonlinear problems it depends on $f_y$ at some point. For systems it depends on eigenvalues of the Jacobian

To find local truncation error for $y_{n+1}=\mathrm{RHS}$:

1. Replace approximations on RHS with exact versions, e.g. $y_n\to y(t_n)$ and $f(t_{n+1},y_{n+1})\to y’(t_{n+1})$.
2. Taylor-expand all RHS quantities about $t_n$.
3. Taylor-expand exact $y(t_{n+1})$ to compare.
4. Compute $y(t_{n+1})-y_{n+1}$. Lowest-degree non-cancelling power of $h$ gives LTE.

For FE,

whereas

so LTE is $\frac{h^2}{2}y’‘(t_n)=O(h^2)$.

For trapezoidal,

Taylor-expanding $y’(t_{n+1})$ gives

so LTE is $O(h^3)$.

For BDF2,

Using

gives

Comparing with Taylor, the difference is $O(h^3)$.

Adaptive Time Stepping

Smaller time steps mean less error but larger computational cost. Vary $h$ with adaptive time stepping:

1. Compute approximate solutions with two schemes of different orders.
2. Estimate error by taking their difference.
3. While error $>$ tolerance, set $h:=h/2$ and recompute.
4. Estimate the error coefficient and predict a good next step $h_{\text{new}}$.
5. Repeat until end time.

At the same time, let

Then

so the dominant term matches method A’s true error.

Estimate

Assume $C$ changes slowly. Set

and solve:

Scale this back by a factor, e.g. $1/2$ or $3/4$, because the estimate will not be exact.

Stiff problem: an interesting aspect of the solution changes slowly, so we want big time steps, but a stability condition forces small time steps. Time-stepping schemes designed for this are called stiff solvers.

Fourier Transforms - Introduction

Fourier analysis:

1. Apply a Fourier transform to discrete data or a continuous function in the time/space domain to get Fourier form in the frequency domain.
2. Process/analyze in the frequency domain, e.g. sharpen an image.
3. Transform back.

For a continuous Fourier series,

$a_k,b_k$ indicate amplitude for period $2\pi/k$ or frequency $k/(2\pi)$.

Orthogonality identities over $[0,2\pi]$:

To determine $a_0$, integrate over $2\pi$:

so $a_0=(2\pi)^{-1}\int_0^{2\pi}f(t)\,dt$, the average value. Multiplying by $\cos(\ell t)$ or $\sin(\ell t)$ and integrating isolates

For $z=a+bi$, $\bar z=a-bi$, $\lvert z\rvert=\sqrt{a^2+b^2}$, and $\arg z=\operatorname{atan2}(b,a)$. Euler’s formula is $e^{i\theta}=\cos\theta+i\sin\theta$, hence

Complex Fourier form:

For $k>0$, $c_k=(a_k-ib_k)/2$, $c_{-k}=(a_k+ib_k)/2$, and $c_0=a_0$. $\lvert c_k\rvert$ is amplitude and $\arg(c_k)$ is phase. Truncate to $-M\leq k\leq M$ for an approximation.

Discrete Fourier Transform

Suppose $f_n=f(t_n)$ for $N$ points, $t_n=2\pi n/N$, with $N$ even. The inverse discrete Fourier transform approximates

The $\omega^k$ are the $N$th roots of unity. Define coefficients periodically, $F_{j+N}=F_j$.

Orthogonality:

For $k=\ell$, the sum is $N$. Otherwise, use the geometric-series identity and $\omega^{N(k-\ell)}=1$ to get zero.

Multiply the inverse relation by $\omega^{-nk}$ and sum over $n$ to isolate the DFT:

DFT properties:

1. ${F_k}$ is doubly infinite and periodic with period $N$: if $k=mN+p$, then $\omega^{-k}=\omega^{-p}$ and $F_k=F_p$.
2. Conjugate symmetry: for real data, $\overline{F_k}=F_{N-k}$, using $\overline{\omega^j}=\omega^{-j}$.

Power spectrum: plot $\lvert F_k\rvert$. A smooth bump indicates one dominant frequency; a rough bump indicates more active frequencies; symmetry indicates real data.

In matrix form $F=Mf$, where the $k$th row is $N^{-1}[1,\omega^{-k},\ldots,\omega^{-(N-1)k}]$. Orthogonality gives $M^TM=I/N$, so $M^{-1}=NM^T$ [with conjugation understood for the complex matrix]. Thus IDFT is $f=M^{-1}F$.

Fast Fourier Transform Derivation

Slow Fourier transform costs $O(N^2)$: for each $k$, sum all $N$ samples.

For the FFT, rewrite one DFT as two half-length DFTs:

Since $\omega^{-Nk/2}=(-1)^k$,

Define $g_n=(f_n+f_{n+N/2})/2$ and $h_n=(f_n-f_{n+N/2})\omega^{-n}/2$. Then $F_{2j}=G_j$ and $F_{2j+1}=H_j$, two DFTs of length $N/2$. Each stage applies even/odd splitting; bit-reverse the output to recover original indices. This works when $N$ is a power of two.

Splitting costs $O(N)$ per stage, and $N=2^m$ gives $m=\log_2N$ stages, hence $O(N\log N)$. The inverse transform has the same form, with signs reversed and the corresponding normalization.

Two-Dimensional FFT and Aliasing

For a one-dimensional grayscale image, take DFT of pixel values. For compression, throw away $F_k$ where $\lvert F_k\rvert<\text{tol}$, then reconstruct with IDFT.

For an $N\times M$ image,

Apply FFT to each of $M$ rows, costing $O(MN\log N)$, then each of $N$ columns, costing $O(NM\log M)$: total $O(MN\log(MN))$. Divide images into blocks for compression because pixels in a block often have similar data.

Aliasing: a high-frequency signal appears as low frequency due to sample spacing. Anti-aliasing filters high frequencies. Sampling rate is $f_s=N/T$ samples/second; Fourier frequency is $k/T$ cycles/second.

Sampling the continuous series at $t_n=nT/N$ gives $f_n=\sum_kc_k\omega^{nk}$. Substituting into the DFT and using orthogonality gives

Thus DFT coefficients sum Fourier-series coefficients including arbitrarily high frequencies.

Fourier Correlation Function

High frequencies outside $[-N/2,N/2]$ alias as low frequencies inside it. Once sampled, there is no way to distinguish aliased frequencies. Partial solutions:

1. Increase sampling rate.
2. Filter before sampling to remove high frequencies, e.g. optical low-pass filter.

For two real periodic data sets $y_i,z_i$, define circular correlation

The maximizing $n$ shifts $y$ into maximum dot product with $z$.

Naively trying all shifts costs $O(N^2)$. With FFT:

1. $Y=\operatorname{FFT}(y)$, $Z=\operatorname{FFT}(z)$: $O(N\log N)$.
2. $\Phi_k=\overline{Y_k}Z_k$ for $k=0,\ldots,N-1$: $O(N)$.
3. $\phi=\operatorname{IFFT}(\Phi)$: $O(N\log N)$.

The notes derive step 2 by substituting the inverse DFTs for $y_{i+n}$ and $z_i$ into the correlation DFT and applying roots-of-unity orthogonality.

Numerical PageRank

Degree: number of links leaving a node. Adjacency matrix: $G_{i,j}=1$ if link $i\to j$ exists and $0$ otherwise. Degree of node $i$ is the sum of entries in its column [under the transition-matrix orientation used here].

Heuristics: nodes with lots of incoming links, and nodes related to nodes with lots of incoming links, are important.

Random-surfer model: Start at some page and follow links at random for $K$ steps. Rank page $i$ by visits$(i)$/total visits. Problems: cycles, dead ends, and large $K$.

Define

Let $d_i=1$ if $\deg(i)=0$, otherwise $0$. Let $R$ be the number of pages and $e=(1,\ldots,1)^T$. Define

so a dead end transitions to all pages with equal chance. Define

The first term follows links without dead ends; the second teleports randomly to avoid cycles. $0\leq M_{i,j}\leq1$ and $\sum_iM_{i,j}=1$: properties of a Markov matrix.

For probability vector $p=(q_1,\ldots,q_R)^T$, $q_i\geq0$, $\sum_iq_i=1$,

is the probability of being on each page after $n+1$ transitions. It remains a probability vector because entries are nonnegative and

We estimate $p^\infty=\lim_{k\to\infty}M^kp^0$; higher probability means greater importance.

p0 = e/R
for k = 1, ...
    pk = M p(k-1)
    if max_i |pk_i - p(k-1)_i| < tol, quit

Compute PageRank once and store it; keyword searches find matching pages and rank by $p^\infty$. Exploit sparsity without building dense $M$:

Since $e^Tp^n=1$, the final term is $(1-\alpha)e/R$.

Eigenvalues/Eigenvectors Review

To incorporate other factors, use

where $v$ places extra weight on preferred sites.

An eigenvalue/eigenvector satisfies $Qx=\lambda x$, or $(\lambda I-Q)x=0$. A nonzero solution requires $\det(\lambda I-Q)=0$.

Properties of Markov matrices:

1. Every Markov matrix $Q$ has $1$ as an eigenvalue.
2. Every eigenvalue satisfies $\lvert\lambda\rvert\leq1$, so $1$ is largest.
3. $Q$ is positive if $Q_{i,j}>0$ for all $i,j$.
4. If $Q$ is positive Markov, there is only one linearly independent eigenvector for $\lambda=1$.

For (1), $Q^Te=e$ because rows of $Q^T$ sum to one. Since $\det Q=\det Q^T$, $1$ is also an eigenvalue of $Q$. For (2), choose the largest-magnitude component $x_k$ of an eigenvector of $Q^T$:

so $\lvert\lambda\rvert\leq1$.

PageRank Convergence

If $M$ is positive Markov, PageRank converges uniquely for any $p^0$. Write $p^0=\sum_ic_ix_i$, with $\lvert\lambda_1\rvert\geq\lvert\lambda_2\rvert\geq\cdots$:

$\lambda_1=1$, and the other terms vanish. Uniqueness follows because two scalar multiples of the unique eigenvector cannot both be distinct probability vectors summing to one. Convergence rate is proportional to $\lvert\lambda_2\rvert$. In PageRank, $\alpha$ approximates $\lvert\lambda_2\rvert$; e.g. $0.85^{144}\approx10^{-8}$. Smaller $\alpha$ converges faster at the cost of accuracy due to more random jumps.

Numerical Linear Algebra

Solve $Ax=b$ by:

1. Factor $PA=LU$, where $U$ is triangular: $O(n^3)$.
2. Solve $Lz=Pb$ by forward substitution: $O(n^2)$.
3. Solve $Ux=z$ by back substitution: $O(n^2)$.

LU factorization core:

for k = 1, ..., n
    for i = k+1, ..., n
        m_ik = a_ik / a_kk
        for j = k+1, ..., n
            a_ij = a_ij - m_ik a_kj

Store the multipliers inside $U$. The right-hand side is unnecessary for factoring; new $b$ values require only $O(n^2)$ triangular solves.

Forward substitution:

for i = 1, ..., n
    z_i = b_i
    for j = 1, ..., i-1
        z_i = z_i - L_ij z_j

Back substitution similarly proceeds from $i=n$ down to $1$ and divides by $U_{ii}$.

Division by zero or near zero causes large factors that magnify floating-point error. A permutation matrix is an identity matrix with rows swapped; partial pivoting swaps rows so the largest-magnitude entry in the current column is in the pivot row. The notes work a $PA=LU$ example and define FLOP as a floating-point addition, subtraction, multiplication, or division. Useful sums include $\sum_{i=1}^ni=n(n+1)/2$ and $\sum_{i=1}^ni^2=n(n+1)(2n+1)/6$.

Ignoring pivots,

so

The inverse of an elimination matrix negates its off-diagonal multiplier. Solving by explicitly computing $A^{-1}$ costs more and incurs more floating-point error: factor $A=LU$, solve $LUx=e_i$ for every identity column, then compute $x=A^{-1}b$.

Norms and Conditioning

Norm means size/magnitude:

and generally $\lVert x\rVert_p=(\sum_i\lvert x_i\rvert^p)^{1/p}$.

Properties: $\lVert x\rVert=0\Leftrightarrow x=0$; $\lVert\alpha x\rVert=\lvert\alpha\rvert\lVert x\rVert$; $\lVert x+y\rVert\leq\lVert x\rVert+\lVert y\rVert$.

Matrix norm:

Matrix-norm properties: positivity, homogeneity, triangle inequality, $\lVert Ax\rVert\leq\lVert A\rVert\lVert x\rVert$, $\lVert AB\rVert\leq\lVert A\rVert\lVert B\rVert$, and $\lVert I\rVert=1$.

For $Ax=b$, ask how perturbing $b$ or $A$ changes $x$.

Conditioning - Perturbing b

Subtract $Ax=b$: $A\delta x=\delta b$, so $\delta x=A^{-1}\delta b$. From

and

obtain

Conditioning - Perturbing A

Expanding and subtracting $Ax=b$ gives

so

Condition number $\kappa$ bounds relative change in $x$ due to relative change in $A$ or $b$. If $\kappa\approx1$, $A$ is well-conditioned; if $\kappa\gg1$, it is ill-conditioned. If unspecified, use the 2-norm. Norm equivalence: $C_1\lVert x\rVert_a\leq\lVert x\rVert_b\leq C_2\lVert x\rVert_a$.

Residual

It stands in for error in $x$ and can be reduced until $r<$ tolerance. If $x_{\text{approx}}=x+\delta x$, then $A(x+\delta x)=b-r$, so

If $Ax=\lambda x$, then $A^{-1}x=\lambda^{-1}x$: eigenvalues of $A^{-1}$ are reciprocals of those of $A$.

Exercises:

1. Prove $\kappa(A^{-1})=\kappa(A)$:

1. Prove $\kappa(A)\geq1$:

1. Decide whether the displayed diagonal/near-diagonal example matrices are ill-conditioned or well-conditioned.