|
/usr/local/perl/lib/site_perl/5.6.1/sun4-solaris/PDL/MatrixOps.pm |
PDL::MatrixOps
$inv = $a->inv;
$det = $a->det;
($lu,$perm,$par) = $a->lu_decomp; $x = lu_backsub($lu,$perm,$b); # solve $a x $x = $b
PDL::Math::MatrixOps contains a bunch of operations for handling matrices -- LU decomposition, inversion, determinant, etc. 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.
It 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. (Contrast this with, e.g., IDL's treatment of matrices -- which appear like their own transposes in that language). If you prefer your matrices indexed in (column, row) order, you can try using the PDL::Matrix object, which includes an implicit xchange 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:
$col_vec2 = $m1 x $mypdl->(0,1); # mypdl used as a 1xn column vector
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 many 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)));
This is intended as a general-purpose linear algebra package. If there is something you want that is not here, please add and document it!
use Carp; use PDL::NiceSlice; use strict;
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. If you ask it to thread then a numerically
unstable (non-pivoting) method is used instead, so avoid thread ing
over largish or near-singular matrices unless precision is not important.
OPTIONS:
$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.
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:
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.
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.
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).
($ev, $e) = eigens($a); # e'vects & e'vals $e = eigens($a); # just eigenvalues
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
}
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 the lu_decomp2 manpage, 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.
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 the lu_decomp manpage, 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 the lu_decomp manpage 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.
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.
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.
Signature: (a(n,n); b(m))
Convert a symmetric square matrix to triangular vector storage.
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.
|
/usr/local/perl/lib/site_perl/5.6.1/sun4-solaris/PDL/MatrixOps.pm |