Numerical Computation Systems of Equations

Gaussian Elimination

We want to solve a system of linear equations:

$$A\mathbf{x}=\mathbf{b}$$

where

$$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right],\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right],\mathbf{b}=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right].$$

---

Forward Elimination

Goal: eliminate all entries below the pivot element $a_{jj}$ in each column.

For each step $j=1,2,\dots ,n-1$:

$$m_{ij}=\frac{a_{ij}}{a_{jj}},\text{for }i=j+1,\dots ,n$$

Then update:

$$a_{ik}=a_{ik}-m_{ij}{a}_{jk},k=j+1,\dots ,n$$$$b_i=b_i-m_{ij}{b}_j$$

After this process, $A$ becomes an upper triangular matrix.

---

Back Substitution

Once $A$ is upper triangular, solve for $x_n,x_{n-1},\dots ,x_1$ by:

$$x_i=\frac{b_i-\sum _{j=i+1}^n{a}_{ij}{x}_j}{a_{ii}},i=n,n-1,\dots ,1$$

---

Result:

The final solution vector is:

$$\mathbf{x}=[x_1,x_2,\dots ,x_n{]}^T$$

obtained by forward elimination followed by back substitution.

import numpy as np

def gauss_naive(A, b):
    """
    Solve A x = b using naive Gaussian elimination (no pivoting).
    A: (n,n) ndarray, b: (n,) or (n,1)
    Returns x shape (n,)
    """
    # Ensure inputs are float type
    A = np.array(A, dtype=float)
    b = np.array(b, dtype=float).reshape(-1,)

    n = A.shape[0]
    # Augment matrix [A | b]
    Ab = np.hstack([A, b.reshape(-1,1)])

    # Forward elimination
    for j in range(n-1):
        # Check for small pivot (numerical stability check)
        if abs(Ab[j, j]) < np.finfo(float).eps:
            raise ValueError(f"Zero pivot encountered at column {j}")
            
        for i in range(j+1, n):
            mult = Ab[i, j] / Ab[j, j]
            # Row operation: Row_i = Row_i - mult * Row_j
            Ab[i, j:] = Ab[i, j:] - mult * Ab[j, j:]

    # Back substitution
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        # x[i] = (b[i] - sum(A[i, j] * x[j])) / A[i, i]
        rhs = Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:n]
        x[i] = rhs / Ab[i, i]
        
    return x

Sources of Error

Vector norm

Formulas:

• Manhattan $1-norm$:

$$\parallel \mathbf{x}{\parallel }_1=\sum _{i=1}^n|x_i|$$

• Euclidean $2-norm$:

$$\parallel \mathbf{x}{\parallel }_2=\sqrt{\sum _{i=1}^n|x_i{|}^2}$$

• Chebyshev $\mathrm{\infty }-norm$:

$$\parallel \mathbf{x}{\parallel }_{\mathrm{\infty }}=\underset{i}{max}|x_i|$$

Matrix norm

Given a vector norm $\parallel \mathbf{x}\parallel$,

the matrix norm induced by it is defined as:

$$\parallel A\parallel =\underset{\mathbf{x}\ne 0}{max}\frac{\parallel A\mathbf{x}\parallel }{\parallel \mathbf{x}\parallel }=\underset{\parallel \mathbf{x}\parallel =1}{max}\parallel A\mathbf{x}\parallel$$

Let $A=(a_{ij})\in {\mathbb{R}}^{n\times n}$, then:

• Manhattan $1-norm$:

$$\parallel A{\parallel }_1=\underset{1\le j\le n}{max}\sum _{i=1}^n|a_{ij}|$$

• Euclidean $2-norm$:

$$\parallel A{\parallel }_2=\sqrt{{\lambda }_{max}(A^{T}A)}$$

• Chebyshev $\mathrm{\infty }-norm$:

$$\parallel A{\parallel }_{\mathrm{\infty }}=\underset{1\le i\le n}{max}\sum _{j=1}^n|a_{ij}|$$

---

Matrix Norm Properties

For matrices $A,B$, vector $\mathbf{x}$, and scalar $c$:

1. $|A|>0$ if $A\ne 0$
2. $|cA|=|c|\cdot |A|$
3. $|A+B|\le |A|+|B|$
4. $|AB|\le |A|\cdot |B|$
5. $|A\mathbf{x}|\le |A|\cdot |\mathbf{x}|$

Condition Number

$$\kappa (A)=\frac{M}{m}$$

where

$$M=\underset{\mathbf{x}\ne 0}{max}\frac{\parallel A\mathbf{x}\parallel }{\parallel \mathbf{x}\parallel }=\parallel A\parallel \text{and}m=\underset{\mathbf{x}\ne 0}{min}\frac{\parallel A\mathbf{x}\parallel }{\parallel \mathbf{x}\parallel }=\frac{1}{\parallel A^{-1}\parallel }$$$$\boxed{{\displaystyle \kappa (A)=\parallel A\parallel \cdot \parallel A^{-1}\parallel }}$$

Meaning:

• $\kappa (A)$ measures how sensitive the solution of $A\mathbf{x}=\mathbf{b}$ is to errors.
• If $\kappa (A)$ is large $\to$ the matrix is ill-conditioned (unstable).
• If $\kappa (A)\approx 1$ $\to$ the matrix is well-conditioned (stable).

The condition number is a good indicator of how close is a matrix to be singular. The large the condition number the closer we are to singularity.

It is also very useful in assessing the accuracy of solutions to linear systems.

In practice we don’t really calculate the condition number, it is merely estimated, to perhaps within an order of magnitude.

Error Magnification

Let $x_a$ be an approximate solution of the linear system $Ax=b$. The residual is the vector $r=b-A{x}_a$.
The backward error is the norm of the residual $\parallel b-A{x}_a{\parallel }_{\mathrm{\infty }}$ The forward error is $\parallel x-x_a{\parallel }_{\mathrm{\infty }}$.

The relative backward error is:

$$\frac{\parallel r{\parallel }_{\mathrm{\infty }}}{\parallel b{\parallel }_{\mathrm{\infty }}},$$

and the relative forward error is:

$$\frac{\parallel x-x_a{\parallel }_{\mathrm{\infty }}}{\parallel x{\parallel }_{\mathrm{\infty }}}.$$

The error magnification factor is:

$$\text{error magnification factor}=\frac{\text{relative forward error}}{\text{relative backward error}}=\frac{\frac{\parallel x-x_a{\parallel }_{\mathrm{\infty }}}{\parallel x{\parallel }_{\mathrm{\infty }}}}{\frac{\parallel r{\parallel }_{\mathrm{\infty }}}{\parallel b{\parallel }_{\mathrm{\infty }}}}.$$

Partial pivoting

In Gaussian elimination, if a pivot element $a_{kk}^{(k)}$ is small compared to an element $a_{jk}^{(k)}$ below, the multiplier

$$m_{jk}=\frac{a_{jk}^{(k)}}{a_{kk}^{(k)}}$$

will be large, resulting in round-off error.

Partial pivoting finds the smallest $p\ge k$ such that

$$\left|a_{pk}^{(k)}\right|=\underset{k\le i\le n}{max}\left|a_{ik}^{(k)}\right|$$

and interchanges the row $(R_k)\leftrightarrow (R_p)$.

Scaled partial pivoting

If there are large variations in magnitude of the elements within a row, scaled partial pivoting can be used.

Define a scale factor $s_i$ for each row:

$$s_i=\underset{1\le j\le n}{max}\left|a_{ij}\right|$$

At step $i$, find $p$ such that

$$\frac{\left|a_{pi}\right|}{s_p}=\underset{i\le k\le n}{max}\frac{\left|a_{kj}\right|}{s_k}$$

and interchange the row $(R_i)\leftrightarrow (R_p)$.

The LU Factorization

DEFINITION: An $m\times n$ matrix $L$ is lower triangular if its entries satisfy $l_{ij}=0$ for $i<j$. An $m\times n$ matrix $U$ is upper triangular if its entries satisfy $u_{ij}=0$ for $i>j$.

EXAMPLE:

$$\left[\begin{array}{cc}1& 1\\ 3& -4\end{array}\right]\to R_2-3R_1\to \left[\begin{array}{cc}1& 1\\ 0& -7\end{array}\right]=U$$

Define $L$ to be the $2\times 2$ lower triangular matrix with 1's on the main diagonal and the multiplier 3 in the (2,1) location:

$$\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right].$$

Check that

$$LU=\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]\left[\begin{array}{cc}1& 1\\ 0& -7\end{array}\right]=\left[\begin{array}{cc}1& 1\\ 3& -4\end{array}\right]=A$$

Derivation of L U Matrix:

$$Ax=b$$$$M_1{M}_2\cdots M_{n}A=U$$$$A=M_n^{-1}\cdots M_2^{-1}{M}_1^{-1}U$$$$L=M_n^{-1}\cdots M_2^{-1}{M}_1^{-1}$$

EXAMPLE:

Find the LU factorization of

$$\left[\begin{array}{ccc}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]$$

SOLUTION:

$$\left[\begin{array}{ccc}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]\to R_2-\mathbf{2}{R}_1\to \left[\begin{array}{ccc}1& 2& -1\\ 0& -3& 0\\ -3& 1& 1\end{array}\right]$$$$\to R_3-(\mathbf{-}\mathbf{3})R_1\to \left[\begin{array}{ccc}1& 2& -1\\ 0& -3& 0\\ 0& 7& -2\end{array}\right]$$$$\to R_3-(\mathbf{-}\frac{\mathbf{7}}{\mathbf{3}})R_2\to \left[\begin{array}{ccc}1& 2& -1\\ 0& -3& 0\\ 0& 0& -2\end{array}\right]=U$$$$L=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ -3& -\frac{7}{3}& 1\end{array}\right],LU=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ -3& -\frac{7}{3}& 1\end{array}\right]\left[\begin{array}{ccc}1& 2& -1\\ 0& -3& 0\\ 0& 0& -2\end{array}\right]=\left[\begin{array}{ccc}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]=A$$

---

• FACT 1:

$$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\to \left[\begin{array}{ccc}1& 0& 0\\ -c& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$$$$=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}-c{a}_{11}& a_{22}-c{a}_{12}& a_{23}-c{a}_{13}\\ a_{31}& a_{32}& a_{33}\end{array}\right].$$

• FACT 2:

$${\left[\begin{array}{ccc}1& 0& 0\\ -c& 1& 0\\ 0& 0& 1\end{array}\right]}^{-1}=\left[\begin{array}{ccc}1& 0& 0\\ c& 1& 0\\ 0& 0& 1\end{array}\right].$$

• FACT 3:

$$\left[\begin{array}{ccc}1& 0& 0\\ c_1& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ c_2& 0& 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& c_3& 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ c_1& 1& 0\\ c_2& c_3& 1\end{array}\right]$$

Forward and Back substitution

Once L and U are known, the problem $Ax=b$ can be written as $LUx=b$. Define a new “auxiliary” vector $c=Ux$. Then back substitution is a two-step procedure: a) Solve $Lc=b$ for $c$. b) Solve $Ux=c$ for $x$.

EXAMPLE:

Solve $\left[\begin{array}{cc}1& 1\\ 3& -4\end{array}\right]x=\left[\begin{array}{c}3\\ 2\end{array}\right]$, using LU factorization.

$$\left[\begin{array}{cc}1& 1\\ 3& -4\end{array}\right]=LU=\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]\left[\begin{array}{cc}1& 1\\ 0& -7\end{array}\right]$$

(a) Solve $\left[\begin{array}{cc}1& 0\\ 3& 1\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}3\\ 2\end{array}\right]$ by forward substitution: $c_1=3,3c_1+c_2=2\Rightarrow c_2=-7$. (b) Solve $\left[\begin{array}{cc}1& 1\\ 0& -7\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}3\\ -7\end{array}\right]$ by back substitution: $x_2=1,x_1+x_2=3\Rightarrow x_1=2$.

EXAMPLE:

Solve $\left[\begin{array}{ccc}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]x=\left[\begin{array}{c}3\\ 3\\ -6\end{array}\right]$, using LU factorization.

$$\left[\begin{array}{ccc}1& 2& -1\\ 2& 1& -2\\ -3& 1& 1\end{array}\right]=LU=\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ -3& -\frac{7}{3}& 1\end{array}\right]\left[\begin{array}{ccc}1& 2& -1\\ 0& -3& 0\\ 0& 0& -2\end{array}\right]$$

(a) Solve $\left[\begin{array}{ccc}1& 0& 0\\ 2& 1& 0\\ -3& -\frac{7}{3}& 1\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}3\\ 3\\ -6\end{array}\right]$ by forward substitution: $c_1=3,2c_1+c_2=3\Rightarrow c_2=-3,-3c_1-\frac{7}{3}{c}_2+c_3=-6\Rightarrow c_3=-4$. (b) Solve $\left[\begin{array}{ccc}1& 2& -1\\ 0& -3& 0\\ 0& 0& -2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}3\\ -3\\ -4\end{array}\right]$ by back substitution: $-2x_3=-4\Rightarrow x_3=2,-3x_2=-3\Rightarrow x_2=1,x_1+2x_2-x_3=-4\Rightarrow x_1=3$.

PA=LU factorization

EXAMPLE: Prove that $A=\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]$ does not have an $LU$ factorization.

SOLUTION:

The factorization must have the form

$$\left[\begin{array}{cc}0& 1\\ 1& 1\end{array}\right]=\left[\begin{array}{cc}1& 0\\ a& 1\end{array}\right]\left[\begin{array}{cc}b& c\\ 0& d\end{array}\right]=\left[\begin{array}{cc}b& c\\ ab& ac+d\end{array}\right].$$

Equating coefficients yields $b=0$ and $ab=1$, a contradiction.

PA = LU factorization is the matrix formulation of elimination with partial pivoting. PA = LU factorization is simply the LU factorization of a row-exchanged version of A.

PA=LU factorization

DEFINITION:

A permutation matrix is an $n\times n$ matrix consisting of all zeros, except for a single 1 in every row and column.

EXAMPLE:

$$\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right],\left[\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right],\dots \dots$$

THEOREM:

Fundamental Theorem of Permutation Matrices. Let $P$ be the $n\times n$ permutation matrix formed by a particular set of row exchanges applied to the identity matrix. Then, for any $n\times n$ matrix $A$, $PA$ is the matrix obtained by applying exactly the same set of row exchanges to $A$.

EXAMPLE:

$$\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right]\left[\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right]=\left[\begin{array}{ccc}a& b& c\\ g& h& i\\ d& e& f\end{array}\right]$$

EXAMPLE: Find the PA=LU factorization of the matrix

$$A=\left[\begin{array}{ccc}2& 1& 5\\ 4& 4& -4\\ 1& 3& 1\end{array}\right]$$$$\left[\begin{array}{ccc}2& 1& 5\\ 4& 4& -4\\ 1& 3& 1\end{array}\right]\to P_1=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right]\to \left[\begin{array}{ccc}4& 4& -4\\ 2& 1& 5\\ 1& 3& 1\end{array}\right].$$$$R_1\leftrightarrow R_2$$$$\left[\begin{array}{ccc}4& 4& -4\\ 2& 1& 5\\ 1& 3& 1\end{array}\right]\to R_2-\frac{\mathbf{1}}{\mathbf{2}}{R}_1\to \left[\begin{array}{ccc}4& 4& -4\\ \frac{\mathbf{1}}{\mathbf{2}}& -1& 7\\ 1& 3& 1\end{array}\right]\to R_3-\frac{\mathbf{1}}{\mathbf{4}}{R}_1\to \left[\begin{array}{ccc}4& 4& -4\\ \frac{\mathbf{1}}{\mathbf{2}}& -1& 7\\ \frac{\mathbf{1}}{\mathbf{4}}& 2& 2\end{array}\right]$$$$\left[\begin{array}{ccc}4& 4& -4\\ \frac{\mathbf{1}}{\mathbf{2}}& -1& 7\\ \frac{\mathbf{1}}{\mathbf{4}}& 2& 2\end{array}\right]\to P_2=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right]\to \left[\begin{array}{ccc}4& 4& -4\\ \frac{\mathbf{1}}{\mathbf{4}}& 2& 2\\ \frac{\mathbf{1}}{\mathbf{2}}& -1& 7\end{array}\right].$$$$R_2\leftrightarrow R_3$$$$\left[\begin{array}{ccc}4& 4& -4\\ \frac{\mathbf{1}}{\mathbf{4}}& 2& 2\\ \frac{\mathbf{1}}{\mathbf{2}}& -1& 7\end{array}\right]\to R_3-(\mathbf{-}\frac{\mathbf{1}}{\mathbf{2}})R_2\to \left[\begin{array}{ccc}4& 4& -4\\ \frac{\mathbf{1}}{\mathbf{4}}& 2& 2\\ \frac{\mathbf{1}}{\mathbf{2}}& \mathbf{-}\frac{\mathbf{1}}{\mathbf{2}}& 8\end{array}\right]$$$$P_2{P}_1=P\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right]\left[\begin{array}{ccc}2& 1& 5\\ 4& 4& -4\\ 1& 3& 1\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ \frac{\mathbf{1}}{\mathbf{4}}& 1& 0\\ \frac{\mathbf{1}}{\mathbf{2}}& \mathbf{-}\frac{\mathbf{1}}{\mathbf{2}}& 1\end{array}\right]\left[\begin{array}{ccc}4& 4& -4\\ 0& 2& 2\\ 0& 0& 8\end{array}\right]$$$$PA=LU$$

Forward and back substitution for PA=LU factorization:

$$Ax=b\Rightarrow PAx=Pb\Rightarrow LUx=Pb$$

Solve

1. $Lc=Pb$ for $c$.
2. $Ux=c$ for $x$.

EXAMPLE: Use the PA=LU factorization to solve the system $Ax=b$, where

$$A=\left[\begin{array}{ccc}2& 1& 5\\ 4& 4& -4\\ 1& 3& 1\end{array}\right],b=\left[\begin{array}{c}5\\ 0\\ 6\end{array}\right]$$

SOLUTION:

1. $Lc=Pb$:

$$\left[\begin{array}{ccc}1& 0& 0\\ \frac{\mathbf{1}}{\mathbf{4}}& 1& 0\\ \frac{\mathbf{1}}{\mathbf{2}}& \mathbf{-}\frac{\mathbf{1}}{\mathbf{2}}& 1\end{array}\right]\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right]\left[\begin{array}{c}5\\ 0\\ 6\end{array}\right]=\left[\begin{array}{c}0\\ 6\\ 5\end{array}\right].$$$$c_1=0$$$$\frac{1}{4}{c}_1+c_2=6\Rightarrow c_2=6$$$$\frac{1}{2}{c}_1-\frac{1}{2}{c}_2+c_3=5\Rightarrow c_3=8$$

2. $Ux=c$:

$$\left[\begin{array}{ccc}4& 4& -4\\ 0& 2& 2\\ 0& 0& 8\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}0\\ 6\\ 8\end{array}\right].$$$$8x_3=8\Rightarrow x_3=1$$$$2x_2+2x_3=6\Rightarrow x_2=2$$$$4x_1+4x_2-4x_3=0\Rightarrow x_1=-1.$$

Therefore, the solution is $x=[-1,2,1]$.

Operation Count

Operation include: $+,-,\times ,÷$

LEMMA:

For any positive integer $n$,

(a) $1+2+3+4+\cdots +n=\frac{n(n+1)}{2}$
(b) $1^2+2^2+3^2+4^2+\cdots +n^2=\frac{n(n+1)(2n+1)}{6}$.

Operations count for Gaussian Elimination:

• For $j=1$,

  • $mult=\frac{a(i,j)}{a(j,j)}$ First a multiplier is computed for each row below the first row. This requires $(n-1)$ multiplies.
  • $a(i,k)=a(i,k)-mult\times a(j,k)$; $b(i)=b(i)-mult\times b(j)$ Then in each row below row 1 the algorithm performs $n$ multiplies and $n$ additions . Thus, there is a total of $(n-1)+(n-1)\times 2\times n=2n^2-n-1$ operations for this step of Gaussian Elimination.

• For $j=2$, we zero out the column below $a_{22}$.

  • There are $n-2$ rows below this pivot, so this takes $2(n-1{)}^2-(n-1)-1$ operations.

• For $j=3$, we would have $2(n-2{)}^2-(n-2)-1$ operations.

• ...

To complete Gaussian Elimination, it will take $I_n$ flops, where

$$I_n=2\sum _{j=1}^n{j}^2-\sum _{j=1}^{n}j-\sum _{j=1}^{n}1$$

Based on Lemma above:

$$\begin{array}{rl}{I}_n& =2\sum _{j=1}^n{j}^2-\sum _{j=1}^{n}j-\sum _{j=1}^{n}1\\ & =\frac{2n(n+1)(2n+1)}{6}-\frac{n(n+1)}{2}-n\\ & =\frac{2}{3}{n}^3+\frac{1}{2}{n}^2-\frac{2}{3}n.\end{array}$$

Counting operation yields

$$1+3+5+\cdots +(2n-1)=\sum _{i=1}^n(2i-1)=2\sum _{i=1}^{n}i-\sum _{i=1}^{n}1=\frac{2n(n+1)}{2}-n=n^2$$

• On large $n$, where lower powers of $n$ become negligible by comparison. In this case, if we ignore the lower order terms in the expressions for the number of multiplication/divisions, we find that elimination takes on the order of $\frac{2}{3}{n}^3$ operations and that back substitution takes on the order of $n^2$.
• We will often use the shorthand terminology of "big-O" to mean "on the order of," saying that elimination is an $O(n^3)$ algorithm and that back substitution is $O(n^2)$. Overall, Gaussian elimination takes $\frac{2}{3}{n}^3+n^2\approx \frac{2}{3}{n}^3$ operations.
• In other words, for large $n$, the lower order terms in the complexity count will not have a large effect on the estimate for running time of the algorithm and can be ignored if only an estimated time is required.

EXAMPLE Estimate the time required to carry out back substitution on a system of 500 equations in 500 unknowns, on a computer where elimination takes 1 second.

SOLUTION

$$\frac{{500}^2}{2(500{)}^3/3}=\frac{t}{1}\Rightarrow t=\frac{3}{2(500)}=0.003\text{ sec.}$$

EXAMPLE On a particular computer, back substitution of a $5000\times 5000$ triangular matrix takes 0.1 seconds. Estimate the time needed to solve a general system of 3000 equations in 3000 unknowns by Gaussian elimination.

SOLUTION${5000}^2$ operations takes 0.1 seconds. A general system of 3000 equations in 3000 unknowns by Gaussian elimination needs $\frac{2}{3}(3000{)}^3$ operations.

$$\frac{2(3000{)}^3/3}{{5000}^2}=\frac{t}{0.1}\Rightarrow t\approx 72\text{ sec.}$$

Operations count for LU factorization:

Now, suppose that we need to solve a number of different problems with the same A and different b. That is, we are presented with the set of problems

$$\begin{array}{rl}Ax& =b_1\\ Ax& =b_2\\ & ⋮\\ Ax& =b_k\end{array}$$

with various right-hand side vectors $b_i$. Classical Gaussian elimination will require approximately $2k{n}^3/3$ operations. With the LU approach, the right-hand-side $b$ doesn't enter the calculations until the elimination (the $A=LU$ factorization) is finished. We can solve the previous set of equations with only one elimination, followed by two back substitutions ($Lc=b,Ux=c$) for each new $b$. $Lc=b,Ux=c$ back substitutions $2n^2$ The approximate number of operations with the LU approach is, therefore, $\frac{2n^3}{3}+2k{n}^2$.

EXAMPLE Assume that it takes one second to factorize the $3000\times 3000$ matrix $A$ into $A=LU$. How many problems $Ax=b_1,\dots ,Ax=b_k$ can be solved in the next second?

SOLUTION

$$\frac{2(3000{)}^3}{3}=2k(3000{)}^2\Rightarrow k=1000$$