Skip site navigation (1)Skip section navigation (2)

FreeBSD Manual Pages


home | help
MatrixOps(3)	      User Contributed Perl Documentation	  MatrixOps(3)

       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)

       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

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

       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.


       o  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

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

       o  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

       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.


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

       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


       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;


       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

       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

	 ($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

       "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

       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

       Copyright (C) 2002 Craig	DeForest (, R.J.R.
       Williams	(,	Karl Glazebrook
       (	 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.32.1			  2021-02-28			  MatrixOps(3)


Want to link to this manual page? Use this URL:

home | help