1MatrixOps(3) User Contributed Perl Documentation MatrixOps(3)
2
6 PDL::MatrixOps -- Some Useful Matrix Operations
7
9 $inv = $x->inv;
10 $det = $x->det;
12 ($lu,$perm,$par) = $x->lu_decomp;
14 $y = lu_backsub($lu,$perm,$z); # solve $x x $y = $z
15
17 PDL::MatrixOps is PDL's built-in matrix manipulation code. It contains
18 utilities for many common matrix operations: inversion, determinant
19 finding, eigenvalue/vector finding, singular value decomposition, etc.
20 PDL::MatrixOps routines are written in a mixture of Perl and C, so that
21 they are reliably present even when there is no FORTRAN compiler or
22 external library available (e.g. PDL::Slatec or any of the PDL::GSL
23 family of modules).
24 Matrix manipulation, particularly with large matrices, is a challenging
26 field and no one algorithm is suitable in all cases. The utilities
27 here use general-purpose algorithms that work acceptably for many cases
28 but might not scale well to very large or pathological (near-singular)
29 matrices.
30 Except as noted, the matrices are PDLs whose 0th dimension ranges over
32 column and whose 1st dimension ranges over row. The matrices appear
33 correctly when printed.
34 These routines should work OK with PDL::Matrix objects as well as with
36 normal PDLs.
37
39 Like most computer languages, PDL addresses matrices in (column,row)
40 order in most cases; this corresponds to (X,Y) coordinates in the
41 matrix itself, counting rightwards and downwards from the upper left
42 corner. This means that if you print a PDL that contains a matrix, the
43 matrix appears correctly on the screen, but if you index a matrix
44 element, you use the indices in the reverse order that you would in a
45 math textbook. If you prefer your matrices indexed in (row, column)
46 order, you can try using the PDL::Matrix object, which includes an
47 implicit exchange of the first two dimensions but should be compatible
48 with most of these matrix operations. TIMTOWDTI.)
49 Matrices, row vectors, and column vectors can be multiplied with the
51 'x' operator (which is, of course, broadcastable):
52 $m3 = $m1 x $m2;
54 $col_vec2 = $m1 x $col_vec1;
55 $row_vec2 = $row_vec1 x $m1;
56 $scalar = $row_vec x $col_vec;
57 Because of the (column,row) addressing order, 1-D PDLs are treated as
59 _row_ vectors; if you want a _column_ vector you must add a dummy
60 dimension:
61 $rowvec = pdl(1,2); # row vector
63 $colvec = $rowvec->slice('*1'); # 1x2 column vector
64 $matrix = pdl([[3,4],[6,2]]); # 2x2 matrix
65 $rowvec2 = $rowvec x $matrix; # right-multiplication by matrix
66 $colvec = $matrix x $colvec; # left-multiplication by matrix
67 $m2 = $matrix x $rowvec; # Throws an error
68 Implicit broadcasting works correctly with most matrix operations, but
70 you must be extra careful that you understand the dimensionality. In
71 particular, matrix multiplication and other matrix ops need nx1 PDLs as
72 row vectors and 1xn PDLs as column vectors. In most cases you must
73 explicitly include the trailing 'x1' dimension in order to get the
74 expected results when you broadcast over multiple row vectors.
75 When broadcasting over matrices, it's very easy to get confused about
77 which dimension goes where. It is useful to include comments with every
78 expression, explaining what you think each dimension means:
79 $x = xvals(360)*3.14159/180; # (angle)
81 $rot = cat(cat(cos($x),sin($x)), # rotmat: (col,row,angle)
82 cat(-sin($x),cos($x)));
83
85 MatrixOps includes algorithms and pre-existing code from several
86 origins. In particular, "eigens_sym" is the work of Stephen Moshier,
87 "svd" uses an SVD subroutine written by Bryant Marks, and "eigens" uses
88 a subset of the Small Scientific Library by Kenneth Geisshirt. They
89 are free software, distributable under same terms as PDL itself.
90
92 This is intended as a general-purpose linear algebra package for small-
93 to-mid sized matrices. The algorithms may not scale well to large
94 matrices (hundreds by hundreds) or to near singular matrices.
95 If there is something you want that is not here, please add and
97 document it!
98
100 identity
101 Signature: (n; [o]a(n,n))
102 Return an identity matrix of the specified size. If you hand in a
104 scalar, its value is the size of the identity matrix; if you hand in a
105 dimensioned PDL, the 0th dimension is the first two dimensions of the
106 matrix, with higher dimensions preserved.
107 stretcher
109 Signature: (a(n); [o]b(n,n))
110 $mat = stretcher($eigenvalues);
112 Return a diagonal matrix with the specified diagonal elements
114 inv
116 Signature: (a(m,m); sv opt )
117 $a1 = inv($a, {$opt});
119 Invert a square matrix.
121 You feed in an NxN matrix in $a, and get back its inverse (if it
123 exists). The code is inplace-aware, so you can get back the inverse in
124 $a itself if you want -- though temporary storage is used either way.
125 You can cache the LU decomposition in an output option variable.
126 "inv" uses "lu_decomp" by default; that is a numerically stable
128 (pivoting) LU decomposition method.
129 OPTIONS:
131 • s
133 Boolean value indicating whether to complain if the matrix is
135 singular. If this is false, singular matrices cause inverse to
136 barf. If it is true, then singular matrices cause inverse to return
137 undef.
138 • lu (I/O)
140 This value contains a list ref with the LU decomposition,
142 permutation, and parity values for $a. If you do not mention the
143 key, or if the value is undef, then inverse calls "lu_decomp". If
144 the key exists with an undef value, then the output of "lu_decomp"
145 is stashed here (unless the matrix is singular). If the value
146 exists, then it is assumed to hold the LU decomposition.
147 • det (Output)
149 If this key exists, then the determinant of $a get stored here,
151 whether or not the matrix is singular.
152 det
154 Signature: (a(m,m); sv opt)
155 $det = det($a,{opt});
157 Determinant of a square matrix using LU decomposition (for large
159 matrices)
160 You feed in a square matrix, you get back the determinant. Some
162 options exist that allow you to cache the LU decomposition of the
163 matrix (note that the LU decomposition is invalid if the determinant is
164 zero!). The LU decomposition is cacheable, in case you want to re-use
165 it. This method of determinant finding is more rapid than recursive-
166 descent on large matrices, and if you reuse the LU decomposition it's
167 essentially free.
168 OPTIONS:
170 • lu (I/O)
172 Provides a cache for the LU decomposition of the matrix. If you
174 provide the key but leave the value undefined, then the LU
175 decomposition goes in here; if you put an LU decomposition here, it
176 will be used and the matrix will not be decomposed again.
177 determinant
179 Signature: (a(m,m))
180 $det = determinant($x);
182 Determinant of a square matrix, using recursive descent
184 (broadcastable).
185 This is the traditional, robust recursive determinant method taught in
187 most linear algebra courses. It scales like "O(n!)" (and hence is
188 pitifully slow for large matrices) but is very robust because no
189 division is involved (hence no division-by-zero errors for singular
190 matrices). It's also broadcastable, so you can find the determinants
191 of a large collection of matrices all at once if you want.
192 Matrices up to 3x3 are handled by direct multiplication; larger
194 matrices are handled by recursive descent to the 3x3 case.
195 The LU-decomposition method "det" is faster in isolation for single
197 matrices larger than about 4x4, and is much faster if you end up
198 reusing the LU decomposition of $a (NOTE: check performance and
199 broadcasting benchmarks with new code).
200 eigens_sym
202 Signature: ([phys]a(m); [o,phys]ev(n,n); [o,phys]e(n))
203 Eigenvalues and -vectors of a symmetric square matrix. If passed an
205 asymmetric matrix, the routine will warn and symmetrize it, by taking
206 the average value. That is, it will solve for 0.5*($a+$a->transpose).
207 It's broadcastable, so if $a is 3x3x100, it's treated as 100 separate
209 3x3 matrices, and both $ev and $e get extra dimensions accordingly.
210 If called in scalar context it hands back only the eigenvalues.
212 Ultimately, it should switch to a faster algorithm in this case (as
213 discarding the eigenvectors is wasteful).
214 The algorithm used is due to J. vonNeumann, which was a rediscovery of
216 Jacobi's Method
217 <http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm> .
218 The eigenvectors are returned in COLUMNS of the returned PDL. That
220 makes it slightly easier to access individual eigenvectors, since the
221 0th dim of the output PDL runs across the eigenvectors and the 1st dim
222 runs across their components.
223 ($ev,$e) = eigens_sym $x; # Make eigenvector matrix
225 $vector = $ev->slice($n); # Select nth eigenvector as a column-vector
226 $vector = $ev->slice("($n)"); # Select nth eigenvector as a row-vector
227 ($ev, $e) = eigens_sym($x); # e-vects & e-values
229 $e = eigens_sym($x); # just eigenvalues
230 eigens_sym ignores the bad-value flag of the input ndarrays. It will
232 set the bad-value flag of all output ndarrays if the flag is set for
233 any of the input ndarrays.
234 eigens
236 Signature: ([phys]a(m); [o,phys]ev(l,n,n); [o,phys]e(l,n))
237 Real eigenvalues and -vectors of a real square matrix.
239 (See also "eigens_sym", for eigenvalues and -vectors of a real,
241 symmetric, square matrix).
242 The eigens function will attempt to compute the eigenvalues and
244 eigenvectors of a square matrix with real components. If the matrix is
245 symmetric, the same underlying code as "eigens_sym" is used. If
246 asymmetric, the eigenvalues and eigenvectors are computed with
247 algorithms from the sslib library. If any imaginary components exist
248 in the eigenvalues, the results are currently considered to be invalid,
249 and such eigenvalues are returned as "NaN"s. This is true for
250 eigenvectors also. That is if there are imaginary components to any of
251 the values in the eigenvector, the eigenvalue and corresponding
252 eigenvectors are all set to "NaN". Finally, if there are any repeated
253 eigenvectors, they are replaced with all "NaN"s.
254 Use of the eigens function on asymmetric matrices should be considered
256 experimental! For asymmetric matrices, nearly all observed matrices
257 with real eigenvalues produce incorrect results, due to errors of the
258 sslib algorithm. If your assymmetric matrix returns all NaNs, do not
259 assume that the values are complex. Also, problems with memory access
260 is known in this library.
261 Not all square matrices are diagonalizable. If you feed in a non-
263 diagonalizable matrix, then one or more of the eigenvectors will be set
264 to NaN, along with the corresponding eigenvalues.
265 "eigens" is broadcastable, so you can solve 100 eigenproblems by
267 feeding in a 3x3x100 array. Both $ev and $e get extra dimensions
268 accordingly.
269 If called in scalar context "eigens" hands back only the eigenvalues.
271 This is somewhat wasteful, as it calculates the eigenvectors anyway.
272 The eigenvectors are returned in COLUMNS of the returned PDL (ie the
274 the 0 dimension). That makes it slightly easier to access individual
275 eigenvectors, since the 0th dim of the output PDL runs across the
276 eigenvectors and the 1st dim runs across their components.
277 ($ev,$e) = eigens $x; # Make eigenvector matrix
279 $vector = $ev->slice($n); # Select nth eigenvector as a column-vector
280 $vector = $ev->slice("($n)"); # Select nth eigenvector as a row-vector
281 DEVEL NOTES:
283 For now, there is no distinction between a complex eigenvalue and an
285 invalid eigenvalue, although the underlying code generates complex
286 numbers. It might be useful to be able to return complex eigenvalues.
287 ($ev, $e) = eigens($x); # e'vects & e'vals
289 $e = eigens($x); # just eigenvalues
290 eigens ignores the bad-value flag of the input ndarrays. It will set
292 the bad-value flag of all output ndarrays if the flag is set for any of
293 the input ndarrays.
294 svd
296 Signature: (a(n,m); [o]u(n,m); [o,phys]z(n); [o]v(n,n))
297 ($u, $s, $v) = svd($x);
299 Singular value decomposition of a matrix.
301 "svd" is broadcastable.
303 Given an m x n matrix $a that has m rows and n columns (m >= n), "svd"
305 computes matrices $u and $v, and a vector of the singular values $s.
306 Like most implementations, "svd" computes what is commonly referred to
307 as the "thin SVD" of $a, such that $u is m x n, $v is n x n, and there
308 are <=n singular values. As long as m >= n, the original matrix can be
309 reconstructed as follows:
310 ($u,$s,$v) = svd($x);
312 $ess = zeroes($x->dim(0),$x->dim(0));
313 $ess->slice("$_","$_").=$s->slice("$_") foreach (0..$x->dim(0)-1); #generic diagonal
314 $a_copy = $u x $ess x $v->transpose;
315 If m==n, $u and $v can be thought of as rotation matrices that convert
317 from the original matrix's singular coordinates to final coordinates,
318 and from original coordinates to singular coordinates, respectively,
319 and $ess is a diagonal scaling matrix.
320 If n>m, "svd" will barf. This can be avoided by passing in the
322 transpose of $a, and reconstructing the original matrix like so:
323 ($u,$s,$v) = svd($x->transpose);
325 $ess = zeroes($x->dim(1),$x->dim(1));
326 $ess->slice($_,$_).=$s->slice($_) foreach (0..$x->dim(1)-1); #generic diagonal
327 $x_copy = $v x $ess x $u->transpose;
328 EXAMPLE
330 The computing literature has loads of examples of how to use SVD.
332 Here's a trivial example (used in PDL::Transform::map) of how to make a
333 matrix less, er, singular, without changing the orientation of the
334 ellipsoid of transformation:
335 { my($r1,$s,$r2) = svd $x;
337 $s++; # fatten all singular values
338 $r2 *= $s; # implicit broadcasting for cheap mult.
339 $x .= $r2 x $r1; # a gets r2 x ess x r1
340 }
341 svd ignores the bad-value flag of the input ndarrays. It will set the
343 bad-value flag of all output ndarrays if the flag is set for any of the
344 input ndarrays.
345 lu_decomp
347 Signature: (a(m,m); [o]lu(m,m); [o]perm(m); [o]parity)
348 LU decompose a matrix, with row permutation
350 ($lu, $perm, $parity) = lu_decomp($x);
352 $lu = lu_decomp($x, $perm, $par); # $perm and $par are outputs!
354 lu_decomp($x->inplace,$perm,$par); # Everything in place.
356 "lu_decomp" returns an LU decomposition of a square matrix, using
358 Crout's method with partial pivoting. It's ported from Numerical
359 Recipes. The partial pivoting keeps it numerically stable but means a
360 little more overhead from broadcasting.
361 "lu_decomp" decomposes the input matrix into matrices L and U such that
363 LU = A, L is a subdiagonal matrix, and U is a superdiagonal matrix. By
364 convention, the diagonal of L is all 1's.
365 The single output matrix contains all the variable elements of both the
367 L and U matrices, stacked together. Because the method uses pivoting
368 (rearranging the lower part of the matrix for better numerical
369 stability), you have to permute input vectors before applying the L and
370 U matrices. The permutation is returned either in the second argument
371 or, in list context, as the second element of the list. You need the
372 permutation for the output to make any sense, so be sure to get it one
373 way or the other.
374 LU decomposition is the answer to a lot of matrix questions, including
376 inversion and determinant-finding, and "lu_decomp" is used by "inv".
377 If you pass in $perm and $parity, they either must be predeclared PDLs
379 of the correct size ($perm is an n-vector, $parity is a scalar) or
380 scalars.
381 If the matrix is singular, then the LU decomposition might not be
383 defined; in those cases, "lu_decomp" silently returns undef. Some
384 singular matrices LU-decompose just fine, and those are handled OK but
385 give a zero determinant (and hence can't be inverted).
386 "lu_decomp" uses pivoting, which rearranges the values in the matrix
388 for more numerical stability. This makes it really good for large and
389 even near-singular matrices. There is a non-pivoting version
390 "lu_decomp2" available which is from 5 to 60 percent faster for typical
391 problems at the expense of failing to compute a result in some cases.
392 Now that the "lu_decomp" is broadcasted, it is the recommended LU
394 decomposition routine. It no longer falls back to "lu_decomp2".
395 "lu_decomp" is ported from Numerical Recipes to PDL. It should probably
397 be implemented in C.
398 lu_decomp2
400 Signature: (a(m,m); [o]lu(m,m))
401 LU decompose a matrix, with no row permutation
403 ($lu, $perm, $parity) = lu_decomp2($x);
405 $lu = lu_decomp2($x,$perm,$parity); # or
407 $lu = lu_decomp2($x); # $perm and $parity are optional
408 lu_decomp($x->inplace,$perm,$parity); # or
410 lu_decomp($x->inplace); # $perm and $parity are optional
411 "lu_decomp2" works just like "lu_decomp", but it does no pivoting at
413 all. For compatibility with "lu_decomp", it will give you a
414 permutation list and a parity scalar if you ask for them -- but they
415 are always trivial.
416 Because "lu_decomp2" does not pivot, it is numerically unstable -- that
418 means it is less precise than "lu_decomp", particularly for large or
419 near-singular matrices. There are also specific types of non-singular
420 matrices that confuse it (e.g. ([0,-1,0],[1,0,0],[0,0,1]), which is a
421 90 degree rotation matrix but which confuses "lu_decomp2").
422 On the other hand, if you want to invert rapidly a few hundred thousand
424 small matrices and don't mind missing one or two, it could be the
425 ticket. It can be up to 60% faster at the expense of possible failure
426 of the decomposition for some of the input matrices.
427 The output is a single matrix that contains the LU decomposition of $a;
429 you can even do it in-place, thereby destroying $a, if you want. See
430 "lu_decomp" for more information about LU decomposition.
431 "lu_decomp2" is ported from Numerical Recipes into PDL.
433 lu_backsub
435 Signature: (lu(m,m); perm(m); b(m))
436 Solve A x = B for matrix A, by back substitution into A's LU
438 decomposition.
439 ($lu,$perm,$par) = lu_decomp($A);
441 $x = lu_backsub($lu,$perm,$par,$A); # or
443 $x = lu_backsub($lu,$perm,$B); # $par is not required for lu_backsub
444 lu_backsub($lu,$perm,$B->inplace); # modify $B in-place
446 $x = lu_backsub(lu_decomp($A),$B); # (ignores parity value from lu_decomp)
448 # starting from square matrix A and columns matrix B, with mathematically
450 # correct dimensions
451 $A = identity(4) + ones(4, 4);
452 $A->slice('2,0') .= 0; # break symmetry to see if need transpose
453 $B = sequence(2, 4);
454 # all these functions take B as rows, interpret as though notional columns
455 # mathematically confusing but can't change as back-compat and also
456 # familiar to Fortran users, so just transpose inputs and outputs
457 # using lu_backsub
459 ($lu,$perm,$par) = lu_decomp($A);
460 $x = lu_backsub($lu,$perm,$par, $B->transpose)->transpose;
461 # or with Slatec LINPACK
463 use PDL::Slatec;
464 gefa($lu=$A->copy, $ipiv=null, $info=null);
465 # 1 = do transpose because Fortran's idea of rows vs columns
466 gesl($lu, $ipiv, $x=$B->transpose->copy, 1);
467 $x = $x->inplace->transpose;
468 # or with LAPACK
470 use PDL::LinearAlgebra::Real;
471 getrf($lu=$A->copy, $ipiv=null, $info=null);
472 getrs($lu, 1, $x=$B->transpose->copy, $ipiv, $info=null); # again, need transpose
473 $x=$x->inplace->transpose;
474 # or with GSL
476 use PDL::GSL::LINALG;
477 LU_decomp(my $lu=$A->copy, my $p=null, my $signum=null);
478 # $B and $x, first dim is because GSL treats as vector, higher dims broadcast
479 # so we transpose in and back out
480 LU_solve($lu, $p, $B->transpose, my $x=null);
481 $x=$x->inplace->transpose;
482 # proof of the pudding is in the eating:
484 print $A x $x;
485 Given the LU decomposition of a square matrix (from "lu_decomp"),
487 "lu_backsub" does back substitution into the matrix to solve "a x = b"
488 for given vector "b". It is separated from the "lu_decomp" method so
489 that you can call the cheap "lu_backsub" multiple times and not have to
490 do the expensive LU decomposition more than once.
491 "lu_backsub" acts on single vectors and broadcasts in the usual way,
493 which means that it treats $y as the transpose of the input. If you
494 want to process a matrix, you must hand in the transpose of the matrix,
495 and then transpose the output when you get it back. that is because
496 pdls are indexed by (col,row), and matrices are (row,column) by
497 convention, so a 1-D pdl corresponds to a row vector, not a column
498 vector.
499 If $lu is dense and you have more than a few points to solve for, it is
501 probably cheaper to find "a^-1" with "inv", and just multiply "x = a^-1
502 b".) in fact, "inv" works by calling "lu_backsub" with the identity
503 matrix.
504 "lu_backsub" is ported from section 2.3 of Numerical Recipes. It is
506 written in PDL but should probably be implemented in C.
507 simq
509 Signature: ([phys]a(n,n); [phys]b(n); [o,phys]x(n); int [o,phys]ips(n); int flag)
510 Solution of simultaneous linear equations, "a x = b".
512 $a is an "n x n" matrix (i.e., a vector of length "n*n"), stored row-
514 wise: that is, "a(i,j) = a[ij]", where "ij = i*n + j".
515 While this is the transpose of the normal column-wise storage, this
517 corresponds to normal PDL usage. The contents of matrix a may be
518 altered (but may be required for subsequent calls with flag = -1).
519 $y, $x, $ips are vectors of length "n".
521 Set "flag=0" to solve. Set "flag=-1" to do a new back substitution for
523 different $y vector using the same a matrix previously reduced when
524 "flag=0" (the $ips vector generated in the previous solution is also
525 required).
526 See also "lu_backsub", which does the same thing with a slightly less
528 opaque interface.
529 simq ignores the bad-value flag of the input ndarrays. It will set the
531 bad-value flag of all output ndarrays if the flag is set for any of the
532 input ndarrays.
533 squaretotri
535 Signature: (a(n,n); b(m))
536 Convert a symmetric square matrix to triangular vector storage.
538 squaretotri does not process bad values. It will set the bad-value
540 flag of all output ndarrays if the flag is set for any of the input
541 ndarrays.
542
544 Copyright (C) 2002 Craig DeForest (deforest@boulder.swri.edu), R.J.R.
545 Williams (rjrw@ast.leeds.ac.uk), Karl Glazebrook
546 (kgb@aaoepp.aao.gov.au). There is no warranty. You are allowed to
547 redistribute and/or modify this work under the same conditions as PDL
548 itself. If this file is separated from the PDL distribution, then the
549 PDL copyright notice should be included in this file.
550perl v5.34.0 2022-02-28 MatrixOps(3)