NAME

       PDL::MatrixOps -- Some Useful Matrix Operations

SYNOPSIS

           $inv = $x->inv;

           $det = $x->det;

           ($lu,$perm,$par) = $x->lu_decomp;
           $y = lu_backsub($lu,$perm,$z); # solve $x x $y = $z

DESCRIPTION

       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.

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, broadcastable):

           $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->slice('*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 broadcasting 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 broadcast over multiple row vectors.

       When broadcasting 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:

               $x = xvals(360)*3.14159/180;        # (angle)
               $rot = cat(cat(cos($x),sin($x)),    # rotmat: (col,row,angle)
                          cat(-sin($x),cos($x)));

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.  They 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 first two dimensions of the matrix, with higher dimensions preserved.

   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.

       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.

   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.

       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($x);

       Determinant of a square matrix, using recursive descent (broadcastable).

       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 broadcastable, 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 broadcasting benchmarks with new code).

   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->transpose).

       It's broadcastable, 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 $x;  # Make eigenvector matrix
           $vector = $ev->slice($n);       # Select nth eigenvector as a column-vector
           $vector = $ev->slice("($n)");     # Select nth eigenvector as a row-vector

           ($ev, $e) = eigens_sym($x); # e-vects & e-values
           $e = eigens_sym($x);        # just eigenvalues

       eigens_sym ignores the bad-value flag of the input ndarrays.  It will set the bad-value
       flag of all output ndarrays if the flag is set for any of the input ndarrays.

   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).

       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 broadcastable, 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 $x;  # Make eigenvector matrix
               $vector = $ev->slice($n);   # Select nth eigenvector as a column-vector
               $vector = $ev->slice("($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($x); # e'vects & e'vals
           $e = eigens($x);        # just eigenvalues

       eigens ignores the bad-value flag of the input ndarrays.  It will set the bad-value flag
       of all output ndarrays if the flag is set for any of the input ndarrays.

   svd
         Signature: (a(n,m); [t]w(wsize); [o]u(n,m); [o,phys]z(n); [o]v(n,n))

        ($u, $s, $v) = svd($x);

       Singular value decomposition of a matrix.

       "svd" is broadcastable.

       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($x);
           $ess = zeroes($x->dim(0),$x->dim(0));
           $ess->slice("$_","$_").=$s->slice("$_") foreach (0..$x->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($x->transpose);
           $ess = zeroes($x->dim(1),$x->dim(1));
           $ess->slice($_,$_).=$s->slice($_) foreach (0..$x->dim(1)-1); #generic diagonal
           $x_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 $x;
             $s++;             # fatten all singular values
             $r2 *= $s;        # implicit broadcasting for cheap mult.
             $x .= $r2 x $r1;  # a gets r2 x ess x r1
           }

       svd ignores the bad-value flag of the input ndarrays.  It will set the bad-value flag of
       all output ndarrays if the flag is set for any of the input ndarrays.

   lu_decomp
         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($x);

         $lu = lu_decomp($x, $perm, $par);  # $perm and $par are outputs!

         lu_decomp($x->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 broadcasting.

       "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 broadcasted, 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.

   lu_decomp2
         Signature: (a(m,m); [o]lu(m,m))

       LU decompose a matrix, with no row permutation

         ($lu, $perm, $parity) = lu_decomp2($x);

         $lu = lu_decomp2($x,$perm,$parity);   # or
         $lu = lu_decomp2($x);                 # $perm and $parity are optional

         lu_decomp($x->inplace,$perm,$parity); # or
         lu_decomp($x->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.

   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,$par) = lu_decomp($A);

         $x = lu_backsub($lu,$perm,$par,$A);  # 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)

         # starting from square matrix A and columns matrix B, with mathematically
         # correct dimensions
         $A = identity(4) + ones(4, 4);
         $A->slice('2,0') .= 0; # break symmetry to see if need transpose
         $B = sequence(2, 4);
         # all these functions take B as rows, interpret as though notional columns
         # mathematically confusing but can't change as back-compat and also
         # familiar to Fortran users, so just transpose inputs and outputs

         # using lu_backsub
         ($lu,$perm,$par) = lu_decomp($A);
         $x = lu_backsub($lu,$perm,$par, $B->transpose)->transpose;

         # or with Slatec LINPACK
         use PDL::Slatec;
         gefa($lu=$A->copy, $ipiv=null, $info=null);
         # 1 = do transpose because Fortran's idea of rows vs columns
         gesl($lu, $ipiv, $x=$B->transpose->copy, 1);
         $x = $x->inplace->transpose;

         # or with LAPACK
         use PDL::LinearAlgebra::Real;
         getrf($lu=$A->copy, $ipiv=null, $info=null);
         getrs($lu, 1, $x=$B->transpose->copy, $ipiv, $info=null); # again, need transpose
         $x=$x->inplace->transpose;

         # or with GSL
         use PDL::GSL::LINALG;
         LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null);
         # $B and $x, first dim is because GSL treats as vector, higher dims broadcast
         # so we transpose in and back out
         LU_solve($lu, $p, $B->transpose, my $x=null);
         $x=$x->inplace->transpose;

         # proof of the pudding is in the eating:
         print $A x $x;

       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 broadcasts in the usual way, which means that it
       treats $y 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.

   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).

       $y, $x, $ips are vectors of length "n".

       Set "flag=0" to solve.  Set "flag=-1" to do a new back substitution for different $y
       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 ndarrays.  It will set the bad-value flag of
       all output ndarrays if the flag is set for any of the input ndarrays.

   squaretotri
         Signature: (a(n,n); [o]b(m))

       Convert a lower-triangular square matrix to triangular vector storage.  Ignores upper half
       of input.

       squaretotri does not process bad values.  It will set the bad-value flag of all output
       ndarrays if the flag is set for any of the input ndarrays.

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.