PDL::MatrixOps -- Some Useful Matrix Operations
$inv = $a->inv;
$det = $a->det;
($lu,$perm,$par) = $a->lu_decomp;
$x = lu_backsub($lu,$perm,$b); # solve $a x $x = $b
PDL::MatrixOps is PDL's built-in matrix manipulation code. It contains utilities
for many common matrix operations: inversion, determinant finding,
eigenvalue/vector finding, singular value decomposition, etc. PDL::MatrixOps
routines are written in a mixture of Perl and C, so that they are reliably
present even when there is no FORTRAN compiler or external library available
(e.g. PDL::Slatec or any of the PDL::GSL family of modules).
Matrix manipulation, particularly with large matrices, is a challenging field
and no one algorithm is suitable in all cases. The utilities here use
general-purpose algorithms that work acceptably for many cases but might not
scale well to very large or pathological (near-singular) matrices.
Except as noted, the matrices are PDLs whose 0th dimension ranges over column
and whose 1st dimension ranges over row. The matrices appear correctly when
printed.
These routines should work OK with PDL::Matrix objects as well as with normal
PDLs.
Like most computer languages, PDL addresses matrices in (column,row) order in
most cases; this corresponds to (X,Y) coordinates in the matrix itself,
counting rightwards and downwards from the upper left corner. This means that
if you print a PDL that contains a matrix, the matrix appears correctly on the
screen, but if you index a matrix element, you use the indices in the reverse
order that you would in a math textbook. If you prefer your matrices indexed
in (row, column) order, you can try using the PDL::Matrix object, which
includes an implicit exchange of the first two dimensions but should be
compatible with most of these matrix operations. TIMTOWDTI.)
Matrices, row vectors, and column vectors can be multiplied with the 'x'
operator (which is, of course, threadable):
$m3 = $m1 x $m2;
$col_vec2 = $m1 x $col_vec1;
$row_vec2 = $row_vec1 x $m1;
$scalar = $row_vec x $col_vec;
Because of the (column,row) addressing order, 1-D PDLs are treated as _row_
vectors; if you want a _column_ vector you must add a dummy dimension:
$rowvec = pdl(1,2); # row vector
$colvec = $rowvec->(*1); # 1x2 column vector
$matrix = pdl([[3,4],[6,2]]); # 2x2 matrix
$rowvec2 = $rowvec x $matrix; # right-multiplication by matrix
$colvec = $matrix x $colvec; # left-multiplication by matrix
$m2 = $matrix x $rowvec; # Throws an error
Implicit threading works correctly with most matrix operations, but you must be
extra careful that you understand the dimensionality. In particular, matrix
multiplication and other matrix ops need nx1 PDLs as row vectors and 1xn PDLs
as column vectors. In most cases you must explicitly include the trailing 'x1'
dimension in order to get the expected results when you thread over multiple
row vectors.
When threading over matrices, it's very easy to get confused about which
dimension goes where. It is useful to include comments with every expression,
explaining what you think each dimension means:
$a = xvals(360)*3.14159/180; # (angle)
$rot = cat(cat(cos($a),sin($a)), # rotmat: (col,row,angle)
cat(-sin($a),cos($a)));
MatrixOps includes algorithms and pre-existing code from several origins. In
particular, "eigens_sym" is the work of Stephen Moshier,
"svd" uses an SVD subroutine written by Bryant Marks, and
"eigens" uses a subset of the Small Scientific Library by Kenneth
Geisshirt. They are free software, distributable under same terms as PDL
itself.
This is intended as a general-purpose linear algebra package for small-to-mid
sized matrices. The algorithms may not scale well to large matrices (hundreds
by hundreds) or to near singular matrices.
If there is something you want that is not here, please add and document it!
Signature: (n; [o]a(n,n))
Return an identity matrix of the specified size. If you hand in a scalar, its
value is the size of the identity matrix; if you hand in a dimensioned PDL,
the 0th dimension is the size of the matrix.
Signature: (a(n); [o]b(n,n))
$mat = stretcher($eigenvalues);
Return a diagonal matrix with the specified diagonal elements
Signature: (a(m,m); sv opt )
$a1 = inv($a, {$opt});
Invert a square matrix.
You feed in an NxN matrix in $a, and get back its inverse (if it exists). The
code is inplace-aware, so you can get back the inverse in $a itself if you
want -- though temporary storage is used either way. You can cache the LU
decomposition in an output option variable.
"inv" uses "lu_decomp" by default; that is a numerically
stable (pivoting) LU decomposition method.
OPTIONS:
- •
- s
Boolean value indicating whether to complain if the matrix is singular. If
this is false, singular matrices cause inverse to barf. If it is true,
then singular matrices cause inverse to return undef.
- •
- lu (I/O)
This value contains a list ref with the LU decomposition, permutation, and
parity values for $a. If you do not mention the key, or if the value is
undef, then inverse calls "lu_decomp". If the key exists with an
undef value, then the output of "lu_decomp" is stashed here
(unless the matrix is singular). If the value exists, then it is assumed
to hold the LU decomposition.
- •
- det (Output)
If this key exists, then the determinant of $a get stored here, whether or
not the matrix is singular.
Signature: (a(m,m); sv opt)
$det = det($a,{opt});
Determinant of a square matrix using LU decomposition (for large matrices)
You feed in a square matrix, you get back the determinant. Some options exist
that allow you to cache the LU decomposition of the matrix (note that the LU
decomposition is invalid if the determinant is zero!). The LU decomposition is
cacheable, in case you want to re-use it. This method of determinant finding
is more rapid than recursive-descent on large matrices, and if you reuse the
LU decomposition it's essentially free.
OPTIONS:
- •
- lu (I/O)
Provides a cache for the LU decomposition of the matrix. If you provide the
key but leave the value undefined, then the LU decomposition goes in here;
if you put an LU decomposition here, it will be used and the matrix will
not be decomposed again.
Signature: (a(m,m))
$det = determinant($a);
Determinant of a square matrix, using recursive descent (threadable).
This is the traditional, robust recursive determinant method taught in most
linear algebra courses. It scales like "O(n!)" (and hence is
pitifully slow for large matrices) but is very robust because no division is
involved (hence no division-by-zero errors for singular matrices). It's also
threadable, so you can find the determinants of a large collection of matrices
all at once if you want.
Matrices up to 3x3 are handled by direct multiplication; larger matrices are
handled by recursive descent to the 3x3 case.
The LU-decomposition method det is faster in isolation for single matrices
larger than about 4x4, and is much faster if you end up reusing the LU
decomposition of $a (NOTE: check performance and threading benchmarks with new
code).
Signature: ([phys]a(m); [o,phys]ev(n,n); [o,phys]e(n))
Eigenvalues and -vectors of a symmetric square matrix. If passed an asymmetric
matrix, the routine will warn and symmetrize it, by taking the average value.
That is, it will solve for 0.5*($a+$a->mv(0,1)).
It's threadable, so if $a is 3x3x100, it's treated as 100 separate 3x3 matrices,
and both $ev and $e get extra dimensions accordingly.
If called in scalar context it hands back only the eigenvalues. Ultimately, it
should switch to a faster algorithm in this case (as discarding the
eigenvectors is wasteful).
The algorithm used is due to J. vonNeumann, which was a rediscovery of Jacobi's
Method <http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm> .
The eigenvectors are returned in COLUMNS of the returned PDL. That makes it
slightly easier to access individual eigenvectors, since the 0th dim of the
output PDL runs across the eigenvectors and the 1st dim runs across their
components.
($ev,$e) = eigens_sym $a; # Make eigenvector matrix
$vector = $ev->($n); # Select nth eigenvector as a column-vector
$vector = $ev->(($n)); # Select nth eigenvector as a row-vector
($ev, $e) = eigens_sym($a); # e-vects & e-values
$e = eigens_sym($a); # just eigenvalues
eigens_sym ignores the bad-value flag of the input piddles. It will set the
bad-value flag of all output piddles if the flag is set for any of the input
piddles.
Signature: ([phys]a(m); [o,phys]ev(l,n,n); [o,phys]e(l,n))
Real eigenvalues and -vectors of a real square matrix.
(See also "eigens_sym", for eigenvalues and -vectors of a real,
symmetric, square matrix).
The eigens function will attempt to compute the eigenvalues and eigenvectors of
a square matrix with real components. If the matrix is symmetric, the same
underlying code as "eigens_sym" is used. If asymmetric, the
eigenvalues and eigenvectors are computed with algorithms from the sslib
library. If any imaginary components exist in the eigenvalues, the results are
currently considered to be invalid, and such eigenvalues are returned as
"NaN"s. This is true for eigenvectors also. That is if there are
imaginary components to any of the values in the eigenvector, the eigenvalue
and corresponding eigenvectors are all set to "NaN". Finally, if
there are any repeated eigenvectors, they are replaced with all
"NaN"s.
Use of the eigens function on asymmetric matrices should be considered
experimental! For asymmetric matrices, nearly all observed matrices with real
eigenvalues produce incorrect results, due to errors of the sslib algorithm.
If your assymmetric matrix returns all NaNs, do not assume that the values are
complex. Also, problems with memory access is known in this library.
Not all square matrices are diagonalizable. If you feed in a non-diagonalizable
matrix, then one or more of the eigenvectors will be set to NaN, along with
the corresponding eigenvalues.
"eigens" is threadable, so you can solve 100 eigenproblems by feeding
in a 3x3x100 array. Both $ev and $e get extra dimensions accordingly.
If called in scalar context "eigens" hands back only the eigenvalues.
This is somewhat wasteful, as it calculates the eigenvectors anyway.
The eigenvectors are returned in COLUMNS of the returned PDL (ie the the 0
dimension). That makes it slightly easier to access individual eigenvectors,
since the 0th dim of the output PDL runs across the eigenvectors and the 1st
dim runs across their components.
($ev,$e) = eigens $a; # Make eigenvector matrix
$vector = $ev->($n); # Select nth eigenvector as a column-vector
$vector = $ev->(($n)); # Select nth eigenvector as a row-vector
DEVEL NOTES:
For now, there is no distinction between a complex eigenvalue and an invalid
eigenvalue, although the underlying code generates complex numbers. It might
be useful to be able to return complex eigenvalues.
($ev, $e) = eigens($a); # e'vects & e'vals
$e = eigens($a); # just eigenvalues
eigens ignores the bad-value flag of the input piddles. It will set the
bad-value flag of all output piddles if the flag is set for any of the input
piddles.
Signature: (a(n,m); [o]u(n,m); [o,phys]z(n); [o]v(n,n))
($u, $s, $v) = svd($a);
Singular value decomposition of a matrix.
"svd" is threadable.
Given an m x n matrix $a that has m rows and n columns (m >= n),
"svd" computes matrices $u and $v, and a vector of the singular
values $s. Like most implementations, "svd" computes what is
commonly referred to as the "thin SVD" of $a, such that $u is m x n,
$v is n x n, and there are <=n singular values. As long as m >= n, the
original matrix can be reconstructed as follows:
($u,$s,$v) = svd($a);
$ess = zeroes($a->dim(0),$a->dim(0));
$ess->slice("$_","$_").=$s->slice("$_") foreach (0..$a->dim(0)-1); #generic diagonal
$a_copy = $u x $ess x $v->transpose;
If m==n, $u and $v can be thought of as rotation matrices that convert from the
original matrix's singular coordinates to final coordinates, and from original
coordinates to singular coordinates, respectively, and $ess is a diagonal
scaling matrix.
If n>m, "svd" will barf. This can be avoided by passing in the
transpose of $a, and reconstructing the original matrix like so:
($u,$s,$v) = svd($a->transpose);
$ess = zeroes($a->dim(1),$a->dim(1));
$ess->slice("$_","$_").=$s->slice("$_") foreach (0..$a->dim(1)-1); #generic diagonal
$a_copy = $v x $ess x $u->transpose;
EXAMPLE
The computing literature has loads of examples of how to use SVD. Here's a
trivial example (used in PDL::Transform::map) of how to make a matrix less,
er, singular, without changing the orientation of the ellipsoid of
transformation:
{ my($r1,$s,$r2) = svd $a;
$s++; # fatten all singular values
$r2 *= $s; # implicit threading for cheap mult.
$a .= $r2 x $r1; # a gets r2 x ess x r1
}
svd ignores the bad-value flag of the input piddles. It will set the bad-value
flag of all output piddles if the flag is set for any of the input piddles.
Signature: (a(m,m); [o]lu(m,m); [o]perm(m); [o]parity)
LU decompose a matrix, with row permutation
($lu, $perm, $parity) = lu_decomp($a);
$lu = lu_decomp($a, $perm, $par); # $perm and $par are outputs!
lu_decomp($a->inplace,$perm,$par); # Everything in place.
"lu_decomp" returns an LU decomposition of a square matrix, using
Crout's method with partial pivoting. It's ported from
Numerical
Recipes. The partial pivoting keeps it numerically stable but means a
little more overhead from threading.
"lu_decomp" decomposes the input matrix into matrices L and U such
that LU = A, L is a subdiagonal matrix, and U is a superdiagonal matrix. By
convention, the diagonal of L is all 1's.
The single output matrix contains all the variable elements of both the L and U
matrices, stacked together. Because the method uses pivoting (rearranging the
lower part of the matrix for better numerical stability), you have to permute
input vectors before applying the L and U matrices. The permutation is
returned either in the second argument or, in list context, as the second
element of the list. You need the permutation for the output to make any
sense, so be sure to get it one way or the other.
LU decomposition is the answer to a lot of matrix questions, including inversion
and determinant-finding, and "lu_decomp" is used by inv.
If you pass in $perm and $parity, they either must be predeclared PDLs of the
correct size ($perm is an n-vector, $parity is a scalar) or scalars.
If the matrix is singular, then the LU decomposition might not be defined; in
those cases, "lu_decomp" silently returns undef. Some singular
matrices LU-decompose just fine, and those are handled OK but give a zero
determinant (and hence can't be inverted).
"lu_decomp" uses pivoting, which rearranges the values in the matrix
for more numerical stability. This makes it really good for large and even
near-singular matrices. There is a non-pivoting version "lu_decomp2"
available which is from 5 to 60 percent faster for typical problems at the
expense of failing to compute a result in some cases.
Now that the "lu_decomp" is threaded, it is the recommended LU
decomposition routine. It no longer falls back to "lu_decomp2".
"lu_decomp" is ported from
Numerical Recipes to PDL. It should
probably be implemented in C.
Signature: (a(m,m); [o]lu(m,m))
LU decompose a matrix, with no row permutation
($lu, $perm, $parity) = lu_decomp2($a);
$lu = lu_decomp2($a,$perm,$parity); # or
$lu = lu_decomp2($a); # $perm and $parity are optional
lu_decomp($a->inplace,$perm,$parity); # or
lu_decomp($a->inplace); # $perm and $parity are optional
"lu_decomp2" works just like lu_decomp, but it does
no pivoting
at all. For compatibility with lu_decomp, it will give you a permutation list
and a parity scalar if you ask for them -- but they are always trivial.
Because "lu_decomp2" does not pivot, it is numerically
unstable
-- that means it is less precise than lu_decomp, particularly for large or
near-singular matrices. There are also specific types of non-singular matrices
that confuse it (e.g. ([0,-1,0],[1,0,0],[0,0,1]), which is a 90 degree
rotation matrix but which confuses "lu_decomp2").
On the other hand, if you want to invert rapidly a few hundred thousand small
matrices and don't mind missing one or two, it could be the ticket. It can be
up to 60% faster at the expense of possible failure of the decomposition for
some of the input matrices.
The output is a single matrix that contains the LU decomposition of $a; you can
even do it in-place, thereby destroying $a, if you want. See lu_decomp for
more information about LU decomposition.
"lu_decomp2" is ported from
Numerical Recipes into PDL.
Signature: (lu(m,m); perm(m); b(m))
Solve a x = b for matrix a, by back substitution into a's LU decomposition.
($lu,$perm,$par) = lu_decomp($a);
$x = lu_backsub($lu,$perm,$par,$b); # or
$x = lu_backsub($lu,$perm,$b); # $par is not required for lu_backsub
lu_backsub($lu,$perm,$b->inplace); # modify $b in-place
$x = lu_backsub(lu_decomp($a),$b); # (ignores parity value from lu_decomp)
Given the LU decomposition of a square matrix (from lu_decomp),
"lu_backsub" does back substitution into the matrix to solve "a
x = b" for given vector "b". It is separated from the
"lu_decomp" method so that you can call the cheap
"lu_backsub" multiple times and not have to do the expensive LU
decomposition more than once.
"lu_backsub" acts on single vectors and threads in the usual way,
which means that it treats $b as the
transpose of the input. If you
want to process a matrix, you must hand in the
transpose of the matrix,
and then transpose the output when you get it back. that is because pdls are
indexed by (col,row), and matrices are (row,column) by convention, so a 1-D
pdl corresponds to a row vector, not a column vector.
If $lu is dense and you have more than a few points to solve for, it is probably
cheaper to find "a^-1" with inv, and just multiply "x = a^-1
b".) in fact, inv works by calling "lu_backsub" with the
identity matrix.
"lu_backsub" is ported from section 2.3 of
Numerical Recipes.
It is written in PDL but should probably be implemented in C.
Signature: ([phys]a(n,n); [phys]b(n); [o,phys]x(n); int [o,phys]ips(n); int flag)
Solution of simultaneous linear equations, "a x = b".
$a is an "n x n" matrix (i.e., a vector of length "n*n"),
stored row-wise: that is, "a(i,j) = a[ij]", where "ij = i*n +
j".
While this is the transpose of the normal column-wise storage, this corresponds
to normal PDL usage. The contents of matrix a may be altered (but may be
required for subsequent calls with flag = -1).
$b, $x, $ips are vectors of length "n".
Set "flag=0" to solve. Set "flag=-1" to do a new back
substitution for different $b vector using the same a matrix previously
reduced when "flag=0" (the $ips vector generated in the previous
solution is also required).
See also lu_backsub, which does the same thing with a slightly less opaque
interface.
simq ignores the bad-value flag of the input piddles. It will set the bad-value
flag of all output piddles if the flag is set for any of the input piddles.
Signature: (a(n,n); b(m))
Convert a symmetric square matrix to triangular vector storage.
squaretotri does not process bad values. It will set the bad-value flag of all
output piddles if the flag is set for any of the input piddles.
Copyright (C) 2002 Craig DeForest (deforest@boulder.swri.edu), R.J.R. Williams
(rjrw@ast.leeds.ac.uk), Karl Glazebrook (kgb@aaoepp.aao.gov.au). There is no
warranty. You are allowed to redistribute and/or modify this work under the
same conditions as PDL itself. If this file is separated from the PDL
distribution, then the PDL copyright notice should be included in this
file.