Difference between revisions of "Linear Algebra"

From Sutherland_wiki
Jump to: navigation, search
m (Matrix-Matrix Product)
m (Matrix-Matrix Product)
 
(19 intermediate revisions by 2 users not shown)
Line 2: Line 2:
 
Let's first make some definitions:
 
Let's first make some definitions:
  
; Vector : A one-dimensional collection of numbers. Examples:
+
; Vector : A one-dimensional collection of numbers. Examples: <math>a=\left( \begin{smallmatrix} 1 & 2 & 3 & 4 \end{smallmatrix} \right), \quad b=\left( \begin{smallmatrix} 1 \\ 2 \\ 3 \\ 4 \end{smallmatrix} \right)</math>
: <math>a=\left( \begin{array}{cccc} 1 & 2 & 3 & 4 \end{array} \right), \quad b=\left( \begin{array}{c} 1 \\ 2 \\ 3 \\ 4 \end{array} \right)</math>
 
  
; Matrix : A two-dimensional collection of numbers. Examples:
+
; Matrix : A two-dimensional collection of numbers. Examples: <math> A = \left[ \begin{smallmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{smallmatrix} \right], \quad B = \left[ \begin{smallmatrix} 9 & 6 \\ 2 & 3 \end{smallmatrix}\right]</math>
: <math> A = \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \end{array} \right], \quad B = \left[ \begin{array}{cc} 9 & 6 \\ 2 & 3 \end{array}\right]</math>
 
  
 
; Array : An ''n''-dimensional collection of numbers.  A ''vector'' is a 1-dimensional array while a ''matrix'' is a 2-dimensional array.
 
; Array : An ''n''-dimensional collection of numbers.  A ''vector'' is a 1-dimensional array while a ''matrix'' is a 2-dimensional array.
  
; Transpose : An operation that involves interchanging the rows and columns of an array.  It is indicated by a <math>^\mathsf{T}</math> superscript. For example:
+
; Transpose : An operation that involves interchanging the rows and columns of an array.  It is indicated by a <math>^\mathsf{T}</math> superscript. For example: <math>a=\left( \begin{smallmatrix} 1 & 2 & 3 & 4 \end{smallmatrix} \right) \; \Rightarrow \; a^\mathsf{T} = \left( \begin{smallmatrix} 1 \\ 2 \\ 3 \\ 4 \end{smallmatrix} \right) </math> for a vector and <math> A = \left[ \begin{smallmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{smallmatrix} \right] \; \Rightarrow \; A^\mathsf{T} = \left[ \begin{smallmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{smallmatrix} \right] </math> for a matrix
: <math>a=\left( \begin{array}{cccc} 1 & 2 & 3 & 4 \end{array} \right) \; \Rightarrow \; a^\mathsf{T} = \left( \begin{array}{c} 1 \\ 2 \\ 3 \\ 4 \end{array} \right) </math>
 
: <math> A = \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 4 & 5 & 6 \end{array} \right] \; \Rightarrow \;
 
  A^\mathsf{T} = \left[ \begin{array}{cc} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{array} \right]  
 
</math>
 
  
; Vector Magnitude, <math>\left|\vec{a}\right|</math> : A measure of the length of a vector, <math>\left|\vec{a}\right| = \sqrt{ \sum_{i=1}^{n} a_i^2 }</math>
+
; Vector Magnitude, <math>\left|\vec{a}\right|</math>
: Example: <math>\vec{a}=\left( \begin{array}{ccc} 2 & 8 & 3 \end{array} \right) \; \Rightarrow \; \left|\vec{a}\right| = \sqrt{ 2^2 + 8^2 + 3^2 }</math>
+
: A measure of the length of a vector, <math>\left|\vec{a}\right| = \sqrt{ \sum_{i=1}^{n} a_i^2 }</math>
: Example: <math>\vec{b}=\left( \begin{array}{cccc} b_1 & b_2 & b_3 & b_4\end{array} \right) \; \Rightarrow \; \left|\vec{b}\right| = \sqrt{ b_1^2 + b_2^2 + b_3^2 +b_4^2 } </math>
+
: <math>\vec{a}=\left( \begin{smallmatrix} 2 & 8 & 3 \end{smallmatrix} \right) \; \Rightarrow \; \left|\vec{a}\right| = \sqrt{ 2^2 + 8^2 + 3^2 }</math>
 +
: <math>\vec{b}=\left( \begin{smallmatrix} b_1 & b_2 & b_3 & b_4\end{smallmatrix} \right) \; \Rightarrow \; \left|\vec{b}\right| = \sqrt{ b_1^2 + b_2^2 + b_3^2 +b_4^2 } </math>
  
 +
; Vector Norm, <math>\left\Vert\vec{a}\right\Vert</math>
 +
: A measure of the size of a vector
 +
{|border="1" cellpadding="5" cellspacing="0" align="center"
 +
|+ Definition of various vector norms
 +
! Name !! Definition
 +
|-
 +
| L<sub>1</sub> norm
 +
|| <math>\left\Vert \vec{a} \right\Vert_1 = \sum_i a_i</math>
 +
|-
 +
| L<sub>2</sub> norm
 +
|| <math>\left\Vert \vec{a} \right\Vert_2 = \sqrt{\sum_i a_i^2}</math>
 +
|-
 +
| L<sub>&infin;</sub> norm
 +
|| <math>\left\Vert \vec{a} \right\Vert_2 = \max{ a_i}</math>
 +
|}
  
 
== Matrix & Vector Algebra ==
 
== Matrix & Vector Algebra ==
Line 29: Line 39:
 
* [[#Matrix-Matrix Product|Matrix-matrix product]]
 
* [[#Matrix-Matrix Product|Matrix-matrix product]]
 
These are each discussed in the following sub-sections.
 
These are each discussed in the following sub-sections.
 +
  
 
=== Vector Dot Product ===
 
=== Vector Dot Product ===
Line 41: Line 52:
 
Occasionally we know the magnitude of the two vectors and the angle between them.  In this case, we can calculate the dot product as  
 
Occasionally we know the magnitude of the two vectors and the angle between them.  In this case, we can calculate the dot product as  
 
<center><math>
 
<center><math>
  |c| = |a|\,|b|\,\cos(\theta),
+
  c = |\vec{a}|\,|\vec{b}|\,\cos(\theta),
 
</math></center>
 
</math></center>
 
where <var>&theta;</var> is the angle between the vectors '''a''' and '''b'''.
 
where <var>&theta;</var> is the angle between the vectors '''a''' and '''b'''.
 +
 +
Note that if '''a''' and '''b''' are both column vectors, then <math>\vec{a}\cdot \vec{b} = a^\mathsf{T}b</math>.
  
 
For more information on the dot product, see [http://en.wikipedia.org/wiki/Dot_product wikipedia's article].
 
For more information on the dot product, see [http://en.wikipedia.org/wiki/Dot_product wikipedia's article].
 +
 +
==== Dot Product in MATLAB ====
 +
There are several ways to calculate the dot product of two vectors in MATLAB:
 +
<source lang="matlab">
 +
  c = dot(a,b);
 +
  c = a' * b;    % if a and b are both column vectors.
 +
</source>
  
 
==== Thought examples ====
 
==== Thought examples ====
Line 51: Line 71:
  
 
* At noon, when the sun is directly overhead, a rocket is launched directly toward the sun at 1000 mph.  How fast is its shadow moving?
 
* At noon, when the sun is directly overhead, a rocket is launched directly toward the sun at 1000 mph.  How fast is its shadow moving?
: Define the rocket's velocity as <math>\vec{R}</math>.
+
*# Define the rocket's velocity as <math>\vec{R}</math>.
: Define the unit normal on the ground (i.e. the direction of the rocket's shadow on the ground) as <math>\vec{s}</math>, where <math>|\vec{s}| = 1</math>.
+
*# Define the unit normal on the ground (i.e. the direction of the rocket's shadow on the ground) as <math>\vec{s}</math>, where <math>|\vec{s}| = 1</math>.
: Intuition tells us that if the it is moving directly toward the sun, then its shadow does not appear to move at all. : The dot product <math>\vec{R}\cdot\vec{s} = |\vec{R}| \, |\vec{s}| \cos(\theta) = 0</math> since cos(90&deg;)=0.
+
*# The dot product between <math>\vec{R}</math> and <math>\vec{s}</math> gives us the component of the velocity in the direction of <math>\vec{s}</math>.  Thus, the speed of the shadow along the ground is given by
 +
*#: <math>v = \vec{R}\cdot\vec{s} = |\vec{R}|\,|\vec{s}| \cos(\theta)</math>
 +
*# Intuition tells us that if the rocket is moving directly toward the sun, then its shadow does not appear to move at all. In this case, the angle between <math>\vec{R}</math> and <math>\vec{s}</math> is 90&deg; and the dot product <math>\vec{R}\cdot\vec{s} = |\vec{R}| \, |\vec{s}| \cos(\theta) = 0</math> since <math>\cos(90^\circ)=0</math>.
  
 
* Consider the same rocket going parallel to the earth's surface.  How fast is its shadow moving?
 
* Consider the same rocket going parallel to the earth's surface.  How fast is its shadow moving?
: If the rocket is going parallel to the earth's surface, our intuition tells us that is shadow is moving at the same speed.  This is confirmed by the mathematics, since the angle between the rocket's path and the ground is 0. Therefore, <math>\vec{R}\cdot\vec{s} = |\vec{R}| \, |\vec{s}| \cos(\theta) = 0</math> since cos(0&deg;)=1.
+
: If the rocket is going parallel to the earth's surface, our intuition tells us that is shadow is moving at the same speed.  This is confirmed by the mathematics, since the angle between the rocket's path and the ground is 0. Therefore, <math>\vec{R}\cdot\vec{s} = |\vec{R}| \, |\vec{s}| \cos(\theta) = |\vec{R}| \, |\vec{s}| = |\vec{R}|</math> since <math>\cos(0^\circ)=1</math> and <math>|\vec{s}|=1</math>.
  
 
* What if the rocket were going at a 45&deg; angle?
 
* What if the rocket were going at a 45&deg; angle?
Line 87: Line 109:
  
 
For more information, see [http://en.wikipedia.org/wiki/Cross_product wikipedia's article on the cross product].
 
For more information, see [http://en.wikipedia.org/wiki/Cross_product wikipedia's article on the cross product].
 +
 +
==== Vector Cross Product in MATLAB ====
 +
<source lang="matlab">
 +
  c = cross(a,b);  % c is a vector!
 +
</source>
 +
  
 
=== Matrix-Vector Product ===
 
=== Matrix-Vector Product ===
  
The matrix-vector product is frequently encountered in solving [[LinearAlgebra#Linear_Systems_of_Equations|linear systems of equations]], although there are many other uses for it.
+
The matrix-vector product is frequently encountered in solving [[Linear_Algebra#Linear_Systems_of_Equations|linear systems of equations]], although there are many other uses for it.
  
 
Consider a matrix '''A''' with ''m'' rows and ''n'' columns and vector '''b''' with ''n'' rows:
 
Consider a matrix '''A''' with ''m'' rows and ''n'' columns and vector '''b''' with ''n'' rows:
Line 149: Line 177:
 
  c=\left(\begin{array}{c} 34 \\ 64 \end{array}\right).
 
  c=\left(\begin{array}{c} 34 \\ 64 \end{array}\right).
 
</math></center>
 
</math></center>
 +
 +
==== Geometric Interpretation of Matrix-Vector Product ====
 +
Here is a great discussion of how [http://www.mathworks.com/moler/exm/chapters/matrices.pdf a matrix-vector product can be interpreted geometrically].
 +
 +
==== Matrix-Vector product in MATLAB ====
 +
<source lang="matlab">
 +
  c = A*b;  % A is a matrix, b and c are vectors.
 +
</source>
 +
  
 
=== Matrix-Matrix Product ===
 
=== Matrix-Matrix Product ===
Line 157: Line 194:
 
# The number of rows in '''C''' is equal to the number of rows in '''A'''.
 
# The number of rows in '''C''' is equal to the number of rows in '''A'''.
 
We can state this differently.  Suppose that '''A''' has <var>m</var> rows and <var>p</var> columns, and that '''B''' has <var>p</var> rows and <var>n</var> columns.  This satisfies rule 1 above.  Then, according to rules 2 and 3, '''C''' has <var>m</var> rows and <var>n</var> columns.
 
We can state this differently.  Suppose that '''A''' has <var>m</var> rows and <var>p</var> columns, and that '''B''' has <var>p</var> rows and <var>n</var> columns.  This satisfies rule 1 above.  Then, according to rules 2 and 3, '''C''' has <var>m</var> rows and <var>n</var> columns.
 +
 +
{|border="1" cellpadding="5" cellspacing="0" align="center"
 +
|+ Examples of valid and invalid matrices for multiplication
 +
! A !! B !! Valid? !! Comments
 +
|-
 +
| <math>\left[\begin{smallmatrix} 1 & 2 \\ 3 & 4 \end{smallmatrix}\right]</math>
 +
||<math>\left[\begin{smallmatrix} 5 & 6 \\ 7 & 8 \end{smallmatrix}\right]</math>
 +
| style="color:green" | Yes
 +
|| Number of columns in A = 2, number or rows in B = 2.
 +
|-
 +
| <math>\left[\begin{smallmatrix} 1 & 2 \end{smallmatrix}\right]</math>
 +
||<math>\left[\begin{smallmatrix} 5 \\ 8 \end{smallmatrix}\right]</math>
 +
| style="color:green" | Yes
 +
||
 +
|-
 +
| <math>\left[\begin{smallmatrix} 1 \\ 2 \end{smallmatrix}\right]</math>
 +
||<math>\left[\begin{smallmatrix} 5 & 8 \end{smallmatrix}\right]</math>
 +
| style="color:green" | Yes
 +
||
 +
|-
 +
| <math>\left[\begin{smallmatrix} 1 & 2 \end{smallmatrix}\right]</math>
 +
||<math>\left[\begin{smallmatrix} 5 & 8 \end{smallmatrix}\right]</math>
 +
| style="color:red" | No
 +
|| Number of columns in A=2, number of rows in B = 1.  Therefore these cannot be multiplied.
 +
|}
  
 
The general formula for matrix multiplication is
 
The general formula for matrix multiplication is
Line 176: Line 238:
 
</math></center>
 
</math></center>
  
We can form '''C'''='''A B''' since rule 1 is satisfied.  From rules 2 and 3, we conclude that the resulting matrix, '''C''' has 2 rows and 3 columns.
+
We can form '''C'''='''A B''' since rule 1 is satisfied.  From the size of '''A''' and '''B''' we define <var>m=2</var> and <var>n=3</var> and from rules 2 and 3, we conclude that the resulting matrix, '''C''' has 2 rows and 3 columns.
  
 
We are now ready to build each element of C.
 
We are now ready to build each element of C.
 
+
{|border="1" cellpadding="5" cellspacing="0" align="center" style="text-align:center"
; <var>i=1</var> and <var>j=1</var>
+
! Row !! Column !! Result
: <math>
+
|-
 +
| 1 || 1
 +
|| <math>
 
   C_{1,1} = \sum_{k=1}^2 A_{1,k} B_{k,1}
 
   C_{1,1} = \sum_{k=1}^2 A_{1,k} B_{k,1}
 
           = A_{1,1} B_{1,1} + A_{1,2} B_{2,1}
 
           = A_{1,1} B_{1,1} + A_{1,2} B_{2,1}
Line 187: Line 251:
 
           = 14.
 
           = 14.
 
  </math>  
 
  </math>  
 
+
|-
; <var>i=2</var> and <var>j=1</var>
+
| 2 || 1
: <math>
+
|| <math>
 
   C_{2,1} = \sum_{k=1}^2 A_{2,k} B_{k,1}
 
   C_{2,1} = \sum_{k=1}^2 A_{2,k} B_{k,1}
 
           = A_{2,1} B_{1,1} + A_{2,2} B_{2,1}
 
           = A_{2,1} B_{1,1} + A_{2,2} B_{2,1}
Line 195: Line 259:
 
           = 36.
 
           = 36.
 
   </math>
 
   </math>
; <var>i=1</var> and <var>j=2</var>
+
|-
: <math>
+
| 1 || 2
 +
||<math>
 
   C_{1,2} = \sum_{k=1}^2 A_{1,k} B_{k,2}  
 
   C_{1,2} = \sum_{k=1}^2 A_{1,k} B_{k,2}  
 
           = A_{1,1} B_{1,2} + A_{1,2} B_{2,2}  
 
           = A_{1,1} B_{1,2} + A_{1,2} B_{2,2}  
Line 202: Line 267:
 
           = 10.
 
           = 10.
 
   </math>
 
   </math>
; <var>i=2</var> and <var>j=2</var>
+
|-
: <math>
+
| 2 || 2
 +
|| <math>
 
   C_{2,2} = \sum_{k=1}^2 A_{2,k} B_{k,2}  
 
   C_{2,2} = \sum_{k=1}^2 A_{2,k} B_{k,2}  
 
           = A_{2,1} B_{1,2} + A_{2,2} B_{2,2}  
 
           = A_{2,1} B_{1,2} + A_{2,2} B_{2,2}  
Line 209: Line 275:
 
           = 27.
 
           = 27.
 
   </math>
 
   </math>
; <var>i=1</var> and <var>j=3</var>
+
|-
: <math>
+
| 1 || 3
 +
|| <math>
 
   C_{1,3} = \sum_{k=1}^2 A_{1,k} B_{k,3}  
 
   C_{1,3} = \sum_{k=1}^2 A_{1,k} B_{k,3}  
 
           = A_{1,1} B_{1,3} + A_{1,2} B_{2,3}  
 
           = A_{1,1} B_{1,3} + A_{1,2} B_{2,3}  
Line 216: Line 283:
 
           = 12.
 
           = 12.
 
   </math>
 
   </math>
; <var>i=2</var> and <var>j=3</var>
+
|-
: <math>
+
| 2 || 3
 +
| <math>
 
   C_{2,3} = \sum_{k=1}^2 A_{2,k} B_{k,3}  
 
   C_{2,3} = \sum_{k=1}^2 A_{2,k} B_{k,3}  
 
           = A_{2,1} B_{1,3} + A_{2,2} B_{2,3}  
 
           = A_{2,1} B_{1,3} + A_{2,2} B_{2,3}  
Line 223: Line 291:
 
           = 35.
 
           = 35.
 
   </math>
 
   </math>
 +
|}
  
 
Now, putting all of these together, we have
 
Now, putting all of these together, we have
Line 231: Line 300:
 
\end{array}\right].
 
\end{array}\right].
 
</math></center>
 
</math></center>
 +
 +
==== Matrix-Matrix product in MATLAB ====
 +
<source lang="matlab">
 +
  C = A*B;  % A, B, and C are all matrices.
 +
</source>
 +
The MATLAB code to implement the above example is
 +
<source lang="matlab">
 +
  A = [ 1 2; 3 5; ];
 +
  B = [ 2 4 10; 6 3 1; ];
 +
  C = A*B;
 +
</source>
  
 
== Linear Systems of Equations ==
 
== Linear Systems of Equations ==
{{stub|section}}
+
Any set of linear systems may be written in the form <math>[A](x)=(b)</math>, where <math>[A]</math> is a matrix and <math>(x)</math> and <math>(b)</math> are vectors.  In general, for a system of equations with <var>n</var> unknowns <var>x<sub>1</sub> ... x<sub>n</sub></var>, can be written as
 +
<center><math>
 +
\begin{align}
 +
  a_{1,1} x_1 + a_{1,2} x_2 + a_{1,3} x_3 + \cdots + a_{1,n} x_n &= b_1, \\
 +
  a_{2,1} x_1 + a_{2,2} x_2 + a_{2,3} x_3 + \cdots + a_{2,n} x_n &= b_2, \\
 +
  &\vdots \\
 +
  a_{n,1} x_1 + a_{n,2} x_2 + a_{n,3} x_3 + \cdots + a_{n,n} x_n &= b_n, \\
 +
\end{align}
 +
</math></center>
 +
where <math>a_{i,j}</math> and <math>b_i</math> are numbers and <math>x_i</math> is the set of unknown variables.  These equations can now be written in matrix format:
 +
<center><math>
 +
\left[ \begin{array}{ccccc}
 +
  a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,n} \\
 +
  a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,n} \\
 +
  \vdots  & \vdots  & \ddots  & \vdots & \vdots \\
 +
  a_{n,1} & a_{n,2} & a_{n,3} & \cdots & a_{n,n}
 +
\end{array} \right]
 +
\left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right)
 +
=
 +
\left( \begin{array}{c} b_1 \\ b_2 \\ \vdots \\ b_n \end{array} \right).
 +
</math></center>
 +
 
 +
=== Example ===
 +
Consider the following equations:
 +
<center><math>
 +
\begin{cases}
 +
  x + y -5z = 9 \\
 +
  6-z = 2y \\
 +
  y-3x = 10z \\
 +
\end{cases}
 +
</math></center>
 +
To put these in matrix format, we first choose the solution vector.  Let's choose to order this as <var>x, y, z</var>.  Now we rewrite the original system so that we have a coefficient multiplying each of these variables, and then we put the "constant" on the right-hand-side:
 +
<center><math>
 +
\begin{align}
 +
    1x + 1y -5z  &=  9 \\
 +
    0x - 2y -1z  &= -6 \\
 +
  -3x + 1y -10z &=  0
 +
\end{align}
 +
</math></center>
 +
We may then identify the '''A''' matrix and '''b''' vectors and write
 +
<center><math>
 +
\left[\begin{array}{ccc}
 +
  1 &  1 & -5  \\
 +
  0 & -2 & -1  \\
 +
-3 &  1 & -10 \\
 +
\end{array} \right]
 +
\left( \begin{array}{c}
 +
  x \\ y \\ z
 +
\end{array}\right)
 +
=
 +
\left( \begin{array}{c}
 +
  9 \\ -6 \\ 0
 +
\end{array}\right)
 +
</math></center>
 +
Note that if we multiply the above matrix and vector, we recover the original set of equations (try it to be sure that you can do it).
 +
 
 +
Note that if we chose a different variable ordering, then it simply rearranges the columns in the matrix.  For example, if we order the variables <var>z, x, y</var> then we have
 +
<center><math>
 +
\begin{align}
 +
    -5z  + 1x + 1y &=  9 \\
 +
    -1z  + 0x - 2y &= -6 \\
 +
    -10z - 3x + 1y &=  0
 +
\end{align}
 +
\quad \Rightarrow \quad
 +
\left[\begin{array}{ccc}
 +
  -5  &  1 &  1 \\
 +
  -1  &  0 & -2 \\
 +
  -10 & -3 &  1 \\
 +
\end{array} \right]
 +
\left( \begin{array}{c}
 +
  z \\ x \\ y
 +
\end{array}\right)
 +
=
 +
\left( \begin{array}{c}
 +
  9 \\ -6 \\ 0
 +
\end{array}\right)
 +
 +
</math></center>
 +
 
 +
 
 +
This page discusses a few methods to solve linear systems of
 +
equations, <math>Ax=b</math>. There are two categories of solution
 +
methods for linear systems of equations:
 +
[[#Direct_Solvers|direct solvers]] and
 +
[[#Iterative_Solvers|iterative solvers]]. These will be discussed
 +
briefly below.
 +
 
 +
 
 +
== Direct Solvers ==
 +
 
 +
Direct solvers solve the linear system to within the numeric precision
 +
of the computer. There are many different algorithms for direct
 +
solution of linear systems, but we will only mention two:
 +
[[#Gaussian_Elimination|Gaussian elimination]] and
 +
[[#The_Thomas_Algorithm (Tridiagonal Systems)|the Thomas Algorithm]].
 +
 
 +
 
 +
=== Gaussian Elimination ===
 +
 
 +
Gaussian elimination is a methododical approach to solving linear
 +
systems of equations. It is essentially what you would do if you were
 +
to solve a system of equations "by hand" but is systematic so that it
 +
can be implemented on the computer.
 +
For information on the algorithm, see this
 +
[http://en.wikipedia.org/wiki/Gaussian_elimination wikipedia article].
 +
 
 +
=== The Thomas Algorithm (Tridiagonal Systems) ===
 +
 
 +
The Thomas Algorithm is used for solving tridiagonal systems of
 +
equations, <math>Ax=b</math>. A tridiagonal matrix is one where the
 +
only non-zero entries in the matrix are on the main diagonal and the
 +
lower and upper diagonals. In general a tridiagonal matrix can be
 +
written as
 +
:<math>A = \left[\begin{array}{ccccc}
 +
  d_1    & u_1    & 0      & \cdots    & 0      \\
 +
  \ell_1 & d_2    & u_2    & 0          & \cdots \\
 +
        & \ddots & \ddots & \ddots    &        \\
 +
  0    & 0      & \cdots & \ell_{n-1} & d_n
 +
\end{array}\right]
 +
</math>
 +
 
 +
Here <math>d_i</math> represents the diagonal of the matrix,
 +
<math>\ell_i</math> represents the lower-diagonal, and
 +
<math>u_i</math> represents the upper-diagonal.
 +
 
 +
In terms of this nomenclature, the Thomas Algorithm can be written as
 +
<code>
 +
''forward pass...''
 +
'''for''' i = 2 ... n
 +
  d(i) = d(i) - ( u(i-1)*ell(i-1) ) / d(i-1)
 +
  b(i) = b(i) - ( b(i-1)*ell(i-1) ) / d(i-1)
 +
'''end'''
 +
x(n) = b(n) / d(n)
 +
''backward pass...''
 +
'''for''' i = n-1 ... 1
 +
  x(i) = ( b(i) - u(i) * x(i+1) ) / d(i)
 +
'''end'''
 +
</code>
 +
 
 +
== Iterative Solvers ==
 +
 
 +
Iterative solvers are fundamentally different than
 +
[[#Direct Solvers|direct solvers]]. They require an initial guess and
 +
then refine the solution until a tolerance is reached.
 +
 
 +
 
 +
=== Jacobi Method ===
 +
 
 +
The Jacobi method for solving linear systems is the simplest iterative
 +
method. In this case, each iteration of the solution involves cycling
 +
through the equations and solving the <var>i<sup>th</sup></var>
 +
equation for the <var>x<sub>i</sub></var>, using values for all other
 +
<var>x<sub>j</sub></var> from the previous iteration.
 +
 
 +
The algorithm for the Jacobi method is:
 +
:<code>
 +
xnew = x;
 +
err = tol+1;
 +
'''while''' err > tol
 +
  '''for''' i=1:1:n
 +
      ''update the i<sup>th</sup> value of x.''
 +
      tmp = 0
 +
      '''for''' j=1:1:i-1
 +
        tmp = tmp + A(i,j) * x(j)
 +
      '''end'''
 +
      '''for''' j=i+1:1:n
 +
        tmp = tmp + A(i,j) * x(j)
 +
      '''end'''
 +
      xnew(i) = x(i) + ( b(j) - tmp ) / A(i,i)
 +
  '''end'''
 +
  x = xnew
 +
  ''check for convergence - set'' '''err'''
 +
'''end'''
 +
</code>
 +
 
 +
<span id="linsys_example">  <!-- anchor for linking -->
 +
For example, assume that we have the following system of equations:
 +
:<math>
 +
\left[ \begin{array}{rrrr}
 +
-3 &  1 &  0 &  0  \\
 +
  1 & -2 &  1 &  0  \\
 +
  0 &  1 & -2 &  1  \\
 +
  0 &  0 &  1 & -3
 +
\end{array} \right]
 +
\left( \begin{array}{c}
 +
x_1 \\ x_2 \\ x_3 \\ x_4
 +
\end{array} \right)
 +
=
 +
\left( \begin{array}{r}
 +
-4 \\ 2 \\ 9 \\ 2
 +
\end{array} \right)
 +
</math>
 +
 
 +
If we set our initial guess as <math>x=\left[ \begin{array}{cccc} 1 & 1 &
 +
1 & 1 \end{array} \right]^T</math> then for the first iteration we
 +
proceed as follows:
 +
 
 +
* Solve the first equation for <var>x<sub>1</sub></var> to find <math> x_1 = \frac{-4 -x_2}{-3} = 5/3.</math>
 +
 
 +
* Solve the second equation for <var>x<sub>2</sub></var> to find <math> x_2 = \frac{2 - x_1-x_3}{-2} = 0.</math>
 +
 
 +
We repeat this for all of the equations. At the end of the first
 +
iteration, we have
 +
<math>x=\left(\begin{array}{cccc} 5/3, & 0, & -7/2, & -1/3 \end{array}\right).</math>
 +
 
 +
We are now ready to start with the second iteration, where we use the
 +
newly computed values for <math>x_i</math>. The
 +
[[#convergencetable|table below]] shows the results for the first few
 +
iterations.
 +
 
 +
The last entry in the table is the analytic solution to this system of
 +
equations. Note that after 40 iterations we have converged to the
 +
exact answer to within an absolute error <var>x-x<sub>exact</sub> =
 +
10<sup>-3</sup></var>.
 +
 
 +
See [[media:jacobi.m|here]] for a Matlab implementation of the Jacobi
 +
method. More information can be found
 +
[http://en.wikipedia.org/wiki/Jacobi_method here].
 +
 
 +
 
 +
=== Gauss-Seidel Method ===
 +
 
 +
The Gauss-Seidel method is a simple extension of the Jacobi method.
 +
The only difference is that rather than using "old" solution values,
 +
we use the most up-to-date solution values in our calculation.
 +
 
 +
The algorithm for the Gauss-Seidel method is:
 +
:<code>
 +
err = tol+1
 +
'''while''' err > tol
 +
  '''for''' i=1:1:n
 +
      ''update x<sub>i</sub> with latest information on all other values of x''
 +
      '''for''' j=1:1:n
 +
        tmp = tmp + A(i,j) * x(j)
 +
      '''end'''
 +
      x(i) = x(i) + ( b(i) - tmp ) / A(i,i)
 +
  '''end'''
 +
  ''check for convergence - set'' '''err'''
 +
'''end'''
 +
</code>
 +
 
 +
If we apply the Gauss-Seidel method to the
 +
[[#linsys_example|example above]], we proceed as follows for the first
 +
iteration:
 +
 
 +
* Solve the first equation for <var>x<sub>1</sub></var> to find <math>x_1 = \frac{-4 -x_2}{-3} = 5/3</math>
 +
 
 +
* Solve the second equation for <var>x<sub>2</sub></var> and use the most recent value of <var>x<sub>1</sub></var> to find <math>x_2 = \frac{2 - x_1-x_3)}{-2} = \frac{2-5/3-1}{-2} = 1/3</math>
 +
 
 +
We repeat this for all of the equations. At the end of the first
 +
iteration, we have
 +
<math>x=\left(\begin{array}{cccc} 5/3, & 1/3, & -23/6, & -35/18, \end{array}\right)</math>
 +
 
 +
The following table shows the convergence results for the Jacobi and
 +
Gauss-Seidel methods. The Gauss-Seidel method converges much faster
 +
than the Jacobi method.
 +
 
 +
<span id="convergencetable">  <!-- anchor for linking -->
 +
 
 +
:{| border="1" cellpadding="3" cellspacing="0" style="text=align:center"
 +
|+ Convergence for the Jacobi and Gauss-Seidel iterative linear solvers
 +
!
 +
! style="text-align:center" colspan="4" | Jacobi Method
 +
!
 +
! style="text-align:center" colspan="4" | Gauss-Seidel Method
 +
|-
 +
| style="text-align:center" | Iteration
 +
| <math>x_1</math> || <math>x_2</math> || <math>x_3</math> || <math>x_4</math>
 +
|
 +
| <math>x_1</math> || <math>x_2</math> || <math>x_3</math> || <math>x_4</math>
 +
|-
 +
! 1
 +
| 5/3  || 0    || -7/2  || -1/3
 +
|
 +
| 5/3  || 1/3  || -23/6 || -35/18
 +
|-
 +
! 2
 +
| 4/3    || -23/12 || -14/3 || -11/6
 +
|
 +
| 1.4444 || -2.1944 || -6.5694 || -2.8565
 +
|-
 +
! 3
 +
| 25/26 || -8/3  || -51/8 || -20/9
 +
|
 +
| 0.6019 || -3.9838 || -7.9201 || -3.3067
 +
|-                                           
 +
! 4
 +
| 0.4444 || -3.8403 || -6.9444 || -2.7917
 +
|
 +
| 0.0054 || -4.9574 || -8.6320 || -3.5440
 +
|-                                           
 +
! 5
 +
|  0.0532 || -4.2500 || -7.8160 || -2.9815
 +
|
 +
| -0.3191 || -5.4756 || -9.0098 || -3.6699
 +
|-
 +
! 10
 +
| -0.5170 || -5.7294 || -9.0648 || -3.6601
 +
|
 +
| -0.6719 || -6.0377 || -9.4194 || -3.8065
 +
|-                                             
 +
! 20
 +
| -0.6803 || -6.0484 || -9.4218 || -3.8061
 +
|
 +
| -0.6875 || -6.0625 || -9.4375 || -3.8125
 +
|-
 +
! &infin;
 +
| -11/16 = -0.6875 || -97/16 = -6.0625 || -151/16 = -9.4375 || -61/16 = -3.8125
 +
|}
 +
 
 +
More information on the Gauss-Seidel method can be found
 +
[http://en.wikipedia.org/wiki/Gauss-Seidel#Convergence here].
 +
 
 +
 
 +
=== Some Caveats ===
 +
 
 +
A [[#Direct Solvers|direct solver]] will always find a solution if the
 +
system is well-posed. However, iterative solvers are more fragile, and
 +
will occasionally not converge. This is particularly true for the
 +
Jacobi and Gauss-Seidel solvers discussed above. In fact, the matrices
 +
must have special properties for these solvers to converge.
 +
Specifically, they must be symmetric and positive-definite or
 +
diagonally dominant. What this means in practice is that in a given
 +
row of the matrix, the diagonal coefficient is larger than the sum of
 +
the off-diagonal coefficients,
 +
:<math>\left| a_{ii} \right| > \sum_{j\ne i} \left| a_{ij} \right|</math>
 +
 
 +
=== Other Iterative Methods ===
 +
 
 +
There are many other iterative solvers that are much more
 +
sophisticated and robust. Among them are the
 +
[http://en.wikipedia.org/wiki/Conjugate_gradient Conjugate-Gradient]
 +
and [http://en.wikipedia.org/wiki/Gmres GMRES] methods.
 +
 
 +
 
 +
=== Convergence ===
  
 +
Iterative solvers require us to decide when to stop iterating. In
 +
other words, we need to decide when the solution is ''converged''. We
 +
do this by comparing an error measure to a tolerance.
  
== Solving Linear Systems of Equations ==
+
There are a few common ways to determine convergence for a linear
{{stub|section}}
+
system:
=== Direct Solvers ===
+
* <math>\epsilon = \left\Vert b-Ax \right\Vert</math>
==== Gaussian Elimination ====
+
* <math>\epsilon = \left\Vert x^{k}-x^{k-1}\right\Vert</math>
==== The Thomas Algorithm (Tridiagonal Systems) ====
+
where the terms <math>\left\Vert a \right\Vert</math> indicate an
 +
appropriate [[Iteration_and_Convergence#Vector_Norms|vector norm]].
  
=== Iterative Solvers ===
+
See the discussion on
==== Jacobi Method ====
+
[[Iteration_and_Convergence|iteration and convergence]] for more
==== Gauss-Seidel Method ====
+
information on this topic.
==== Other Methods ====
 
* Conjugate-Gradient
 
* GMRES
 

Latest revision as of 12:23, 26 August 2016

Basics

Let's first make some definitions:

Vector 
A one-dimensional collection of numbers. Examples: a=\left( \begin{smallmatrix} 1 & 2 & 3 & 4 \end{smallmatrix} \right), \quad b=\left( \begin{smallmatrix} 1 \\ 2 \\ 3 \\ 4 \end{smallmatrix} \right)
Matrix 
A two-dimensional collection of numbers. Examples:  A = \left[ \begin{smallmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{smallmatrix} \right], \quad B = \left[ \begin{smallmatrix} 9 & 6 \\ 2 & 3 \end{smallmatrix}\right]
Array 
An n-dimensional collection of numbers. A vector is a 1-dimensional array while a matrix is a 2-dimensional array.
Transpose 
An operation that involves interchanging the rows and columns of an array. It is indicated by a ^\mathsf{T} superscript. For example: a=\left( \begin{smallmatrix} 1 & 2 & 3 & 4 \end{smallmatrix} \right) \; \Rightarrow \; a^\mathsf{T} = \left( \begin{smallmatrix} 1 \\ 2 \\ 3 \\ 4 \end{smallmatrix} \right) for a vector and  A = \left[ \begin{smallmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{smallmatrix} \right] \; \Rightarrow \; A^\mathsf{T} = \left[ \begin{smallmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \end{smallmatrix} \right] for a matrix
Vector Magnitude, \left|\vec{a}\right|
A measure of the length of a vector, \left|\vec{a}\right| = \sqrt{ \sum_{i=1}^{n} a_i^2 }
\vec{a}=\left( \begin{smallmatrix} 2 & 8 & 3 \end{smallmatrix} \right) \; \Rightarrow \; \left|\vec{a}\right| = \sqrt{ 2^2 + 8^2 + 3^2 }
\vec{b}=\left( \begin{smallmatrix} b_1 & b_2 & b_3 & b_4\end{smallmatrix} \right) \; \Rightarrow \; \left|\vec{b}\right| = \sqrt{ b_1^2 + b_2^2 + b_3^2 +b_4^2 }
Vector Norm, \left\Vert\vec{a}\right\Vert
A measure of the size of a vector
Definition of various vector norms
Name Definition
L1 norm \left\Vert \vec{a} \right\Vert_1 = \sum_i a_i
L2 norm \left\Vert \vec{a} \right\Vert_2 = \sqrt{\sum_i a_i^2}
L norm \left\Vert \vec{a} \right\Vert_2 = \max{ a_i}

Matrix & Vector Algebra

There are many common operations involving matrices and vectors including:

These are each discussed in the following sub-sections.


Vector Dot Product

The dot product of two vectors produces a scalar. Physically, the dot product of a and b represents the projection of a onto b.

Given two vectors a and b, their dot product is formed as


 c = \vec{a} \cdot \vec{b} = \sum_{i=1}^n a_i b_i,

where a_i and b_i are the components in vectors a and b respectively. This is most useful when we know the components of the vectors a and b.

Occasionally we know the magnitude of the two vectors and the angle between them. In this case, we can calculate the dot product as


 c = |\vec{a}|\,|\vec{b}|\,\cos(\theta),

where θ is the angle between the vectors a and b.

Note that if a and b are both column vectors, then \vec{a}\cdot \vec{b} = a^\mathsf{T}b.

For more information on the dot product, see wikipedia's article.

Dot Product in MATLAB

There are several ways to calculate the dot product of two vectors in MATLAB:

  c = dot(a,b);
  c = a' * b;    % if a and b are both column vectors.

Thought examples

Consider the cartoon shown to the right.
DotProdExample.png
  • At noon, when the sun is directly overhead, a rocket is launched directly toward the sun at 1000 mph. How fast is its shadow moving?
    1. Define the rocket's velocity as \vec{R}.
    2. Define the unit normal on the ground (i.e. the direction of the rocket's shadow on the ground) as \vec{s}, where |\vec{s}| = 1.
    3. The dot product between \vec{R} and \vec{s} gives us the component of the velocity in the direction of \vec{s}. Thus, the speed of the shadow along the ground is given by
      v = \vec{R}\cdot\vec{s} = |\vec{R}|\,|\vec{s}| \cos(\theta)
    4. Intuition tells us that if the rocket is moving directly toward the sun, then its shadow does not appear to move at all. In this case, the angle between \vec{R} and \vec{s} is 90° and the dot product \vec{R}\cdot\vec{s} = |\vec{R}| \, |\vec{s}| \cos(\theta) = 0 since \cos(90^\circ)=0.
  • Consider the same rocket going parallel to the earth's surface. How fast is its shadow moving?
If the rocket is going parallel to the earth's surface, our intuition tells us that is shadow is moving at the same speed. This is confirmed by the mathematics, since the angle between the rocket's path and the ground is 0. Therefore, \vec{R}\cdot\vec{s} = |\vec{R}| \, |\vec{s}| \cos(\theta) =  |\vec{R}| \, |\vec{s}| = |\vec{R}| since \cos(0^\circ)=1 and |\vec{s}|=1.
  • What if the rocket were going at a 45° angle?
Our intuition tells us that the shadow will appear to move, but it will not be moving as fast as the rocket is moving. Mathematically, we have v_s = \vec{R}\cdot\vec{s} = |\vec{R}| |\vec{s}| \cos(45^\circ) = 500\sqrt{2} = 707.1 \, \textrm{mph}.

Vector Cross Product

The cross product of two vectors produces a vector perpendicular to the original two. The common right-hand rule is used to determine the direction of the resulting vector.

CrossProd.png

Assume that we have vectors a and b defined as


\begin{align}
 \vec{a} &= a_x \hat{i} + a_y \hat{j} + a_z \hat{k} \\
 \vec{b} &= b_x \hat{i} + b_y \hat{j} + b_z \hat{k}
\end{align}

where \hat{i}, \hat{j}, and \hat{k} represent the unit-normal vectors in the x, y, and z directions, respectively. The cross product of a and b is then defined as


 \vec{c} = \vec{a} \times \vec{b} = \left( a_y b_z -a_z b_y \right) \hat{i} + \left(a_z b_x -a_x b_z\right) \hat{j} + \left(a_x b_y -a_y b_x\right) \hat{k}

If we know the magnitude of a and b and the angle between them (θ), then the cross-product is given as


 \vec{a}\times\vec{b} = |\vec{a}|\,|\vec{b}|\,\sin(\theta) \, \vec{n},

where n is the unit-normal vector perpendicular to the plane defined by the vectors a and b.

For more information, see wikipedia's article on the cross product.

Vector Cross Product in MATLAB

   c = cross(a,b);  % c is a vector!


Matrix-Vector Product

The matrix-vector product is frequently encountered in solving linear systems of equations, although there are many other uses for it.

Consider a matrix A with m rows and n columns and vector b with n rows:


A=\left[\begin{array}{cccc}
 a_{11} & a_{12} & \cdots & a_{1n}\\
 a_{21} & a_{22} & \cdots & a_{2n}\\
 \vdots & \vdots & \ddots & \vdots\\
 a_{m1} & a_{m2} & \cdots & a_{mn}\end{array}\right],
\quad
b=
\left(\begin{array}{c}
 b_{1}\\ b_{2}\\ \vdots\\ b_{n}
\end{array}\right)

We may multiply A and b if the number of columns in A is equal to the number of rows in b. The product, c is a vector with m rows. The ith entry in c is given as


  c_i=\sum_{j=1}^n A_{i,j}b_j

For example, let's consider the following case:


A = \left[\begin{array}{ccc}
1 & 3 & 2 \\
2 & 5 & 4
\end{array}\right],
\quad
b = \left(\begin{array}{c}
 10 \\ 4 \\ 6
\end{array} \right)

Here we have m=2 and n=3. We define the product, c=Ab, as follows:

  1. Be sure that A and b are compatible for multiplication.
    • Since A has three columns and b has three rows, they are compatible.
  2. Determine how many entries are in c.
    • The c vector will have as many rows as A has.
    • Since A has two rows, c will have two rows.
  3. Determine the elements in c
    • Applying the above formula for i=1 we find
      
 \begin{align}
   c_1 &= \sum_{j=1}^{3} A_{1,j} b_j \\
       &= A_{1,1} b_1 + A_{1,2} b_2 + A_{1,3} b_3 \\
       &= 1 \cdot 10 + 3\cdot 4 + 2 \cdot 6 \\
       &= 34.
 \end{align}
    • Applying the above formula for i=2 we find
      
 \begin{align}
   c_2 &= \sum_{j=1}^{3} A_{2,j} b_j \\
       &= A_{2,1} b_1 + A_{2,2} b_2 + A_{2,3} b_3 \\
       &= 2 \cdot 10 + 4\cdot 4 + 4 \cdot 6 \\
       &= 64.
 \end{align}

Therefore, we have


 c=\left(\begin{array}{c} 34 \\ 64 \end{array}\right).

Geometric Interpretation of Matrix-Vector Product

Here is a great discussion of how a matrix-vector product can be interpreted geometrically.

Matrix-Vector product in MATLAB

   c = A*b;   % A is a matrix, b and c are vectors.


Matrix-Matrix Product

The matrix-matrix product is a generalization of the matrix-vector product. Consider the general case where we want to multiply two matrices, C = A B, where A and B are matrices. The following rules apply:

  1. The number of columns in A must be equal to the number of rows in B.
  2. The number of columns in C is equal to the number of columns in B.
  3. The number of rows in C is equal to the number of rows in A.

We can state this differently. Suppose that A has m rows and p columns, and that B has p rows and n columns. This satisfies rule 1 above. Then, according to rules 2 and 3, C has m rows and n columns.

Examples of valid and invalid matrices for multiplication
A B Valid? Comments
\left[\begin{smallmatrix} 1 & 2 \\ 3 & 4 \end{smallmatrix}\right] \left[\begin{smallmatrix} 5 & 6 \\ 7 & 8 \end{smallmatrix}\right] Yes Number of columns in A = 2, number or rows in B = 2.
\left[\begin{smallmatrix} 1 & 2 \end{smallmatrix}\right] \left[\begin{smallmatrix} 5 \\ 8 \end{smallmatrix}\right] Yes
\left[\begin{smallmatrix} 1 \\ 2 \end{smallmatrix}\right] \left[\begin{smallmatrix} 5 & 8 \end{smallmatrix}\right] Yes
\left[\begin{smallmatrix} 1 & 2 \end{smallmatrix}\right] \left[\begin{smallmatrix} 5 & 8 \end{smallmatrix}\right] No Number of columns in A=2, number of rows in B = 1. Therefore these cannot be multiplied.

The general formula for matrix multiplication is


  C_{i,j} = \sum_{k=1}^n A_{i,k} B_{k,j}

Let's show how this works for a simple example:


 A = \left[ \begin{array}{cc}
   1 & 2 \\ 3 & 5 
  \end{array} \right],
 \quad
 B = \left[ \begin{array}{ccc}
   2 & 4 & 10 \\ 6 & 3 & 1
  \end{array} \right].

We can form C=A B since rule 1 is satisfied. From the size of A and B we define m=2 and n=3 and from rules 2 and 3, we conclude that the resulting matrix, C has 2 rows and 3 columns.

We are now ready to build each element of C.

Row Column Result
1 1 
  C_{1,1} = \sum_{k=1}^2 A_{1,k} B_{k,1}
          = A_{1,1} B_{1,1} + A_{1,2} B_{2,1}
          = 1 \cdot 2 + 2 \cdot 6
          = 14.
2 1 
  C_{2,1} = \sum_{k=1}^2 A_{2,k} B_{k,1}
          = A_{2,1} B_{1,1} + A_{2,2} B_{2,1}
          = 3 \cdot 2 + 5 \cdot 6
          = 36.
1 2 
  C_{1,2} = \sum_{k=1}^2 A_{1,k} B_{k,2} 
          = A_{1,1} B_{1,2} + A_{1,2} B_{2,2} 
          = 1 \cdot 4 + 2 \cdot 3 
          = 10.
2 2 
  C_{2,2} = \sum_{k=1}^2 A_{2,k} B_{k,2} 
          = A_{2,1} B_{1,2} + A_{2,2} B_{2,2} 
          = 3 \cdot 4 + 5 \cdot 3 
          = 27.
1 3 
  C_{1,3} = \sum_{k=1}^2 A_{1,k} B_{k,3} 
          = A_{1,1} B_{1,3} + A_{1,2} B_{2,3} 
          = 1 \cdot 10 + 2 \cdot 1 
          = 12.
2 3 
  C_{2,3} = \sum_{k=1}^2 A_{2,k} B_{k,3} 
          = A_{2,1} B_{1,3} + A_{2,2} B_{2,3} 
          = 3 \cdot 10 + 5 \cdot 1 
          = 35.

Now, putting all of these together, we have


 C=\left[\begin{array}{ccc}
  14 & 10 & 12 \\
  36 & 27 & 35
\end{array}\right].

Matrix-Matrix product in MATLAB

   C = A*B;   % A, B, and C are all matrices.

The MATLAB code to implement the above example is

   A = [ 1 2; 3 5; ];
   B = [ 2 4 10; 6 3 1; ];
   C = A*B;

Linear Systems of Equations

Any set of linear systems may be written in the form [A](x)=(b), where [A] is a matrix and (x) and (b) are vectors. In general, for a system of equations with n unknowns x1 ... xn, can be written as


 \begin{align}
   a_{1,1} x_1 + a_{1,2} x_2 + a_{1,3} x_3 + \cdots + a_{1,n} x_n &= b_1, \\
   a_{2,1} x_1 + a_{2,2} x_2 + a_{2,3} x_3 + \cdots + a_{2,n} x_n &= b_2, \\
   &\vdots \\
   a_{n,1} x_1 + a_{n,2} x_2 + a_{n,3} x_3 + \cdots + a_{n,n} x_n &= b_n, \\
 \end{align}

where a_{i,j} and b_i are numbers and x_i is the set of unknown variables. These equations can now be written in matrix format:


 \left[ \begin{array}{ccccc}
   a_{1,1} & a_{1,2} & a_{1,3} & \cdots & a_{1,n} \\
   a_{2,1} & a_{2,2} & a_{2,3} & \cdots & a_{2,n} \\
   \vdots  & \vdots  & \ddots  & \vdots & \vdots \\
   a_{n,1} & a_{n,2} & a_{n,3} & \cdots & a_{n,n}
 \end{array} \right]
 \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right)
 =
 \left( \begin{array}{c} b_1 \\ b_2 \\ \vdots \\ b_n \end{array} \right).

Example

Consider the following equations:


 \begin{cases}
   x + y -5z = 9 \\
   6-z = 2y \\
   y-3x = 10z \\
 \end{cases}

To put these in matrix format, we first choose the solution vector. Let's choose to order this as x, y, z. Now we rewrite the original system so that we have a coefficient multiplying each of these variables, and then we put the "constant" on the right-hand-side:


 \begin{align}
    1x + 1y -5z  &=  9 \\
    0x - 2y -1z  &= -6 \\
   -3x + 1y -10z &=  0
 \end{align}

We may then identify the A matrix and b vectors and write


 \left[\begin{array}{ccc}
  1 &  1 & -5  \\
  0 & -2 & -1  \\
 -3 &  1 & -10 \\
 \end{array} \right]
 \left( \begin{array}{c}
  x \\ y \\ z
 \end{array}\right)
 =
 \left( \begin{array}{c}
  9 \\ -6 \\ 0
 \end{array}\right)

Note that if we multiply the above matrix and vector, we recover the original set of equations (try it to be sure that you can do it).

Note that if we chose a different variable ordering, then it simply rearranges the columns in the matrix. For example, if we order the variables z, x, y then we have


 \begin{align}
    -5z  + 1x + 1y &=  9 \\
    -1z  + 0x - 2y &= -6 \\
    -10z - 3x + 1y &=  0
 \end{align}
 \quad \Rightarrow \quad
 \left[\begin{array}{ccc}
  -5  &  1 &  1 \\
  -1  &  0 & -2 \\
  -10 & -3 &  1 \\
 \end{array} \right]
 \left( \begin{array}{c}
  z \\ x \\ y
 \end{array}\right)
 =
 \left( \begin{array}{c}
  9 \\ -6 \\ 0
 \end{array}\right)


This page discusses a few methods to solve linear systems of equations, Ax=b. There are two categories of solution methods for linear systems of equations: direct solvers and iterative solvers. These will be discussed briefly below.


Direct Solvers

Direct solvers solve the linear system to within the numeric precision of the computer. There are many different algorithms for direct solution of linear systems, but we will only mention two: Gaussian elimination and the Thomas Algorithm.


Gaussian Elimination

Gaussian elimination is a methododical approach to solving linear systems of equations. It is essentially what you would do if you were to solve a system of equations "by hand" but is systematic so that it can be implemented on the computer. For information on the algorithm, see this wikipedia article.

The Thomas Algorithm (Tridiagonal Systems)

The Thomas Algorithm is used for solving tridiagonal systems of equations, Ax=b. A tridiagonal matrix is one where the only non-zero entries in the matrix are on the main diagonal and the lower and upper diagonals. In general a tridiagonal matrix can be written as

A = \left[\begin{array}{ccccc}
  d_1    & u_1    & 0      & \cdots     & 0      \\ 
  \ell_1 & d_2    & u_2    & 0          & \cdots \\
         & \ddots & \ddots & \ddots     &        \\
   0     & 0      & \cdots & \ell_{n-1} & d_n
\end{array}\right]

Here d_i represents the diagonal of the matrix, \ell_i represents the lower-diagonal, and u_i represents the upper-diagonal.

In terms of this nomenclature, the Thomas Algorithm can be written as

forward pass...
for i = 2 ... n
  d(i) = d(i) - ( u(i-1)*ell(i-1) ) / d(i-1)
  b(i) = b(i) - ( b(i-1)*ell(i-1) ) / d(i-1)
end
x(n) = b(n) / d(n)
backward pass...
for i = n-1 ... 1
  x(i) = ( b(i) - u(i) * x(i+1) ) / d(i)
end

Iterative Solvers

Iterative solvers are fundamentally different than direct solvers. They require an initial guess and then refine the solution until a tolerance is reached.


Jacobi Method

The Jacobi method for solving linear systems is the simplest iterative method. In this case, each iteration of the solution involves cycling through the equations and solving the ith equation for the xi, using values for all other xj from the previous iteration.

The algorithm for the Jacobi method is:

xnew = x;
err = tol+1;
while err > tol
  for i=1:1:n
     update the ith value of x.
     tmp = 0
     for j=1:1:i-1
        tmp = tmp + A(i,j) * x(j)
     end
     for j=i+1:1:n
        tmp = tmp + A(i,j) * x(j)
     end
     xnew(i) = x(i) + ( b(j) - tmp ) / A(i,i)
  end
  x = xnew
  check for convergence - set err
end

For example, assume that we have the following system of equations:


\left[ \begin{array}{rrrr}
 -3 &  1 &  0 &  0  \\
  1 & -2 &  1 &  0  \\
  0 &  1 & -2 &  1  \\
  0 &  0 &  1 & -3 
\end{array} \right]
\left( \begin{array}{c}
 x_1 \\ x_2 \\ x_3 \\ x_4
\end{array} \right)
=
\left( \begin{array}{r}
 -4 \\ 2 \\ 9 \\ 2
\end{array} \right)

If we set our initial guess as x=\left[ \begin{array}{cccc} 1 & 1 &
1 & 1 \end{array} \right]^T then for the first iteration we proceed as follows:

  • Solve the first equation for x1 to find  x_1 = \frac{-4 -x_2}{-3} = 5/3.
  • Solve the second equation for x2 to find  x_2 = \frac{2 - x_1-x_3}{-2} = 0.

We repeat this for all of the equations. At the end of the first iteration, we have x=\left(\begin{array}{cccc} 5/3, & 0, & -7/2, & -1/3 \end{array}\right).

We are now ready to start with the second iteration, where we use the newly computed values for x_i. The table below shows the results for the first few iterations.

The last entry in the table is the analytic solution to this system of equations. Note that after 40 iterations we have converged to the exact answer to within an absolute error x-xexact = 10-3.

See here for a Matlab implementation of the Jacobi method. More information can be found here.


Gauss-Seidel Method

The Gauss-Seidel method is a simple extension of the Jacobi method. The only difference is that rather than using "old" solution values, we use the most up-to-date solution values in our calculation.

The algorithm for the Gauss-Seidel method is:

err = tol+1
while err > tol
  for i=1:1:n
     update xi with latest information on all other values of x
     for j=1:1:n
        tmp = tmp + A(i,j) * x(j)
     end
     x(i) = x(i) + ( b(i) - tmp ) / A(i,i)
  end
  check for convergence - set err
end

If we apply the Gauss-Seidel method to the example above, we proceed as follows for the first iteration:

  • Solve the first equation for x1 to find x_1 = \frac{-4 -x_2}{-3} = 5/3
  • Solve the second equation for x2 and use the most recent value of x1 to find x_2 = \frac{2 - x_1-x_3)}{-2} = \frac{2-5/3-1}{-2} = 1/3

We repeat this for all of the equations. At the end of the first iteration, we have x=\left(\begin{array}{cccc} 5/3, & 1/3, & -23/6, & -35/18, \end{array}\right)

The following table shows the convergence results for the Jacobi and Gauss-Seidel methods. The Gauss-Seidel method converges much faster than the Jacobi method.

Convergence for the Jacobi and Gauss-Seidel iterative linear solvers
Jacobi Method Gauss-Seidel Method
Iteration x_1 x_2 x_3 x_4 x_1 x_2 x_3 x_4
1 5/3 0 -7/2 -1/3 5/3 1/3 -23/6 -35/18
2 4/3 -23/12 -14/3 -11/6 1.4444 -2.1944 -6.5694 -2.8565
3 25/26 -8/3 -51/8 -20/9 0.6019 -3.9838 -7.9201 -3.3067
4 0.4444 -3.8403 -6.9444 -2.7917 0.0054 -4.9574 -8.6320 -3.5440
5 0.0532 -4.2500 -7.8160 -2.9815 -0.3191 -5.4756 -9.0098 -3.6699
10 -0.5170 -5.7294 -9.0648 -3.6601 -0.6719 -6.0377 -9.4194 -3.8065
20 -0.6803 -6.0484 -9.4218 -3.8061 -0.6875 -6.0625 -9.4375 -3.8125
-11/16 = -0.6875 -97/16 = -6.0625 -151/16 = -9.4375 -61/16 = -3.8125

More information on the Gauss-Seidel method can be found here.


Some Caveats

A direct solver will always find a solution if the system is well-posed. However, iterative solvers are more fragile, and will occasionally not converge. This is particularly true for the Jacobi and Gauss-Seidel solvers discussed above. In fact, the matrices must have special properties for these solvers to converge. Specifically, they must be symmetric and positive-definite or diagonally dominant. What this means in practice is that in a given row of the matrix, the diagonal coefficient is larger than the sum of the off-diagonal coefficients,

\left| a_{ii} \right| > \sum_{j\ne i} \left| a_{ij} \right|

Other Iterative Methods

There are many other iterative solvers that are much more sophisticated and robust. Among them are the Conjugate-Gradient and GMRES methods.


Convergence

Iterative solvers require us to decide when to stop iterating. In other words, we need to decide when the solution is converged. We do this by comparing an error measure to a tolerance.

There are a few common ways to determine convergence for a linear system:

  • \epsilon = \left\Vert b-Ax \right\Vert
  • \epsilon = \left\Vert x^{k}-x^{k-1}\right\Vert

where the terms \left\Vert a \right\Vert indicate an appropriate vector norm.

See the discussion on iteration and convergence for more information on this topic.