MatrixOps(3) User Contributed Perl Documentation MatrixOps(3)NAMEPDL::MatrixOps-- Some Useful Matrix Operations
SYNOPSIS
$inv = $a->inv;
$det = $a->det;
($lu,$perm,$par) = $a->lu_decomp;
$x = lu_backsub($lu,$perm,$b); # solve $a x $x = $b
DESCRIPTIONPDL::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 no FORTRAN compiler or
external library available (e.g. PDL::Slatec or PDL::GSL).
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.
TIPS ON MATRIX OPERATIONS
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)));
ACKNOWLEDGEMENTS
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. All are
free software, distributable under same terms as PDL itself.
NOTES
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!
FUNCTIONS
identity
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.
stretcher
Signature: (a(n); [o]b(n,n))
$mat = stretcher($eigenvalues);
Return a diagonal matrix with the specified diagonal elements
inv
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. If you ask it to thread then a
numerically unstable (non-pivoting) method is used instead, so avoid
threading over collections of large (say, more than 4x4) or near-
singular matrices unless precision is not important.
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. In the threading case, no checking for singularity is
performed, if any of the matrices in your threaded collection are
singular, they receive NaN entries.
· 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.
det
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.
If you ask det to thread (by giving it a 3-D or higher dim piddle) then
lu_decomp drops you through to lu_decomp2, which is numerically
unstable (and hence not useful for very large matrices) but quite fast.
If you want to use threading on a matrix that's less than, say, 10x10,
and might be near singular, then you might want to use determinant,
which is a more robust (but slower) determinant finder, instead.
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.
determinant
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, but does not thread well.
eigens_sym
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'vals
$e = eigens_sym($a); # just eigenvalues
eigens
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).
PLEASE NOTE: THE EIGENS FUNCTION HAS BEEN DISABLED FOR ASYMMETRIC
MATRICES PENDING FIX TO THE CODE.
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 (anyone know what algorithm is
used?). 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
svd
Signature: (a(n,m); [o]u(n,m); [o,phys]z(n); [o]v(n,n))
($r1, $s, $r2) = svd($a);
Singular value decomposition of a matrix.
"svd" is threadable.
$r1 and $r2 are rotation matrices that convert from the original
matrix's singular coordinates to final coordinates, and from original
coordinates to singular coordinates, respectively. $s is the diagonal
of the singular value matrix, so that, if $a is square, then you can
make an expensive copy of $a by saying:
$ess = zeroes($r1); $ess->diagonal(0,1) .= $s;
$a_copy .= $r2 x $ess x $r1;
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
}
lu_decomp
Signature: (a(m,m); [o]b(n); [o]c; [o]lu)
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 defeats efficient
threading, so if you have a few matrices to decompose accurately, you
should use lu_decomp, but if you have a million matrices to decompose
and don't mind a higher error budget you probably want to use
lu_decomp2, which doesn't do the pivoting (and hence gives wrong
answers for near-singular or large matrices), but does do 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 inverse.
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, but makes it unable to properly vectorize
threaded operation. If you have a LOT of small matrices to invert
(like, say, a 3x3x1000000 PDL) you should use lu_decomp2, which doesn't
pivot and is therefore threadable (and, of course, works in-place).
If you ask lu_decomp to thread (by having a nontrivial third dimension
in the matrix) then it will call lu_decomp2 instead. That is a
numerically unstable (non-pivoting) method that is mainly useful for
smallish, not-so-singular matrices but is threadable.
lu_decomp is ported from _Numerical_Recipes to PDL. It should probably
be implemented in C.
lu_decomp2
Signature: (a(m,m); [0]lu(n)
LU decompose a matrix, with no row permutation (threadable!)
($lu, $perm, $parity) = lu_decomp2($a);
$lu = lu_decomp2($a,[$perm,$par]);
lu_decomp($a->inplace,[$perm,$par]);
"lu_decomp2" works just like lu_decomp, but it does no pivoting at all
and hence can be usefully threaded. 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's just the ticket.
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. If lu_decomp
were implemented in C, then lu_decomp2 might become unnecessary.
lu_backsub
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) = lu_decomp($a);
$x = lu_backsub($lu,$perm,$par,$b);
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 inverse, and just multiply "X =
A^-1 B".) In fact, inverse 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.
simq
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.
squaretotri
Signature: (a(n,n); b(m))
Convert a symmetric square matrix to triangular vector storage.
AUTHOR
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.
perl v5.10.0 2008-08-29 MatrixOps(3)
