Guide to Computational Linear Algebra

[jemidiah]

0.1 Introduction

I decided to write this guide in the spirit of VBAhack's Cubic Spline Tutorial. Overall it attempts to be a self-contained (where possible) guide for implementing basic operations in linear algebra. Part 1 is an introduction to vectors and matrices. It can be skipped, but only if you're good at pattern matching or have done some linear algebra before. Part 2 covers the basic uses of linear algebra with an eye towards computer graphics. This guide is generally meant to be used in pieces--for instance, 2.3 and 1.5 together tell you how to rotate points in 2D using matrix multiplication. It would likely be too much information to read straight through if you are completely unfamiliar with linear algebra, though anyone interested is welcome to try.

This guide is meant to be different things to different people. It's meant to be a conceptual introduction to vectors and matrices, with many of the details left out for other sources to fill in. It's meant to be a search result when you search the forums for, say, "dot products", where you'll get a formula for computing dot products as well as an explanation of what they are.

This guide tries to assume very little about your mathematical background. You should be familiar with basic algebra--addition, multiplication, and so on, using symbols. Aside from that, your background is much less important than your ability to reason and deal with some symbolic notation, since (if you're reading this) the ideas should be (largely) new anyway. Enjoy!


0.2 Table of Contents

PART 1: Background and Definitions

1. What is Linear Algebra?
2. From Numbers to Vectors
3. Basic Algebra with Vectors
4. From Vectors to Matrices
5. Basic Algebra with Matrices

PART 2: Basic Computer Graphics with Linear Algebra

1. Dot Products
2. Cross Products
3. 2D Rotation Matrix
4. 3D Rotation Matrices

[jemidiah]

PART 1: Background and Definitions

1.1 What is Linear Algebra?

Without getting very far into details, linear algebra is a branch of math that usually comes right after Calculus. It deals with points in higher dimensions, and so-called "linear" transformations of those points. It turns out that rotations can be represented as linear transformations, which is one of several reasons linear algebra is used in programming. For instance, I've heard graphics cards described as giant matrix multipliers.

Some of the buzzwords of linear algebra are: vector, matrix, vector space, subspace, eigenvalue, eigenvector, eigenspace, orthogonal vectors, linear (in)dependence, dot product, cross product, inner product, and ice cream (:P). The first two in that list are probably the most important, so this guide builds up each in turn, starting from basic numbers.

For more details on the scope of linear algebra, see the Wikipedia article, or find a book on it.


1.2 From Numbers to Vectors

A number can be considered as a single point on an infinitely long line. This is useful, but we can do more by combining multiple numbers. For instance, two Cartesian Coordinates let us find points in a plane (see first attachment).

In two dimensions, we usually use the notation p = (x, y) for a point p. The green point has p = (2, 3). Vectors can be considered as points, with a direction. By convention, vectors point from the origin--(0, 0) in 2D--towards the point in question. Converting p from a point to a vector, we have p = <2, 3>. The vector <-3, 1> is drawn in red, with its direction given by the arrow head. In general, a vector is just a list of numbers that specify a single point in space, which points from the origin to that point.

In programming, vectors are data structures which have components that are numbers. The above vector p has x-component 2, y-component 3. After 3 dimensions, naming each component gets tiresome, so we usually just number them.

[Note: In this guide, all numbers are real numbers unless otherwise specified. Complex matrix theory can be fascinating, but I don't discuss it.]


1.3 Basic Algebra with Vectors

With single numbers, we can do arithmetic. This turns out to be useful, for instance, in figuring out how much your dinner bill was last night. We can do arithmetic with vectors as well. The following each have intuitive meaning that I usually won't get in to; create it yourself if you're interested, or look it up.

For simplicity, I'll use object-oriented programming notation. The following assumes you have a Vector object, with as many member coordinates as your application needs. For instance, the following definitions use 2D vectors p=<p.x, p.y>, q=<q.x, q.y>; 3D vectors r=<r.x, r.y, r.z>, s=<s.x, s.y, s.z>; and general vectors of the same, arbitrary dimension u and v. I also use the arbitrary number c, which is not a vector.

Equality
2D: p=q means p.x=q.x and p.y=q.y
3D: r=s means r.x=s.x, r.y=s.y, and r.z=s.z
nD: u=v means, for each coordinate i, u.i=v.i

Addition
2D: p+q = <p.x+q.x, p.y+q.y>
3D: r+s = <r.x+s.x, r.y+s.y, r.z+s.z>
nD: u+v = <u.1+v.1, u.2+v.2, ..., u.n+v.n>

Negation
2D: -p = <-p.x, -p.y>
3D: -r = <-r.x, -r.y, -r.z>
nD: -u = <-u.1, ..., -u.n>

Subtraction
Vectors are implicitly centered at the origin, <0, 0, ..., 0>. Subtracting a vector p from a vector o shifts the origin to vector o. That is, (p-o) can be thought of as a vector pointing from o to p.
2D: p-q = <p.x-q.x, p.y-q.y>
3D: r-s = <r.x-s.x, r.y-s.y, r.z-s.z>
nD: u-v = u+(-v) = <u.1-v.1, ..., u.n-v.n>

Scaling
2D: cp = c*p = <c*p.x, c*p.y>
3D: cr = c*r = <c*r.x, c*r.y, c*r.z>
nD: cu = c*u = <c*u.1, ..., c*u.n>
[Note: for c != 0, scalar division is multiplying by a reciprocal:
2D: p/c = <p.x/c, p.y/c>
3D: r/c = <r.x/c, r.y/c, r.z/c>
nD: u/c = (1/c)u = <u.1/c, ..., u.n/c>]

Multiplication and Division
These are strange cases, and can be done in several ways. For now, suffice it to say that both "dot products" and "cross products" are ways of multiplying two vectors, and that division of vectors depends on the choice of multiplication, and even then may not make sense using a traditional notion of division as the opposite of multiplication.

Length
The Pythagorean Theorem gives us the length of vectors in 2D directly. We can use the theorem several times to give length in higher dimensions. There are many notations for length; I give several of the common ones. The length of a vector is also called its norm.
2D: Length of p = |p| = ||p|| = p = Sqrt(p.x² + p.y²)
3D: Length of r = |r| = ||r|| = r = Sqrt(r.x² + r.y² + r.z²)
nD: Length of u = Sqrt(Sum_i(u.i²))

Normalizing
Unit vectors are vectors of length 1. All non-zero vectors can be scaled up or down to have such a length while keeping their direction constant. The common notation for a unit vector is to put a "hat" on a vector symbol. I don't want to type that, so I'll use the made-up function "unit()", which accepts a non-zero vector and gives a unit vector version of the input. This output is also called a normalized vector, and this operation is also called norming or normalizing. Informally, a unit vector keeps track of only direction information. Also, if |u| = 1, then unit(u) = u.
2D: unit(p) = p/|p| = 1/Sqrt(p.x² + p.y²) * <p.x, p.y>
3D: unit(r) = r/|r| = 1/Sqrt(r.x² + r.y² + r.z²) * <r.x, r.y, r.z>
nD: unit(u) = u/|u| = 1/Sqrt(Sum_i(u.i²)) * <u.1, ..., u.n>
[Note: this operation is one application of a magical Quake III function]

[jemidiah]

1.4 From Vectors to Matrices

We'll find in Section 2 how useful vectors can be. Extending the idea of a vector as a point in space, we can consider a matrix as an array of points in space. This generates a 2D array of numbers, with some number of rows, call it m, and some number of columns, call it n. The dimensions of such a matrix are mxn. The following is a 3x4 matrix, which can be considered to be either three vectors of dimension four (one vector to a row), or four vectors of dimension three (one vector to a column):

┌                 ┐
│  0   2   0   0  │
│  1 0.5   4 1/9  │
│ -1   1   2   3  │
└                 ┘

I have trouble remembering if mxn means "m rows, n columns" or "m columns, n rows". Just memorize "rows by columns", and you can figure it out from there.

Usually, you'll see matrices written with capital letters. You access the elements of a matrix, A, with two indicies, usually written as i and j. i corresponds to the row index, j corresponds to the column index. So, A_i is a row vector, and A_i,j is an element. Supposing the above example is the matrix A, A_2,3 = 4. In general, for an mxn matrix, A_1,1 is the upper-left entry, and A_m,n is the lower-right entry.

In math, vectors are actually either column vectors or row vectors. A column vector is an mx1 matrix, while a row vector is a 1xn matrix. By convention, unless otherwise specified, vectors are almost always column vectors. In particular, the vectors in this guide are column vectors. You may turn a column vector into a row vector (and vice-versa) by transposing it (see Section 1.5). In programming, you have two options of how to mix vectors and matrices in code. You can implement them separately, having both Vector objects (or 1D arrays) and Matrix objects (or 2D arrays), or you can make all of your vectors mx1 Matrix objects corresponding to column vectors.


1.5 Basic Algebra with Matrices

As with vectors, defining various operations on matrices is useful. This section uses an object-oriented Matrix object. This object can be implemented in several ways. Regardless, I'll simply use row-column indecies, with A(i,j) giving the element at row i, column j. This section shows how to compute an element at row i, column j using the given operations. To extend this to a matrix, simply loop over i and j, as the example pseudocode shows. Again, each operation has an intuitive meaning, though it's irrelevant in this context. Very often, you'll be working with square matrices, which somewhat simplifies the following.

Most operations are defined component-wise, require one or two matrices of the same size as inputs, and return a matrix of the same size as output. The exceptions are noted. Whenever an explicit formula or algorithm is not given, it is either too complex or tedious to include in this introductory explanation. See the appropriate links for more.

The following definitions use several common variables:
Let m, n, and p be positive whole numbers.
Let A and B be mxn matrices.
Let Z be an nxp matrix.
Let S be an nxn square matrix.
Let c be an arbitrary number.
Let u be a column vector of dimension n.

Equality
A = B means, for all indecies i and j, A(i,j) = B(i,j)

Addition
If M = A+B, then
M(i,j) = A(i,j) + B(i,j)

Negation
If M = -A, then
M(i,j) = -A(i,j)

Subtraction
If M = A-B, then
M(i,j) = A(i,j) - B(i,j)

Scaling
If M = cA = c*A, then
M(i,j) = c*A(i,j)

Transposition
Transposing an mxn matrix results in an nxm matrix that is "flipped"--the rows and columns are inverted. While this seems almost silly, it's remarkably useful.
If M = A^T, then
M(i,j) = A(j,i)

Multiplication
Multiplication requires two matrices of "compatible" dimensions: the first matrix must have the same number of columns as the second matrix has rows. This is because it's most naturally defined in terms of dot products, which require vectors of the same dimension. However, I won't get in to an intuitive explanation of matrix multiplication here, since it's unnecessary when writing a program. Multiplying an mxn matrix by an nxp matrix results in an mxp matrix. See the example for more.
If M = AZ = A*Z, then
M(i,j) = Sum_{k=1 to n} A(i,k)*B(k,j)

Multiplication by a Column Vector
Multiplication of a matrix by a (column) vector is very common, and is a special case of multiplication by a matrix. Using the rules above, it returns a column vector, possibly of a different length, which is given by the number of rows in the matrix. In particular, using the notation of Section 1.3 for a Vector object, the following formula arises:
v = Au = A*u = <Sum_k=1...n A(1,k)*u.k, Sum_k=1...n A(2,k)*u.k, ..., Sum_k=1...n A(m,k)*u.k>
Since A is mxn, and u is nx1, v is mx1.

Determinant
The determinant of a matrix gives important information about how it transforms sets of (column) vectors under multiplication. It also determines when a matrix is invertible (see below). However, any serious explanation of the determinant is too complex to include here. It may be calculated using the Laplace Expansion (also known as Cofactor Expansion). This is grossly inefficient on large matrices, and can be improved in a number of ways.

Inversion
Matrix inversion is difficult, slow, and may not be possible. It requires a square matrix. In general, matrix inversion tries to "undo" the operation of multiplying a matrix by a vector. For instance, suppose Su = v. Ideally we would like an inverse of S, denoted by S^-1 (also an nxn matrix), to reverse this equation. That is, we want u = S^-1v. As it turns out, a matrix S is invertible if and only if the determinant of S is not zero, which of course has intuitive meaning that's too complicated to include here. There are several algorithms for matrix inversion. 2x2 matrices have particularly simple inverses. For larger matrices, one may use Gauss-Jordan Elimination (the basic "industry standard") or Cramer's Rule (commonly considered slow for n>3 or so). A matrix may also be inverted recursively using the blockwise formula for the matrix inverse, which uses block matrices.

Example 1:

function MatrixAddition(A as mxn Matrix, B as mxn Matrix) returns mxn Matrix:
    Make Sum a new mxn Matrix
    For i = 1 to m     'loop over Sum's rows
        For j = 1 to n 'loop over Sum's columns
            'Compute component-wise sum
            Sum(i,j) = A(i,j) + B(i,j)
        Next j
    Next i
    return Sum
end function

Example 2:

function MatrixMultiplication(A as mxn Matrix, Z as nxp Matrix) returns mxp Matrix:
    Make Prod a new mxp matrix
    For i = 1 to m     'loop over Prod's rows
        For j = 1 to p 'loop over Prod's columns
            'Compute dot product
            dot = 0
            For k = 1 to n
                dot = dot + A(i,k)*B(k,j)
            Next k
            
            Prod(i,j) = dot
        Next j
    Next i
    return Prod
end function

[jemidiah]

PART 2: Basic Computer Graphics with Linear Algebra

2.1 Dot Products

The dot product lets us find angles between vectors easily. To be explicit, take the following example. Suppose we have a turret in a 2D game. Say that the turret can only turn through 90 degrees: 45 degrees left from its starting position, and 45 degrees right. We can remember which way the turret was pointing initially through the use of a unit vector (see Section 1.3), which would probably be set during level design. We can also get a vector pointing from the turret to an enemy easily with vector subtraction. Now, we need to know what the angle is between the starting position vector and the enemy, to see if the angle is too large. The dot product lets us find this angle easily without resorting to drawing out the situation to analyze with geometry and trigonometry.

The dot product of two vectors u and v of the same length is given as follows (the . should be centered in the middle vertically as well as horizontally, but a period is easier for typesetting this):
2D: u . v = <u.x*v.x, u.y*v.y>
3D: u . v = <u.x*v.x, u.y*v.y, u.z*v.z>
nD: u . v = <u.1*v.1, ..., u.n*v.n>

The defining property of the dot product can be derived from basic geometry and trigonometry. It relates the lengths of each of the input vectors and their dot product to the angle between the vectors. There are actually two angles between vectors--one small one and one large one, which add up to 360 degrees. The dot product deals with the smaller angle. This angle is traditionally called "theta", after the Greek letter. The relation between these quantities is, using the notation of Section 1.3:

u . v = |u|*|v|*cos(theta)

If vectors p and q are unit vectors, then by definition their lengths are 1, and the angle can be found very easily:

theta = arccos(p . q)

Using the turret example, the following pseudocode would tell us if the target was in range or not:

'Fire the turret, if enemy is within turning range. Return the new enemy health.
function TurretAttack(InitialDirection as Unit 2DVector, TurretPosition as 2DVector,
                       EnemyPosition as 2DVector, AttackStrength as Number,
                       EnemyHealth as Number) returns Number:

    'Use dot product to find angle from turret to enemy
    TurretToEnemy = unit(EnemyPosition - TurretPosition) 'see Section 1.3 for unit()
    theta = Math.acos(InitialDirection . TurretToEnemy)
    
    'Attack if possible
    If |theta| <= 45 degrees Then
        'In range
        return EnemyHealth - AttackStrength
    Else
        'Out of range
        return EnemyHealth
end function

2.2 Cross Products

The cross product is only defined in 3D (and, incidentally, 7D, though this description doesn't consider the 7D version, because honestly you likely couldn't care less). It lets us find a vector mutually perpendicular to two other vectors. This is mostly useful to find a vector normal to a surface in 3D.

As a motivating example, suppose you have three points in 3D which form a triangle. Take the points as the vectors p, q, and r. These points must be in some plane (the triangle is in that plane as well). We want a vector pointing directly out of this plane. Such a vector is useful, for instance, in collision physics, since it tells us in which direction an object bounces off the triangle.

The cross product takes two vectors that lie in a plane and returns a vector pointing directly out of that plane. In this example, two such vectors are u = (p-q) and v = (p-r) (see Section 1.3, vector subtraction, for more). Their cross product, which points out of the plane, is given by

u x v = <u.y*v.z-u.z*v.y, u.z*v.x-u.x*v.z, u.x*v.y-u.y*v.x>

There are a few edge cases to this example. If u and v point in the same (or directly opposite) direction, their cross product will be the zero vector. If either vector is the zero vector, their cross product is again the zero vector. Often, you only want direction information out of the cross product. By convention, a normal vector is usually normalized (see Section 1.3 and the unit() function). The length of the cross product does carry some information, though, and is similar to the dot product. The length of the input vectors and the angle theta between them is related to the length of their cross product by the following:

|u x v| = |u|*|v|*sin(theta)

This can be interpreted as the area of the parallelogram the two input vectors determine.

[jemidiah]

2.3 2D Rotation Matrix

Suppose you're simulating a 2D analog clock. You want the hour, minute, and second hands to move around the face of the clock. There are several ways to do this, though a 2D rotation matrix works well. You start with, say, the minute hand pointing straight up. This is represented by the 2D vector m, which contains both the length and initial direction information for the minute hand. Using Section 1.4, we'll say m is a column vector. Using Section 1.5, we can multiply a matrix and a column vector to get back a column vector. It turns out that multiplying m by a certain matrix returns a new vector, which is just m rotated by an angle, called theta (which is the traditional Greek letter name). The length is kept constant. This matrix is

           ┌                         ┐
           │ cos(theta)  -sin(theta) │
R(theta) = │ sin(theta)   cos(theta) │
           └                         ┘

That is, for any initial vector m and any angle theta, the vector v = R(theta)*m is a copy of m, rotated by angle theta clockwise or counterclockwise, depending on the coordinate system. Often in programming, the y-axis increases going down instead of up. In this case, R(theta) performs a rotation clockwise as shown in the attached Fig 2.3 (source). This is the type of rotation we want on the face of a clock. The following pseudocode would then give the current position vector for the minute hand:

function getMinuteHand(minutes as Number, InitialPosition as 2D Vector) returns 2D Vector:
    angle = minutes * 6 '6 degrees per minute with 360 degrees on a clock face
    return 2DRotation(angle) * InitialPosition 'matrix-vector multiplication;
                                               '2DRotation(theta) returns the 2x2 matrix above.
                                               'Theta is in degrees here, though most computers
                                               'need radians as input to trig functions
end function

2.4 3D Rotation Matrices

Rotations in 3D are much more complicated than in 2D, since the extra dimension allows an infinite number of rotation types. First, you have to choose an axis to rotate around. You then have to pick how far you want to rotate, by choosing an angle. See Fig 2.4 (modified from here). Simple rotations can be done around the x, y, or z axis. These can be done with the following three matrices, which rotate about the indicated axis in a counterclockwise direction (as drawn in Fig 2.4). This is sometimes called a right-handed rotation, since you can point your right thumb along the axis and curl your fingers in the direction of the rotation.

Rotation about the x-axis by angle theta:

            ┌                                       ┐
            │ 1             0            0          │
Rx(theta) = │ 0             cos(theta)   sin(theta) │
            │ 0            -sin(theta)   cos(theta) │
            └                                       ┘

Rotation about the y-axis by angle theta:

            ┌                                       ┐
            │  cos(theta)   0           -sin(theta) │
Ry(theta) = │  0            1            0          │
            │  sin(theta)   0            cos(theta) │
            └                                       ┘

Rotation about the z-axis by angle theta:

            ┌                                       ┐
            │  cos(theta)   sin(theta)   0          │
Rz(theta) = │ -sin(theta)   cos(theta)   0          │
            │  0            0            1          │
            └                                       ┘

These can be done in sequence through matrix multiplication. That is, to rotate the (column) vector r₀ around the x-axis 45 degrees, and to then rotate the result around the y-axis 13 degrees, you would compute the product Ry(13 deg)*Rx(45 deg)*r₀. You may combine matrices in other, similar ways to simulate roll-pitch-yaw or more generally Euler angles.

You may also rotate about an arbitrary axis, represented by the vector n, using a (significantly more hideous) matrix. This formula requires that n is a unit vector (see Section 1.3). It rotates in the same counterclockwise/right-handed direction as the above matrices.

Rotation about an arbitrary axis, n, by angle theta:

Let c = cos(theta)
Let s = sin(theta)

            ┌                                                                   ┐
            │  n.x^2+(1-n.x^2)*c    n.x*n.y*(1-c)-n.z*s   n.x*n.z*(1-c)+n.y*s   │
Tn(theta) = │  n.x*n.y*(1-c)+n.z*s  n.y^2+(1-n.y^2)*c     n.y*n.z*(1-c)-n.x*s   │
            │  n.x*n.z*(1-c)-n.y*s  n.y*n.z*(1-c)+n.x*s   n.z^2+(1-n.z^2)*c     │
            └                                                                   ┘

[jemidiah]

[Reserved for possible expansion]