1 of 90

Numerical Linear Algebra

CS 740

Instructor: Ajit Rajwade

1

2 of 90

What is numerical linear algebra?

  • Branch of mathematics dealing with operations on vectors and matrices.

  • Vast applications in: signal and image processing, computer graphics, computer vision, statistics, machine learning and data mining, economics and finance, many other fields.

2

3 of 90

What is a matrix?

  • A matrix is a 2D array of size m x n containing values defined over a “field”.
  • A “field” is a set of numbers – most common being “real” or “complex”.

3

Examples:

Square matrix: m = n

4 of 90

Row and Column Vectors

  • A matrix having only one row is called as a row-vector.

  • A matrix having only one column is called as a column-vector.

4

5 of 90

Some notation

  • The space of real numbers in n dimensions is denoted as Rn.
  • The space of complex numbers in n dimensions is denoted as Cn.
  • If A is a real/complex matrix of size m x n, we say that

5

6 of 90

Some notation

  • By convention, scalars are denoted in non-bold italic font.
  • Vectors and matrices are denoted in bold font, with matrices in upper case and vectors in lower case.

6

7 of 90

Common types of matrices

  • Diagonal matrix: a matrix with all zeros except along the diagonal, i.e. for every entry where i=j.

  • Identity matrix – square diagonal matrix with value of 1 along the diagonal – commonly denoted as Inxn.
  • Symmetric matrix: square matrix which is equal to its transpose.
  • We will see many more types.

7

8 of 90

Common matrix operations

  • Addition/subtraction
  • Scalar multiplication
  • Multiplication of a row vector with a matrix.

  • Multiplication of a matrix with a column vector.

8

9 of 90

Common matrix operations: multiplication

  • Multiplication of matrix A having size m1 x n1, with matrix B shaving size m2 x n2 (possible only when n1 = m2) produces a matrix of size m1 x n2.

9

Note:

ABBA

Matrix multiplication is non-commutative (or does not commute).

10 of 90

Common matrix operations: multiplication

  • Matrix multiplication, while not commutative, is associative, i.e. A(BC) = (AB)C for any matrices A, B, C of appropriate dimension.

  • Multiplication of a row vector (1 x n matrix) and a column vector (n x 1 matrix) produces a scalar (1 x 1 matrix!). Such a multiplication is called the inner product (dot product) of two vectors.

  • Multiplication of a column vector (n x 1 matrix) and a row vector (1 x n matrix) produces an n x n matrix. Such a multiplication is called the outer product of two vectors!

10

11 of 90

Common matrix operations

  • Transpose of a m x n real matrix A is the matrix of size n x m – denoted as AT or A’ – such that the rows of AT are columns of A, and columns of AT are rows of A.

11

12 of 90

Common matrix operations

  • For a complex matrix A, the transpose is defined slightly differently – it is the conjugate transpose and denoted A* or AH :

12

13 of 90

Common matrix operations

  • Trace of a square matrix = sum of values along the diagonal.
  • Matrix multiplication is not commutative, but the trace of the product of two matrices does not depend upon the order of the matrices in the product: trace(AB) = trace(BA).
  • Trace of a matrix equals that of its transpose.

13

14 of 90

Common matrix operations

  • Determinant of a square matrix is a scalar quantity defined recursively as follows:

14

Matrix of size n x n

Matrix of size (n-1) x (n-1) created by dropping the i-th row and j-th column of A

Element from the i-th row and j-th column of A

This formula is true for any integer i. Note that replacing i by some i’ is equivalent to swapping rows i and i’. A swap of any two rows of a matrix causes the sign of its determinant to flip. This formula here essentially involves exchanging the i-th row and the first row.

15 of 90

Matrix Inverse

  • A square matrix A is called invertible or non-singular if there exists a matrix B such that AB = I where I is the identity matrix.
  • Matrix B is called the inverse of A. Also, A is the inverse of B.
  • If a matrix has an inverse, it is unique (just like if a real number has a reciprocal, it is unique).
  • If A has no inverse, it is said to be non-invertible or singular. Such a matrix has determinant 0.
  • An invertible matrix (and its inverse) have strictly non-zero determinants.

15

16 of 90

Rank of a square matrix: linearly (in)dependent vectors

  • Consider a matrix A of size n x n. Let ai be the i-th column vector of A.
  • The set of n vectors {a1, a2,…, an} is said to be linearly dependent if:

Otherwise the vectors are said to be linearly independent.

16

  • u and v are linearly independent (they define a plane)
  • u, v, w are linearly dependent as all three are coplanar
  • u and j are anti-parallel and hence dependent
  • u, v, k are linearly independent

17 of 90

Rank of a square matrix: linearly (in)dependent vectors

  • If the vectors are linearly dependent, then there exists a vector that can be expressed as a linear combination of the other vectors:

  • If two vectors are exactly the same or proportional to each other, then the set of vectors is necessarily linearly dependent.

17

18 of 90

Rank of a square matrix: linearly (in)dependent vectors

  • The set of n vectors {a1, a2,…, an} will be linearly dependent if any one or more of these n vectors equal the zero vector.
  • Reason: Suppose ai = 0. Then 0 a1 + 0 a2 + … 1 ai + …0an-1 + 0an = 0. The weighting coefficients are not all zero.

18

19 of 90

Rank of a square matrix

  • If the column vectors of a matrix are linearly dependent, then it can shown that the row vectors will also be linearly dependent (and conversely), and the matrix will have determinant 0.
  • The rank of a matrix is equal to the number of linearly independent column vectors (also equal to the number of linearly independent row vectors).
  • If all n rows (equivalently all n columns) of a matrix of size n x n are linearly independent, then it is said to have full rank. Its determinant will be strictly non-zero, and hence it will be invertible.

19

20 of 90

Rank of a square matrix

  • If the rank of a matrix of size n x n is less than n, it is said to be a rank-deficient (or low-rank) matrix. Such a matrix has determinant 0, and hence is non-invertible.
  • Examples:

20

Non-invertible

Invertible

21 of 90

What about rank of a non-square matrix?

  • In this case, we define row-rank and column-rank.
  • Row-rank = number of linearly independent row vectors.
  • Column-rank = number of linearly independent column vectors.
  • It turns out that for any matrix, row-rank = column-rank (we will prove this identity later).
  • An m x n matrix has full row-rank if its row-rank = m, and full column-rank if its column-rank = n.

21

22 of 90

Multiplying a matrix with a vector – Geometric meaning

  • A point in 2D is represented as a column vector v with 2 elements.
  • A 2 x 2 matrix M multiplied with a 2 x 1 column vector v induces a geometric transformation on the point underlying the vector. The transformed point is Mv.
  • A shape (polygon) can be represented by an ordered sequence of k points – which can be represented as a 2 x k matrix Vk.

22

23 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Each point of the transformed shape is represented as a column of the matrix MVk.
  • Example 1:

23

Reflection across Y axis

Reflection across X axis

These matrices are called reflection matrices.

24 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Example 2:

24

  • This matrix is called a scaling matrix.
  • Scaling by a factor of ‘a’ along X axis, and a factor of ‘b’ along the Y axis.
  • If a > 1, it induces expansion along X axis, otherwise if 0 < a < 1, it causes shrinking along X axis.
  • If b > 1, it induces expansion along Y axis, otherwise if 0 < b < 1, it causes shrinking along Y axis.
  • If a= b, it is called isotropic scaling, otherwise it is called anisotropic scaling.

25 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Example 3:

25

This matrix is called as a rotation matrix – more specifically it is a rotation matrix for anti-clockwise rotation through angle Θ about the origin. The rotation matrix for clockwise rotation through angle Θ about the origin is given as follows:

26 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Any 2 x 2 matrix induces a geometric transformation – called as a linear transformation.
  • Rotation, scaling and reflection transformations are all special cases of a linear transformation.
  • Why linear transformation?

26

27 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Any transformation of vectors – denoted as T(.) – is said to be linear if:

T(ax + by) = aT(x) + bT(y) for any vector x and y, and any scalars a, b (for now, we assume all these are defined over real numbers).

  • Notice that:

A(ax + by) = aAx + bAy

27

28 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Consider the matrix

  • It can viewed as a matrix that

transforms the points (0,0),

(0,1),(1,0),(1,1) into (0,0),

(c,d),(a+c,b+d),(a,b).

28

The vectors represented by a 2-by-2 matrix correspond to the sides of a unit square transformed into a parallelogram.

29 of 90

Multiplying a matrix with a vector – Geometric meaning

  • What happens if the transformation matrix M is singular?
  • It will shrink shapes into a single line, or even a single point. In other words, it will transform the original shape into one with area 0.
  • There is an interesting geometric property of the determinant:

29

30 of 90

Multiplying a matrix with a vector – Geometric meaning

  • Look at the rotation matrix. Rotation is an operation that does not change the area of the shape (this is called as a rigid transformation).
  • Notice that the determinant of the rotation matrix is 1.
  • For a reflection matrix, the determinant is -1: notice that reflection flips the directions of the vectors!

30

31 of 90

Multiple geometric transformations

  • Multiple geometric transformations can be performed by sequential multiplication of their matrices with the original vector(s).
  • This is called composition of transformations.
  • Example:

  • Matrix multiplication generally does not commute – hence the changing the order of the matrices will change the effective geometric transformation!

31

32 of 90

Inverse geometric transformations

  • Some geometric transformations are inverses of one another. Their corresponding transformation matrices also turn out to be inverses of each other!
  • Example: anti-clockwise and clockwise rotation about the origin through angle theta

  • Any more examples?

32

33 of 90

More geometric transformations

  • We have seen that 2 x 2 matrices represent linear geometric transformation in 2D space. Likewise, 3 x 3 matrices represent linear geometric transformation in 3D space.
  • Bear in mind: this geometric interpretation is also valid for matrices of arbitrary dimension (> 3), but we do not pursue it here further! It is very useful in pattern recognition and statistics.
  • Likewise, n x n matrices represent linear transformations in n-D space.

33

34 of 90

Motion Models

  • Let us consider points (x1,y1) and (x2,y2).
  • Translation:

  • Rotation about point (0,0) anti-clockwise through angle θ

34

2D Rotation matrix (orthonormal matrix)

35 of 90

Motion Models

  • Rotation about point (xc,yc) anti-clockwise through angle θ

35

-Perform translation such that (xc,yc) coincides with the origin (0,0).

-Rotate about the new origin.

-Translate back.

-The extra ones (third row) are called homogeneous coordinates – they facilitate using matrix multiplication to represent translations.

36 of 90

Motion Models

  • Rotation and translation:

  • Affine transformation: (rotation, scaling and shearing) besides translation

36

Assumption: the 2 x 2 sub-matrix A is NOT rank deficient, otherwise it will transform two-dimensional figures into a line or a point

37 of 90

Systems of Linear Equations

  • Equations of the following type arise in many areas of scientific computing:

37

Unknowns are x1,x2,…,xn – the coefficients {bi} and {aij} are all known.

38 of 90

Systems of Linear Equations

  • Equations of the following type arise in many areas of scientific computing:

38

Unknown: x (size n x 1),

Known: A (size n x n),b (size n x 1).

There are systems of linear equations where A has size m x n (m not equal to n), but we will take up those cases later!

39 of 90

Systems of Linear Equations

  • Not all such systems of equations have a solution.
  • A solution exists if and only if A is invertible. In this case, the solution is unique and given by x = A-1b.
  • If A is singular, then there may exist no solution or more than one solutions. This depends on the vector b.

39

40 of 90

Systems of Linear Equations

  • Example 1: Unique solution

  • Example 2: No solution

  • Example 3: Multiple solutions

40

41 of 90

Systems of Linear Equations: Geometric Interpretation

  • Let’s stick to 2D (i.e. 2 unknowns) for now.
  • The system of equations has the form:

  • Geometrically, we are dealing with two lines: ax+by = e, cx+dy=f.
  • If these lines intersect, they either do so in one and only one point, or the two lines are coincident.
  • First case: unique solution, second case: infinitely many solutions all lying on a line.

41

42 of 90

Systems of Linear Equations: Geometric Interpretation

  • If these lines do not intersect, they will be parallel. In that case, there is no solution.
  • Try to interpret the earlier three examples in this geometric context.

42

43 of 90

Systems of Linear Equations: Geometric Interpretation

  • The same geometric treatment extends to 3D (i.e. 3 unknowns).
  • In this case, each of the 3 equations has the form ax+by+cz=d, which is the equation of a plane.
  • If the three planes intersect, they either do so in a point (unique solution), or all three planes are coincident (infinitely many solutions all lying on a plane), or they intersect in a line (infinitely many solutions all lying on a line).
  • If the 3 planes do not intersect, they are parallel (no solution).

43

44 of 90

Systems of Linear Equations: Geometric Interpretation

  • The same geometric treatment extends to n-D (i.e. n unknowns, n = 4 or more).
  • In this case, each of the n equations is a hyper-plane.
  • If the n hyper-planes intersect, may do so in a point (unique solution), or in an r-dimensional “flat” where 1 <= r <= n -1 (infinitely many solutions).
  • If the n hyper-planes do not intersect, they are parallel (no solution).

44

45 of 90

Homogeneous system of equations

  • The system Ax = 0 is called a homogenous system – a system in which the constant terms are all zero.
  • Such a system always has at least one solution, i.e. x = 0 – often called the trivial solution.
  • If A has full rank, the only solution is the trivial solution.
  • In particular note that if A is not invertible, the system Ax = 0 will have many solutions x which are non-zero. Can you construct such an example yourself?

45

46 of 90

Homogeneous system of equations

  • If A is singular, this system has infinitely many solutions.
  • If u and v are two solutions to this system, then u+v and ku are also solutions (where k is a scalar).
  • All such solutions constitute what is called the null-space or kernel of matrix A.

46

47 of 90

Homogeneous system of equations

  • Example:

  • For a full-rank matrix, the null-space contains the vector 0 only and it is called the trivial nullspace.
  • For a low-rank matrix, the null-space will contain other vectors (besides 0) as well: it is called the non-trivial nullspace.

  • How do we solve Ax = 0 for the other non-zero vectors? We will study later!

47

The null-space of this matrix consists of all points lying on the line x=0,y+z=0

48 of 90

Range of a matrix

  • All those vectors which can be expressed as linear combinations of the columns of A constitute the range (also called column-space) of A.
  • The column vectors are said to span the range of A.
  • For full rank matrix A of size n x n, the range is the entire Rn.

48

49 of 90

Example: Range

  • Any vector with 3 elements can be expressed as a linear combination of the columns of:

  • But some vectors cannot be expressed as linear combinations of the columns of:

49

50 of 90

Recap: singular matrices

  • An n x n matrix A is singular if any one of the following equivalent conditions is true:
  • Det(A) = 0
  • Rank(A) < n
  • Inverse of A does not exist
  • There exists some non-zero vector y such that Ay = 0 (which means that the columns of A are not linearly independent).

50

51 of 90

System of linear equations with m≠n

  • Generally, if m (number of equations) < n (number of unknowns), the system has infinitely many solutions. It is called an under-determined system.
  • Generally, if m > n, the system has no solution. It is called an over-determined system.

51

52 of 90

System of linear equations with m≠n

  • Example of under-determined system:

  • Example of over-determined system:

52

53 of 90

Solving an over-determined system

53

  • The inverse of A is not defined as A is not a square matrix. But ATA is a square matrix of size n x n, and it will have a well-defined inverse if A has full column-rank.

  • (ATA)-1AT is called the pseudo-inverse of A and is often denoted as A.

  • The x thus obtained is not a true solution – it will not exactly satisfy all (or any of!) the equations. The solution derivation is on the next slide.

54 of 90

Solving an over-determined system

54

See the board-work to understand how the derivatives were obtained.

55 of 90

Solving an over-determined system

55

  • Over-determined systems can be solved in MATLAB using the following commands :

x = inv(A’*A)*A’*b; % slower, less accurate

x = A\b; % this is the back-slash operator and it implements the

% pseudo-inverse of A efficiently and accurately

% In case of either method, the MATLAB system will produce a %warning if the matrix is low-rank or ill-conditioned (i.e. close to %being a low-rank matrix)

56 of 90

Application of over-determined systems

  • Image alignment is a popular application in image processing.
  • Consider two images – where one image is a linear transformed version of the other.

56

57 of 90

Application of over-determined systems

  • Let’s say I mark out some N salient points {(x1i,y1i)} in the left image, and their corresponding points {(x2i,y2i)} in the right image – see figure on next slide.
  • I can write the equation: P2 = A P1 where P1 is a 2 X N matrix, whose i-th row contains the x and y coordinates of the i-th point in the first set (likewise for P2). A is the unknown linear transformation.
  • The solution for A: .

57

58 of 90

58

59 of 90

Application of over-determined systems

  • This is a least squares solution.
  • Note: In practice, the equation P2 = A P1 will never be satisfied exactly. We are humans after all, and there will always be some usually small error in marking the point positions in P2 (we will assume P1 contains no errors – it provides a reference set of points after all!)
  • Hence P2 = A P1 + Ε is a more appropriate form of the equation, where Ε is a 2 x N matrix containing small values representing errors – often called as noise.
  • The least squares solution is also applicable in such a situation, with some very mild restrictions on E that are usually satisfied in practice.

59

60 of 90

Application of over-determined systems

  • What will happen if the points you marked out in any one of the images were coincident or collinear? Answer: If points are coincident or they lie on a line passing through the origin, system is undetermined. If points lie on any other line, the system is not undetermined.
  • If all points are non-collinear, what is the minimum number of points (N) so that the system of equations has a solution? Answer is N = 2 as these provide 4 knowns (2 X coordinates, 2 Y coordinates). We have only 4 unknowns here.

60

61 of 90

Orthogonal Matrices

  • An orthogonal matrix is a square matrix (over the real field) whose inverse is equal to its transpose:

  • For square matrices over the complex field, the transpose is replaced by the conjugate transpose. Such a matrix is called unitary.

61

62 of 90

Properties of orthogonal matrices

  • Have determinant +1 or -1 (rotation or reflection matrices respectively).
  • The magnitude of any row or column vector is 1.
  • The dot-product of any two different row-vectors or any two different column-vectors is 0.
  • Orthogonal transformations always preserve the magnitude of a vector, and the dot product (hence also angle between) two vectors.

  • The identity matrix is a simple example of an orthogonal matrix.

62

63 of 90

Applications of orthogonal matrices

  • Discrete signals (represented vectors with n elements) are often expressed as linear combinations of the columns of some orthogonal matrices.
  • Common examples of such orthogonal matrices are the Discrete Cosine Transform (DCT) matrix, the Discrete Fourier Transform (DFT) matrix and the Haar Wavelet matrix.

63

64 of 90

Applications of orthogonal matrices

  • The afore-mentioned representation is widely used in signal compression.
  • The DCT matrix is used for compressing audio (in MP3) and images (in JPEG).
  • Why? Because it is observed that only few of the coefficients in the vector Θ (for audio or images) have large magnitudes.
  • The rest have magnitudes equal to or close to 0 – and need not be stored.
  • This is a property of audio signals and natural images – and NOT general vectors.

64

65 of 90

Applications of orthogonal matrices

  • The audio signals or images reconstructed from just these large coefficients (and setting all the rest to 0) are often indistinguishable from the original signals!
  • So you need to store fewer coefficients!
  • There is of course an error incurred in this process – but that error is not (or barely) noticeable to the eye/ear.

65

66 of 90

Ajit Rajwade

66

512 x 512 Barbara image

Image reconstructed from top 80,000 largest DCT coefficients

67 of 90

67

DCT coefficients of images/audio signals decay very rapidly. The top plot shows the magnitudes of the DCT coefficients in descending order of absolute value. The plot below shows the magnitudes of all except the 20 largest coefficients (in terms of absolute value).

68 of 90

How does the DCT matrix look like?

68

It can be shown that Q is orthonormal. Q-1 is called as the DCT matrix. Q provides a mapping from the signal to its DCT coefficients, whereas its inverse provides the mapping from the DCt coefficients to the signal values.

Here is an 8 x 8 DCT matrix

N x N DCT matrix

69 of 90

Discrete Cosine Transform (DCT) in 1D

69

In this slide and the ones that follow, F is the vector of DCT coefficients (i.e. 𝛉 in slide 63), f is the original signal (i.e. x in previous slides), n is the time index (for the original signal) and u is the frequency index.

70 of 90

Discrete Cosine Transform (DCT) in 1D

70

n

u

u

n

71 of 90

DCT

  • Expresses a signal as a linear combination of cosine bases (as opposed to the complex exponentials as in the Fourier transform).
  • The coefficients of this linear combination are called DCT coefficients.
  • Is real-valued unlike the Fourier transform!
  • Discovered by Ahmed, Natarajan and Rao (1974)

71

72 of 90

72

u

n

  • DCT basis matrix is orthonormal. The dot product of any row (or column) with itself is 1. The dot product of any two different rows (or two different columns) is 0. The inverse is equal to the transpose.

  • Being orthonormal, it preserves the squared norm, i.e.

  • DCT is NOT the real part of the Fourier!

  • DCT basis matrix is NOT symmetric.

  • Columns of the DCT matrix are called the DCT basis vectors.

73 of 90

Digression: matrix view of a discrete orthonormal transform (Fourier transform used as example here)

  • Remember:

  • In matrix form, we write:

73

Fourier matrix: in any row, the value of x is fixed, the value of u ranges from 0 to M-1

74 of 90

74

75 of 90

Vector Norms

  • A norm is a function that assigns a real-valued length to a vector x.
  • It is denoted by ‖x‖.
  • To conform to a notion of length, the norm must satisfy the following for vectors x and y:

75

triangle inequality

76 of 90

Vector Norms

  • The most common norm is the 2-norm (or the L2 norm) or the usual notion of Euclidean length.
  • But there are many other types of norms and this additional flexibility is useful in many applications.
  • The most common class of norms is the p-norm (or Lp-norm) defined as follows:

76

77 of 90

Weighted Norms

  • In a weighted p-norm, different coordinates of a vector can be given different weights.
  • This is defined as follows:

77

In this definition, W is a diagonal matrix. But a weighted norm may be defined with a non-diagonal matrix as well.

78 of 90

Matrix Norms

  • If a m x n matrix A is interpreted as a vector in mn-dimensional space, we can use an Lp-norm from the previous slides.
  • For example, there is the Frobenius norm:

78

79 of 90

79

80 of 90

Matrix Norms

  • However there are certain norms which take into account that the matrix is a 2D entity and not a reshaped vector.
  • They are called induced norms.
  • The domain of a m x n matrix A is the set of n-dimensional vectors.
  • The range is the set of m-dimensional vectors of the form Ax where x is in the domain. The range is a subspace of the m-dimensional space – it may not include all of Rm (why?).
  • Thus Range(A) = {Ax| x lies in domain(A)}.

80

81 of 90

Matrix Norms

  • Definition of an induced norms:

Given norms ∥.∥(m) and ∥.∥(n) on the range and domain of A, the induced norm ∥A(m,n)is the smallest number C such that the following is true for all vectors x in R^n:

81

This is the maximum factor by which a vector can stretch under the action of A. For the purposes of this course, sup and max are (roughly) equivalent.

82 of 90

82

83 of 90

83

84 of 90

Matrix Norm Inequalities

  • All matrix norms (induced or not) satisfy the following:

  • Induced norms additionally satisfy the following conditions:

84

The last inequality follows from the very definition of ||A|| – why?

85 of 90

Why study matrix norms?

  • Matrix norms are very useful in high-dimensional statistics.
  • We will see a more immediate application called the condition number.
  • The condition number of a non-singular square matrix A with respect to a given norm is defined as cond(A) = ∥A∥∥A-1∥ (most commonly, the 2-norm is used).
  • If A is singular, then cond(A) is infinity.
  • If A equals the identity matrix, then cond(A) = 1.
  • The larger the condition number, the more “singular” is the matrix and the closer it is to 1, the “better behaved” it is.
  • cond(A) measures the ratio of the maximum relative stretching to the maximum relative shrinking by A to any nonzero vector x (see examples on next slide).

85

86 of 90

86

source: chatgpt (there is a mistake on this slide: what is it?)

87 of 90

Why study matrix norms?

  • In general, cond(A) ≥ 1 (why?).
  • For a diagonal matrix A, cond(A) = max |Aii|/ min(|Aii|).
  • If c is a scalar cond(cA) = cond(A). Why?
  • So far, we have seen that the condition number is an indicator of singularity or near-singularity of a matrix A.
  • But it has other important applications: consider the solution x to a linear system Ax = b where A is invertible and has size n x n, and x,b have size n x 1.
  • cond(A) provides a bound for the error in the computed solution if b is perturbed by some 𝛥b.
  • Let the solution to Ax = b+𝛥b be x’. Let x’ = x + 𝛥x.
  • Then Ax’=b+𝛥b A(x + 𝛥x)=b+𝛥b A𝛥x=𝛥b since we know Ax = b.
  • Then 𝛥x=A-1𝛥b → ‖𝛥x‖=‖A-1𝛥b‖ ≤ ‖A-1‖‖𝛥b
  • We also similarly have ‖b‖=‖Ax‖≤ ‖A‖‖x‖ → x‖ ≥ ‖b‖/‖A.
  • Combining these two inequalities we have:

𝛥x‖/‖x‖≤‖A-1‖‖𝛥b‖ ‖A‖/‖b‖=cond(A)‖𝛥b‖/‖b

87

88 of 90

Why study matrix norms?

  • The previous inequality tells us that for a matrix A with smaller cond(A), the effect of a small perturbation of b on the solution is small.
  • The condition number is like an amplification factor that amplifies the relative error in b.
  • A similar result holds if the entries in the matrix A are perturbed.
  • Let Ax = b and let (A+E)x’ = b.
  • Then 𝛥x=x’-x=A-1Ax’-A-1b = A-1(Ax’-b) = -A-1Ex’.
  • Combining with x‖ ≥ ‖b‖/‖A‖, we get

𝛥x‖/‖x’‖≤‖A-1‖‖E‖ = cond(A) ‖E‖/‖A

  • This result is not so useful due to x’ instead of x in the denominator.
  • We will derive a stronger result on the next slide.

88

89 of 90

89

Assume

If the condition number is high, then only a smaller perturbation is accommodated

90 of 90

References

  • Scientific Computing, Michael Heath
  • Numerical Linear Algebra, Treftehen and Bau
  • http://en.wikipedia.org/wiki/Matrix_(mathematics)

90