Linear Algebra

Intro

Basic idea behind linear algebra: solving large number of linear algebraic equations simultaneously.

$ a_{i,1} x_1 + a_{i,2} x_2 + a_{i,3} x_3 + \dots + a_{i,N} x_{N} = b_{N} \qquad i=1 \dots M $

There are N unknowns and M equations, forming a linear system from matrices:

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

where lowercase letters denote vectors, and the dot represents multiplication of two matrices or a matrix and a vector, or dot product of two vectors.

Nonsingular vs Singular

We are not guaranteed to find a solution to a given linear system Ax=b

• row degeneracy - one or more of the M equations is a linear combination of other equations
• column degeneracy - one or more of the N variables are interchangeable and cannot be distinguished (are non-unique)
• for square matrices, these two conditions are equivalent

Having singular system of equations means we can't find good solutions.

Other things preventing good solution:

• Equations being equivalent to within roundoff error (approximately degenerate)
• Accumulation of roundoff error during solution can swamp solution

Many linear equation solving packages dedicate a good chunk of their code to solving these problems

Sense of scale:

• Linear equations with 20-50 unknowns are routine, can almost always be solved in single precision with straightforward routines
• Linear equations with up to 1,000 unknowns require double precision
• Larger sets of thousands or millions of equations can be solved with coefficients are sparse (mostly zero) by using specialty sparse methods

Tasks of Computational Linear Algebra

Beyond simply solving Ax=b for the unknowns x, we may want to do other things:

• Solve more than one matrix equation, $ \mathbf{A \cdot x}_j = \mathbf{b}_j $, where the matrix A does not change but unknown/RHS vectors do change
• Calculating the inverse of the matrix
• Calculating the determinant of a square matrix
• Calculating the condition number of a matrix

If M < N (row degeneracy), or if M = N and equations are degenerate, there are fewer equations than unknowns. If fewer equations than unknowns:

• Either no solution vector, or
• More than one solution vector
• If more than one solution, solution space consists of particular solution added to linear combination of N-M vectors (which are in the nullspace of the matrix A)
• SVD can be used to solve this problem, by find these vectors (SVD finds the nullspace of a matrix A)

If M > N (column degeneracy), there are fewer unknowns than equations. If fewer unknowns than equations:

• Think least squares - more data points than unknowns
• Looking for best fit/compromise solution
• Linear least-squares problem can be used to solve overdetermined linear system

Linear Algebra Software

Software:

• Blas is short for "basic linear algebra subroutines" and provides the linear algebra operations that compose Lapack
• Lapack is the most common linear algebra package (Linpack is old and was replaced by Lapack)
• ScaLapack is version of Lapack for parallel architectures
• Lapack routines divided by precision, algorithm, simplifications (tridiag., sparse, dense, etc.)
• Atlas is version of Lapack that automatically tunes itself to run as fast as possible for your system's architecture when it is compiled

Direct vs Iterative:

• Major dividing line between routines
• Direct - algorithms that execute in a predictable number of operations
• Iterative - attempt to converge on the desired answer to within a desired tolerance, taking as many steps as necessary

Gauss-Jordan Elimination

Classic technique for solving linear equations: combine equations to cancel out unknowns, turning A into the identity matrix

Gauss-Jordan elimination produces the solution of the equations and the matrix inverse $ \mathbf{A}^{-1} $ - this is an expensive extra thing to produce, not always desirable to have it

Deficiencies:

• All right hand side values must be stored and manipulated at once (not feasible for very large systems of equations)
• Gauss-Jordan is three times slower than alternative techniques that do not produce the matrix inverse

Advantages:

• Straightforward
• Solid
• Pivoting

Pivoting:

• Pivoting is just interchanging rows (partial pivoting) or rows and columns (full pivoting) to put desired elements on matrix diagonals
• Partial pivoting is easier, full pivoting requires keeping track of permutations, and un-doing these when solution vector achieved

Storage requirements of method:

• The matrix inverse of A is gradually built up in A, as original A is overwritten
• Likewise solution vector x can gradually replace right-hand side vector b (solution vector is overwritten one row element at a time, each time an element is overwritten it will not be used again)

On the subject of row vs. column elimination strategy:

• if we perform row operations, row transformations build right to left, so right hand side is transformed at each stage from one vector to another
• if we perform column operations, the column operations build from left to right, so we have to remember each transformation at each stage, and apply them all in reverse order at the end
• this means column operations are much less useful