Try our Free Online Math Solver!

Numerical Linear Algebra
3.1 – Overview
Numerical linear algebra includes a class of problems which involves multiplying
vectors and matrices, solving linear equations, computing the inverse of a
matrix, and
computing various matrix decompositions. It is one of the areas where knowledge
of
appropriate numerical methods can lead to rather substantial speed improvements.
3.2 – Types of Matrices
In numerical linear algebra, efficient algorithms will take advantage of the
special
structure of a matrix. In addition to the rectangular matrix, special matrix
types include
symmetric matrices, symmetric positive definite matrices, banddiagonal
matrices, blockdiagonal
matrices, lower and upper triangular matrices, lower and upper band diagonal
matrices, Toeplitz matrices, and sparse matrices.
Each of these matrices is stored in a different way ,
generally in an array in
memory. A general matrix is usually stored as one row after another in memory. A
symmetric matrix can be stored in either packed or unpacked storage, and we may
store
either the lower triangle, the upper triangle, or both. Suppose that n is the
dimension of a
symmetric matrix. Unpacked storage requires n^{2} locations in memory, while packed
storage requires n(n +1) / 2 units. In both cases, we typically store symmetric
matrices
row after row. Lower and upper triangular matrices may also be stored in both
packed
and unpacked form.
Figures 3.1 and 3.2 illustrate the storage of lower, upper, and symmetric
matrices
in packed and unpacked storage. Bold number illustrate a location in memory that
contains an element of the matrix and nonbold elements contain a location in
memory
which contains garbage. The advantage of packed storage is, of course, efficient
use of
memory. The advantage of unpacked storage is that this storage usually allows
for faster
computation. This occurs for a number of reasons, but one important reasons is
that block
versions of common algorithms can be applied only if matrices are stored in
unpacked
form.
Band diagonal matrices are stored in a different fashion. Let
denote the
number
of lower bands, let denote the number of upper bands, and let
. A band
diagonal matrix is stored in an nd element arrays. Each d units in the array
correspond
to the nonzero elements of a row of A . For example, consider the following
example of a
5x5 band diagonal matrix with and
. Figure 3.3 illustrates the
storage of
such a matrix. The advantage of storing the matrix this way is both efficient
use of
memory and computational speed.
Figure 3.1: Storage of Lower, Upper, and Symmetric
Matrices
(Unpacked Storage)
Figure 3.2: Storage of Lower, Upper, and Symmetric
Matrices
(Packed Storage)
Figure 3.3: Storage of Band Matrix
A sparse matrix is a large matrix that contains few
non zero elements . Sparse
matrices are most often stored in compressed row format. Compressed row format
requires 3 arrays, which we label r , c , and v . The array v contains the
values of all of
the nonzero elements in the matrix. The array c contains column in which each
corresponding nonzero element is located. The array r points to the location in
c which
contains the elements that correspond to a row.
For example, consider the following sparse matrix,
We could store this matrix as
The final type of matrix (which is often useful in
timeseries econometrics
applications) is a Toeplitz matrix. A Toeplitz matrix has the form,
A Toeplitz matrix can be stored as,
.
3.3 – Elementary Linear Algebra and the BLAS
Elementary linear algebra operations include operations
such as y← Ax ,
C ← A^{T} B , and y← L^{1}x . Often, these operations can be implemented in a few
lines of
code. There are, however, more efficient ways to implement these types of
algorithms.
This differences most important in matrixmatrix computation (e.g. C ← A^{T} B ),
but can
occasionally be useful in matrixvector computations (e.g. y← Ax , y← L^{1}x ).
The Basic Linear Algebra Subroutines (BLAS) contains efficient implementations
of these and many other algorithms. The BLAS are a set of FORTRAN subroutines,
but
CBLAS ports these subroutines to C using an automatic converter called
Fortran^{2}C. In
order to apply the CBLAS, you must be relatively familiar with pointers.
Consider
GEMM (which is the routine for multiplying two rectangular matrices). The
routine
computes C ←α op(A)op(B) +β C , where op(A) = A or op(A) = A^{T} .
Notice that every variable here is a pointer. The variable transa determines
where
op(A) = A or op(A) = A^{T} . The variables m and n are the dimensions of C . The
variable
k is the number of columns or A , which is also equal to the number of rows in B
. The
variables lda, ldb, and ldc, store the “leading dimension” of A , B , and C .
This option
allows the algorithm to refer to the following types of matrices,
In this case, we have a 4 by 5 matrix, with a leading
dimension of 8, starting at location 3.
Notice that this allows us the call the BLAS algorithms on sub matrices without
reallocating
memory.
BLAS is divided in to level 1, level2, and level3 algorithms. Level1
algorithms
perform operations on vectors (and are O(n) ), level2 algorithms operate on a
vector and
a matrix (and are O(n^{2} ) ), and level3 algorithms operate on two matrices (and
are
O(n^{3} ) ).
Prev  Next 