1Math::MatrixReal(3) User Contributed Perl Documentation Math::MatrixReal(3)
2
6 Math::MatrixReal - Matrix of Reals
7 Implements the data type "matrix of reals" (and consequently also "vec‐
9 tor of reals").
10
12 Implements the data type "matrix of reals", which can be used almost
13 like any other basic Perl type thanks to OPERATOR OVERLOADING, i.e.,
14 $product = $matrix1 * $matrix2;
16 does what you would like it to do (a matrix multiplication).
18 Also features many important operations and methods: matrix norm,
20 matrix transposition, matrix inverse, determinant of a matrix, order
21 and numerical condition of a matrix, scalar product of vectors, vector
22 product of vectors, vector length, projection of row and column vec‐
23 tors, a comfortable way for reading in a matrix from a file, the key‐
24 board or your code, and many more.
25 Allows to solve linear equation systems using an efficient algorithm
27 known as "L-R-decomposition" and several approximative (iterative)
28 methods.
29 Features an implementation of Kleene's algorithm to compute the minimal
31 costs for all paths in a graph with weighted edges (the "weights" being
32 the costs associated with each edge).
33
35 Constructor Methods And Such
36 "use Math::MatrixReal;"
38 Makes the methods and overloaded operators of this module available
40 to your program.
41 · "$new_matrix = new Math::MatrixReal($rows,$columns);"
43 The matrix object constructor method. A new matrix of size $rows by
45 $columns will be created, with the value 0.0 for all elements.
46 Note that this method is implicitly called by many of the other
48 methods in this module.
49 · "$new_matrix = $some_matrix->""new($rows,$columns);"
51 Another way of calling the matrix object constructor method.
53 Matrix "$some_matrix" is not changed by this in any way.
55 · "$new_matrix = $matrix->new_from_cols( [ $column_vec‐
57 tor⎪$array_ref⎪$string, ... ] )"
58 Creates a new matrix given a reference to an array of any of the
60 following:
61 * column vectors ( n by 1 Math::MatrixReal matrices )
63 * references to arrays
64 * strings properly formatted to create a column with
65 Math::MatrixReal's "new_from_string" command
66 You may mix and match these as you wish. However, all must be of
68 the same dimension--no padding happens automatically. Example:
69 my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
71 print $matrix;
72 will print
74 [ 1.000000000000E+00 3.000000000000E+00 ]
76 [ 2.000000000000E+00 4.000000000000E+00 ]
77 · "new_from_rows( [ $row_vector⎪$array_ref⎪$string, ... ] )"
79 Creates a new matrix given a reference to an array of any of the
81 following:
82 * row vectors ( 1 by n Math::MatrixReal matrices )
84 * references to arrays
85 * strings properly formatted to create a row with
86 Math::MatrixReal's "new_from_string" command
87 You may mix and match these as you wish. However, all must be of
89 the same dimension--no padding happens automatically. Example:
90 my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
92 print $matrix;
93 will print
95 [ 1.000000000000E+00 2.000000000000E+00 ]
97 [ 3.000000000000E+00 4.000000000000E+00 ]
98 · "$new_matrix = Math::MatrixReal->new_diag( $array_ref );"
100 This method allows you to create a diagonal matrix by only specify‐
102 ing the diagonal elements. Example:
103 $matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
105 print $matrix;
106 will print
108 [ 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
110 [ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
111 [ 0.000000000000E+00 0.000000000000E+00 3.000000000000E+00 0.000000000000E+00 ]
112 [ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 4.000000000000E+00 ]
113 · "$new_matrix = Math::MatrixReal->""new_from_string($string);"
115 This method allows you to read in a matrix from a string (for
117 instance, from the keyboard, from a file or from your code).
118 The syntax is simple: each row must start with ""[ "" and end with
120 "" ]\n"" (""\n"" being the newline character and "" "" a space or
121 tab) and contain one or more numbers, all separated from each other
122 by spaces or tabs.
123 Additional spaces or tabs can be added at will, but no comments.
125 Examples:
127 $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
129 $matrix = Math::MatrixReal->new_from_string($string);
130 print "$matrix";
131 By the way, this prints
133 [ 1.000000000000E+00 2.000000000000E+00 3.000000000000E+00 ]
135 [ 2.000000000000E+00 2.000000000000E+00 -1.000000000000E+00 ]
136 [ 1.000000000000E+00 1.000000000000E+00 1.000000000000E+00 ]
137 But you can also do this in a much more comfortable way using the
139 shell-like "here-document" syntax:
140 $matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
142 [ 1 0 0 0 0 0 1 ]
143 [ 0 1 0 0 0 0 0 ]
144 [ 0 0 1 0 0 0 0 ]
145 [ 0 0 0 1 0 0 0 ]
146 [ 0 0 0 0 1 0 0 ]
147 [ 0 0 0 0 0 1 0 ]
148 [ 1 0 0 0 0 0 -1 ]
149 MATRIX
150 You can even use variables in the matrix:
152 $c1 = 2 / 3;
154 $c2 = -2 / 5;
155 $c3 = 26 / 9;
156 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
158 [ 3 2 0 ]
160 [ 0 3 2 ]
161 [ $c1 $c2 $c3 ]
162 MATRIX
164 (Remember that you may use spaces and tabs to format the matrix to
166 your taste)
167 Note that this method uses exactly the same representation for a
169 matrix as the "stringify" operator "": this means that you can con‐
170 vert any matrix into a string with "$string = "$matrix";" and read
171 it back in later (for instance from a file!).
172 Note however that you may suffer a precision loss in this process
174 because only 13 digits are supported in the mantissa when printed!!
175 If the string you supply (or someone else supplies) does not obey
177 the syntax mentioned above, an exception is raised, which can be
178 caught by "eval" as follows:
179 print "Please enter your matrix (in one line): ";
181 $string = <STDIN>;
182 $string =~ s/\\n/\n/g;
183 eval { $matrix = Math::MatrixReal->new_from_string($string); };
184 if ($@)
185 {
186 print "$@";
187 # ...
188 # (error handling)
189 }
190 else
191 {
192 # continue...
193 }
194 or as follows:
196 eval { $matrix = Math::MatrixReal->new_from_string(<<"MATRIX"); };
198 [ 3 2 0 ]
199 [ 0 3 2 ]
200 [ $c1 $c2 $c3 ]
201 MATRIX
202 if ($@)
203 # ...
204 Actually, the method shown above for reading a matrix from the key‐
206 board is a little awkward, since you have to enter a lot of "\n"'s
207 for the newlines.
208 A better way is shown in this piece of code:
210 while (1)
212 {
213 print "\nPlease enter your matrix ";
214 print "(multiple lines, <ctrl-D> = done):\n";
215 eval { $new_matrix =
216 Math::MatrixReal->new_from_string(join('',<STDIN>)); };
217 if ($@)
218 {
219 $@ =~ s/\s+at\b.*?$//;
220 print "${@}Please try again.\n";
221 }
222 else { last; }
223 }
224 Possible error messages of the "new_from_string()" method are:
226 Math::MatrixReal::new_from_string(): syntax error in input string
228 Math::MatrixReal::new_from_string(): empty input string
229 If the input string has rows with varying numbers of columns, the
231 following warning will be printed to STDERR:
232 Math::MatrixReal::new_from_string(): missing elements will be set to zero!
234 If everything is okay, the method returns an object reference to
236 the (newly allocated) matrix containing the elements you specified.
237 · "$new_matrix = $some_matrix->shadow();"
239 Returns an object reference to a NEW but EMPTY matrix (filled with
241 zero's) of the SAME SIZE as matrix "$some_matrix".
242 Matrix "$some_matrix" is not changed by this in any way.
244 · "$matrix1->copy($matrix2);"
246 Copies the contents of matrix "$matrix2" to an ALREADY EXISTING
248 matrix "$matrix1" (which must have the same size as matrix
249 "$matrix2"!).
250 Matrix "$matrix2" is not changed by this in any way.
252 · "$twin_matrix = $some_matrix->clone();"
254 Returns an object reference to a NEW matrix of the SAME SIZE as
256 matrix "$some_matrix". The contents of matrix "$some_matrix" have
257 ALREADY BEEN COPIED to the new matrix "$twin_matrix". This is the
258 method that the operator "=" is overloaded to when you type "$a =
259 $b", when $a and $b are matrices.
260 Matrix "$some_matrix" is not changed by this in any way.
262 Matrix Row, Column and Element operations
264 · "$row_vector = $matrix->row($row);"
266 This is a projection method which returns an object reference to a
268 NEW matrix (which in fact is a (row) vector since it has only one
269 row) to which row number "$row" of matrix "$matrix" has already
270 been copied.
271 Matrix "$matrix" is not changed by this in any way.
273 · "$column_vector = $matrix->column($column);"
275 This is a projection method which returns an object reference to a
277 NEW matrix (which in fact is a (column) vector since it has only
278 one column) to which column number "$column" of matrix "$matrix"
279 has already been copied.
280 Matrix "$matrix" is not changed by this in any way.
282 · "$matrix->assign($row,$column,$value);"
284 Explicitly assigns a value "$value" to a single element of the
286 matrix "$matrix", located in row "$row" and column "$column",
287 thereby replacing the value previously stored there.
288 · "$value = $matrix->""element($row,$column);"
290 Returns the value of a specific element of the matrix "$matrix",
292 located in row "$row" and column "$column".
293 · "$new_matrix = $matrix->each( \&function )";
295 Creates a new matrix by evaluating a code reference on each element
297 of the given matrix. The function is passed the element, the row
298 index and the column index, in that order. The value the function
299 returns ( or the value of the last executed statement ) is the
300 value given to the corresponding element in $new_matrix.
301 Example:
303 # add 1 to every element in the matrix
305 $matrix = $matrix->each ( sub { (shift) + 1 } );
306 Example:
308 my $cofactor = $matrix->each( sub { my(undef,$i,$j) = @_;
310 ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
311 : -1*$matrix->minor($i,$j)->det();
312 } );
313 This code needs some explanation. For each element of $matrix, it
315 throws away the actual value and stores the row and column indexes
316 in $i and $j. Then it sets element [$i,$j] in $cofactor to the
317 determinant of "$matrix->minor($i,$j)" if it is an "even" element,
318 or "-1*$matrix->minor($i,$j)" if it is an "odd" element.
319 · "$new_matrix = $matrix->each_diag( \&function )";
321 Creates a new matrix by evaluating a code reference on each diago‐
323 nal element of the given matrix. The function is passed the ele‐
324 ment, the row index and the column index, in that order. The value
325 the function returns ( or the value of the last executed statement
326 ) is the value given to the corresponding element in $new_matrix.
327 · "$matrix->swap_col( $col1, $col2 );"
329 This method takes two one-based column numbers and swaps the values
331 of each element in each column. "$matrix->swap_col(2,3)" would
332 replace column 2 in $matrix with column 3, and replace column 3
333 with column 2.
334 · "$matrix->swap_row( $row1, $row2 );"
336 This method takes two one-based row numbers and swaps the values of
338 each element in each row. "$matrix->swap_row(2,3)" would replace
339 row 2 in $matrix with row 3, and replace row 3 with row 2.
340 Matrix Operations
342 · "$det = $matrix->det();"
344 Returns the determinant of the matrix, without going through the
346 rigamarole of computing a LR decomposition. This method should be
347 much faster than LR decomposition if the matrix is diagonal or tri‐
348 angular. Otherwise, it is just a wrapper for "$matrix->decom‐
349 pose_LR->det_LR". If the determinant is zero, there is no inverse
350 and vice-versa. Only quadratic matrices have determinants.
351 · "$inverse = $matrix->inverse();"
353 Returns the inverse of a matrix, without going through the rigama‐
355 role of computing a LR decomposition. If no inverse exists, undef
356 is returned and an error is printed via "carp()". This is nothing
357 but a wrapper for "$matrix->decompose_LR->invert_LR".
358 · "($rows,$columns) = $matrix->dim();"
360 Returns a list of two items, representing the number of rows and
362 columns the given matrix "$matrix" contains.
363 · "$norm_one = $matrix->norm_one();"
365 Returns the "one"-norm of the given matrix "$matrix".
367 The "one"-norm is defined as follows:
369 For each column, the sum of the absolute values of the elements in
371 the different rows of that column is calculated. Finally, the maxi‐
372 mum of these sums is returned.
373 Note that the "one"-norm and the "maximum"-norm are mathematically
375 equivalent, although for the same matrix they usually yield a dif‐
376 ferent value.
377 Therefore, you should only compare values that have been calculated
379 using the same norm!
380 Throughout this package, the "one"-norm is (arbitrarily) used for
382 all comparisons, for the sake of uniformity and comparability,
383 except for the iterative methods "solve_GSM()", "solve_SSM()" and
384 "solve_RM()" which use either norm depending on the matrix itself.
385 · "$norm_max = $matrix->norm_max();"
387 Returns the "maximum"-norm of the given matrix $matrix.
389 The "maximum"-norm is defined as follows:
391 For each row, the sum of the absolute values of the elements in the
393 different columns of that row is calculated. Finally, the maximum
394 of these sums is returned.
395 Note that the "maximum"-norm and the "one"-norm are mathematically
397 equivalent, although for the same matrix they usually yield a dif‐
398 ferent value.
399 Therefore, you should only compare values that have been calculated
401 using the same norm!
402 Throughout this package, the "one"-norm is (arbitrarily) used for
404 all comparisons, for the sake of uniformity and comparability,
405 except for the iterative methods "solve_GSM()", "solve_SSM()" and
406 "solve_RM()" which use either norm depending on the matrix itself.
407 · "$norm_sum = $matrix->norm_sum();"
409 This is a very simple norm which is defined as the sum of the abso‐
411 lute values of every element.
412 · $p_norm = $matrix->norm_p($n);>
414 This function returns the "p-norm" of a vector. The argument $n
416 must be a number greater than or equal to 1 or the string "Inf".
417 The p-norm is defined as (sum(x_i^p))^(1/p). In words, it raised
418 each element to the p-th power, adds them up, and then takes the
419 p-th root of that number. If the string "Inf" is passed, the
420 "infinity-norm" is computed, which is really the limit of the
421 p-norm as p goes to infinity. It is defined as the maximum element
422 of the vector. Also, note that the familiar Euclidean distance
423 between two vectors is just a special case of a p-norm, when p is
424 equal to 2.
425 Example:
427 $a = Math::MatrixReal->new_from_cols([[1,2,3]]);
428 $p1 = $a->norm_p(1);
429 $p2 = $a->norm_p(2);
430 $p3 = $a->norm_p(3);
431 $pinf = $a->norm_p("Inf");
432 print "(1,2,3,Inf) norm:\n$p1\n$p2\n$p3\n$pinf\n";
434 $i1 = $a->new_from_rows([[1,0]]);
436 $i2 = $a->new_from_rows([[0,1]]);
437 # this should be sqrt(2) since it is the same as the
439 # hypotenuse of a 1 by 1 right triangle
440 $dist = ($i1-$i2)->norm_p(2);
442 print "Distance is $dist, which should be " . sqrt(2) . "\n";
443 Output:
445 (1,2,3,Inf) norm:
447 6
448 3.74165738677394139
449 3.30192724889462668
450 3
451 Distance is 1.41421356237309505, which should be 1.41421356237309505
453 · $frob_norm = "$matrix->norm_frobenius();"
455 This norm is similar to that of a p-norm where p is 2, except it
457 acts on a matrix, not a vector. Each element of the matrix is
458 squared, this is added up, and then a square root is taken.
459 · "$matrix->spectral_radius();"
461 Returns the maximum value of the absolute value of all eigenvalues.
463 Currently this computes all eigenvalues, then sifts through them to
464 find the largest in absolute value. Needless to say, this is very
465 inefficient, and in the future an algorithm that computes only the
466 largest eigenvalue may be implemented.
467 · "$matrix1->transpose($matrix2);"
469 Calculates the transposed matrix of matrix $matrix2 and stores the
471 result in matrix "$matrix1" (which must already exist and have the
472 same size as matrix "$matrix2"!).
473 This operation can also be carried out "in-place", i.e., input and
475 output matrix may be identical.
476 Transposition is a symmetry operation: imagine you rotate the
478 matrix along the axis of its main diagonal (going through elements
479 (1,1), (2,2), (3,3) and so on) by 180 degrees.
480 Another way of looking at it is to say that rows and columns are
482 swapped. In fact the contents of element "(i,j)" are swapped with
483 those of element "(j,i)".
484 Note that (especially for vectors) it makes a big difference if you
486 have a row vector, like this:
487 [ -1 0 1 ]
489 or a column vector, like this:
491 [ -1 ]
493 [ 0 ]
494 [ 1 ]
495 the one vector being the transposed of the other!
497 This is especially true for the matrix product of two vectors:
499 [ -1 ]
501 [ -1 0 1 ] * [ 0 ] = [ 2 ] , whereas
502 [ 1 ]
503 * [ -1 0 1 ]
505 [ -1 ] [ 1 0 -1 ]
506 [ 0 ] * [ -1 0 1 ] = [ -1 ] [ 1 0 -1 ] = [ 0 0 0 ]
507 [ 1 ] [ 0 ] [ 0 0 0 ] [ -1 0 1 ]
508 [ 1 ] [ -1 0 1 ]
509 So be careful about what you really mean!
511 Hint: throughout this module, whenever a vector is explicitly
513 required for input, a COLUMN vector is expected!
514 · "$trace = $matrix->trace();"
516 This returns the trace of the matrix, which is defined as the sum
518 of the diagonal elements. The matrix must be quadratic.
519 · "$minor = $matrix->minor($row,$col);"
521 Returns the minor matrix corresponding to $row and $col. $matrix
523 must be quadratic. If $matrix is n rows by n cols, the minor of
524 $row and $col will be an (n-1) by (n-1) matrix. The minor is
525 defined as crossing out the row and the col specified and returning
526 the remaining rows and columns as a matrix. This method is used by
527 "cofactor()".
528 · "$cofactor = $matrix->cofactor();"
530 The cofactor matrix is constructed as follows:
532 For each element, cross out the row and column that it sits in.
534 Now, take the determinant of the matrix that is left in the other
535 rows and columns. Multiply the determinant by (-1)^(i+j), where i
536 is the row index, and j is the column index. Replace the given
537 element with this value.
538 The cofactor matrix can be used to find the inverse of the matrix.
540 One formula for the inverse of a matrix is the cofactor matrix
541 transposed divided by the original determinant of the matrix.
542 The following two inverses should be exactly the same:
544 my $inverse1 = $matrix->inverse;
546 my $inverse2 = ~($matrix->cofactor)->each( sub { (shift)/$matrix->det() } );
547 Caveat: Although the cofactor matrix is simple algorithm to compute
549 the inverse of a matrix, and can be used with pencil and paper for
550 small matrices, it is comically slower than the native "inverse()"
551 function. Here is a small benchmark:
552 # $matrix1 is 15x15
554 $det = $matrix1->det;
555 timethese( 10,
556 {'inverse' => sub { $matrix1->inverse(); },
557 'cofactor' => sub { (~$matrix1->cofactor)->each ( sub { (shift)/$det; } ) }
558 } );
559 Benchmark: timing 10 iterations of LR, cofactor, inverse...
561 inverse: 1 wallclock secs ( 0.56 usr + 0.00 sys = 0.56 CPU) @ 17.86/s (n=10)
562 cofactor: 36 wallclock secs (36.62 usr + 0.01 sys = 36.63 CPU) @ 0.27/s (n=10)
563 · "$adjoint = $matrix->adjoint();"
565 The adjoint is just the transpose of the cofactor matrix. This
567 method is just an alias for " ~($matrix->cofactor)".
568 Arithmetic Operations
570 · "$matrix1->add($matrix2,$matrix3);"
572 Calculates the sum of matrix "$matrix2" and matrix "$matrix3" and
574 stores the result in matrix "$matrix1" (which must already exist
575 and have the same size as matrix "$matrix2" and matrix
576 "$matrix3"!).
577 This operation can also be carried out "in-place", i.e., the output
579 and one (or both) of the input matrices may be identical.
580 · "$matrix1->subtract($matrix2,$matrix3);"
582 Calculates the difference of matrix "$matrix2" minus matrix
584 "$matrix3" and stores the result in matrix "$matrix1" (which must
585 already exist and have the same size as matrix "$matrix2" and
586 matrix "$matrix3"!).
587 This operation can also be carried out "in-place", i.e., the output
589 and one (or both) of the input matrices may be identical.
590 Note that this operation is the same as
592 "$matrix1->add($matrix2,-$matrix3);", although the latter is a lit‐
593 tle less efficient.
594 · "$matrix1->multiply_scalar($matrix2,$scalar);"
596 Calculates the product of matrix "$matrix2" and the number
598 "$scalar" (i.e., multiplies each element of matrix "$matrix2" with
599 the factor "$scalar") and stores the result in matrix "$matrix1"
600 (which must already exist and have the same size as matrix
601 "$matrix2"!).
602 This operation can also be carried out "in-place", i.e., input and
604 output matrix may be identical.
605 · "$product_matrix = $matrix1->multiply($matrix2);"
607 Calculates the product of matrix "$matrix1" and matrix "$matrix2"
609 and returns an object reference to a new matrix "$product_matrix"
610 in which the result of this operation has been stored.
611 Note that the dimensions of the two matrices "$matrix1" and
613 "$matrix2" (i.e., their numbers of rows and columns) must harmonize
614 in the following way (example):
615 [ 2 2 ]
617 [ 2 2 ]
618 [ 2 2 ]
619 [ 1 1 1 ] [ * * ]
621 [ 1 1 1 ] [ * * ]
622 [ 1 1 1 ] [ * * ]
623 [ 1 1 1 ] [ * * ]
624 I.e., the number of columns of matrix "$matrix1" has to be the same
626 as the number of rows of matrix "$matrix2".
627 The number of rows and columns of the resulting matrix "$prod‐
629 uct_matrix" is determined by the number of rows of matrix
630 "$matrix1" and the number of columns of matrix "$matrix2", respec‐
631 tively.
632 · "$matrix1->negate($matrix2);"
634 Calculates the negative of matrix "$matrix2" (i.e., multiplies all
636 elements with "-1") and stores the result in matrix "$matrix1"
637 (which must already exist and have the same size as matrix
638 "$matrix2"!).
639 This operation can also be carried out "in-place", i.e., input and
641 output matrix may be identical.
642 · "$matrix_to_power = $matrix1->exponent($integer);"
644 Raises the matrix to the $integer power. Obviously, $integer must
646 be an integer. If it is zero, the identity matrix is returned. If a
647 negative integer is given, the inverse will be computed (if it
648 exists) and then raised the the absolute value of $integer. The
649 matrix must be quadratic.
650 Boolean Matrix Operations
652 · "$matrix->is_quadratic();"
654 Returns a boolean value indicating if the given matrix is quadratic
656 (also know as "square" or "n by n"). A matrix is quadratic if it
657 has the same number of rows as it does columns.
658 · "$matrix->is_square();"
660 This is an alias for "is_quadratic()".
662 · "$matrix->is_symmetric();"
664 Returns a boolean value indicating if the given matrix is symmet‐
666 ric. By definition, a matrix is symmetric if and only if
667 (M[i,j]=M[j,i]). This is equivalent to "($matrix == ~$matrix)" but
668 without memory allocation. Only quadratic matrices can be symmet‐
669 ric.
670 Notes: A symmetric matrix always has real eigenvalues/eigenvectors.
672 A matrix plus its transpose is always symmetric.
673 · "$matrix->is_skew_symmetric();"
675 Returns a boolean value indicating if the given matrix is skew sym‐
677 metric. By definition, a matrix is symmetric if and only if
678 (M[i,j]=-M[j,i]). This is equivalent to "($matrix == -(~$matrix))"
679 but without memory allocation. Only quadratic matrices can be skew
680 symmetric.
681 · "$matrix->is_diagonal();"
683 Returns a boolean value indicating if the given matrix is diagonal,
685 i.e. all of the nonzero elements are on the main diagonal. Only
686 quadratic matrices can be diagonal.
687 · "$matrix->is_tridiagonal();"
689 Returns a boolean value indicating if the given matrix is tridiago‐
691 nal, i.e. all of the nonzero elements are on the main diagonal or
692 the diagonals above and below the main diagonal. Only quadratic
693 matrices can be tridiagonal.
694 · "$matrix->is_upper_triangular();"
696 Returns a boolean value indicating if the given matrix is upper
698 triangular, i.e. all of the nonzero elements not on the main diago‐
699 nal are above it. Only quadratic matrices can be upper triangular.
700 Note: diagonal matrices are both upper and lower triangular.
701 · "$matrix->is_lower_triangular();"
703 Returns a boolean value indicating if the given matrix is lower
705 triangular, i.e. all of the nonzero elements not on the main diago‐
706 nal are below it. Only quadratic matrices can be lower triangular.
707 Note: diagonal matrices are both upper and lower triangular.
708 · "$matrix->is_orthogonal();"
710 Returns a boolean value indicating if the given matrix is orthogo‐
712 nal. An orthogonal matrix is has the property that the transpose
713 equals the inverse of the matrix. Instead of computing each and
714 comparing them, this method multiplies the matrix by it's trans‐
715 pose, and returns true if this turns out to be the identity matrix,
716 false otherwise. Only quadratic matrices can orthogonal.
717 · "$matrix->is_binary();"
719 Returns a boolean value indicating if the given matrix is binary.
721 A matrix is binary if it contains only zeroes or ones.
722 · "$matrix->is_gramian();"
724 Returns a boolean value indicating if the give matrix is Gramian.
726 A matrix $A is Gramian if and only if there exists a square matrix
727 $B such that "$A = ~$B*$B". This is equivalent to checking if $A is
728 symmetric and has all nonnegative eigenvalues, which is what
729 Math::MatrixReal uses to check for this property.
730 · "$matrix->is_LR();"
732 Returns a boolean value indicating if the matrix is an LR decompo‐
734 sition matrix.
735 · "$matrix->is_positive();"
737 Returns a boolean value indicating if the matrix contains only pos‐
739 itive entries. Note that a zero entry is not positive and will
740 cause "is_positive()" to return false.
741 · "$matrix->is_negative();"
743 Returns a boolean value indicating if the matrix contains only neg‐
745 ative entries. Note that a zero entry is not negative and will
746 cause "is_negative()" to return false.
747 · "$matrix->is_periodic($k);"
749 Returns a boolean value indicating if the matrix is periodic with
751 period $k. This is true if "$matrix ** ($k+1) == $matrix". When
752 "$k == 1", this reduces down to the "is_idempotent()" function.
753 · "$matrix->is_idempotent();"
755 Returns a boolean value indicating if the matrix is idempotent,
757 which is defined as the square of the matrix being equal to the
758 original matrix, i.e "$matrix ** 2 == $matrix".
759 · "$matrix->is_row_vector();"
761 Returns a boolean value indicating if the matrix is a row vector.
763 A row vector is a matrix which is 1xn. Note that the 1x1 matrix is
764 both a row and column vector.
765 · "$matrix->is_col_vector();"
767 Returns a boolean value indicating if the matrix is a col vector.
769 A col vector is a matrix which is nx1. Note that the 1x1 matrix is
770 both a row and column vector.
771 Eigensystems
773 · "($l, $V) = $matrix->sym_diagonalize();"
775 This method performs the diagonalization of the quadratic symmet‐
777 ric matrix M stored in $matrix. On output, l is a column vector
778 containing all the eigenvalues of M and V is an orthogonal matrix
779 which columns are the corresponding normalized eigenvectors. The
780 primary property of an eigenvalue l and an eigenvector x is of
781 course that: M * x = l * x.
782 The method uses a Householder reduction to tridiagonal form fol‐
784 lowed by a QL algoritm with implicit shifts on this tridiagonal.
785 (The tridiagonal matrix is kept internally in a compact form in
786 this routine to save memory.) In fact, this routine wraps the
787 householder() and tri_diagonalize() methods described below when
788 their intermediate results are not desired. The overall algo‐
789 rithmic complexity of this technique is O(N^3). According to sev‐
790 eral books, the coefficient hidden by the 'O' is one of the best
791 possible for general (symmetric) matrixes.
792 · "($T, $Q) = $matrix->householder();"
794 This method performs the Householder algorithm which reduces the
796 n by n real symmetric matrix M contained in $matrix to tridiago‐
797 nal form. On output, T is a symmetric tridiagonal matrix (only
798 diagonal and off-diagonal elements are non-zero) and Q is an
799 orthogonal matrix performing the tranformation between M and T
800 ("$M == $Q * $T * ~$Q").
801 · "($l, $V) = $T->tri_diagonalize([$Q]);"
803 This method diagonalizes the symmetric tridiagonal matrix T. On
805 output, $l and $V are similar to the output values described for
806 sym_diagonalize().
807 The optional argument $Q corresponds to an orthogonal transforma‐
809 tion matrix Q that should be used additionally during V (eigen‐
810 vectors) computation. It should be supplied if the desired eigen‐
811 vectors correspond to a more general symmetric matrix M previ‐
812 ously reduced by the householder() method, not a mere tridiago‐
813 nal. If T is really a tridiagonal matrix, Q can be omitted (it
814 will be internally created in fact as an identity matrix). The
815 method uses a QL algorithm (with implicit shifts).
816 · "$l = $matrix->sym_eigenvalues();"
818 This method computes the eigenvalues of the quadratic symmetric
820 matrix M stored in $matrix. On output, l is a column vector con‐
821 taining all the eigenvalues of M. Eigenvectors are not computed
822 (on the contrary of "sym_diagonalize()") and this method is more
823 efficient (even though it uses a similar algorithm with two
824 phases). However, understand that the algorithmic complexity of
825 this technique is still also O(N^3). But the coefficient hidden
826 by the 'O' is better by a factor of..., well, see your benchmark,
827 it's wiser.
828 This routine wraps the householder_tridiagonal() and tri_eigen‐
830 values() methods described below when the intermediate tridiago‐
831 nal matrix is not needed.
832 · "$T = $matrix->householder_tridiagonal();"
834 This method performs the Householder algorithm which reduces the
836 n by n real symmetric matrix M contained in $matrix to tridiago‐
837 nal form. On output, T is the obtained symmetric tridiagonal
838 matrix (only diagonal and off-diagonal elements are non-zero).
839 The operation is similar to the householder() method, but poten‐
840 tially a little more efficient as the transformation matrix is
841 not computed.
842 · "$l = $T->tri_eigenvalues();"
844 This method computesthe eigenvalues of the symmetric tridiagonal
846 matrix T. On output, $l is a vector containing the eigenvalues
847 (similar to "sym_eigenvalues()"). This method is much more effi‐
848 cient than tri_diagonalize() when eigenvectors are not needed.
849 Miscellaneous
851 · "$matrix->zero();"
853 Assigns a zero to every element of the matrix "$matrix", i.e.,
855 erases all values previously stored there, thereby effectively
856 transforming the matrix into a "zero"-matrix or "null"-matrix, the
857 neutral element of the addition operation in a Ring.
858 (For instance the (quadratic) matrices with "n" rows and columns
860 and matrix addition and multiplication form a Ring. Most prominent
861 characteristic of a Ring is that multiplication is not commutative,
862 i.e., in general, ""matrix1 * matrix2"" is not the same as
863 ""matrix2 * matrix1""!)
864 · "$matrix->one();"
866 Assigns one's to the elements on the main diagonal (elements (1,1),
868 (2,2), (3,3) and so on) of matrix "$matrix" and zero's to all oth‐
869 ers, thereby erasing all values previously stored there and trans‐
870 forming the matrix into a "one"-matrix, the neutral element of the
871 multiplication operation in a Ring.
872 (If the matrix is quadratic (which this method doesn't require,
874 though), then multiplying this matrix with itself yields this same
875 matrix again, and multiplying it with some other matrix leaves that
876 other matrix unchanged!)
877 · "$latex_string = $matrix->as_latex( align=> "c", format => "%s",
879 name => "" );"
880 This function returns the matrix as a LaTeX string. It takes a hash
882 as an argument which is used to control the style of the output.
883 The hash element "align" may be "c","l" or "r", corresponding to
884 center, left and right, respectively. The "format" element is a
885 format string that is given to "sprintf" to control the style of
886 number format, such a floating point or scientific notation. The
887 "name" element can be used so that a LaTeX string of "$name = " is
888 prepended to the string.
889 Example:
891 my $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
893 print $a->as_latex( ( format => "%.2f", align => "l",name => "A" ) );
894 Output:
896 $A = $ $
897 \left( \begin{array}{ll}
898 1.23&1.00 \\
899 5.68&2.00 \\
900 9.10&3.00
901 \end{array} \right)
902 $
903 · "$yacas_string = $matrix->as_yacas( format => "%s", name => "",
905 semi => 0 );"
906 This function returns the matrix as a string that can be read by
908 Yacas. It takes a hash as an an argument which controls the style
909 of the output. The "format" element is a format string that is
910 given to "sprintf" to control the style of number format, such a
911 floating point or scientific notation. The "name" element can be
912 used so that "$name = " is prepended to the string. The <semi> ele‐
913 ment can be set to 1 to that a semicolon is appended (so Matlab
914 does not print out the matrix.)
915 Example:
917 $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
919 print $a->as_yacas( ( format => "%.2f", align => "l",name => "A" ) );
920 Output:
922 A := {{1.23,1.00},{5.68,2.00},{9.10,3.00}}
924 · "$matlab_string = $matrix->as_matlab( format => "%s", name => "",
926 semi => 0 );"
927 This function returns the matrix as a string that can be read by
929 Matlab. It takes a hash as an an argument which controls the style
930 of the output. The "format" element is a format string that is
931 given to "sprintf" to control the style of number format, such a
932 floating point or scientific notation. The "name" element can be
933 used so that "$name = " is prepended to the string. The <semi> ele‐
934 ment can be set to 1 to that a semicolon is appended (so Matlab
935 does not print out the matrix.)
936 Example:
938 my $a = Math::MatrixReal->new_from_rows([[ 1.234, 5.678, 9.1011],[1,2,3]] );
940 print $a->as_matlab( ( format => "%.3f", name => "A",semi => 1 ) );
941 Output:
943 A = [ 1.234 5.678 9.101;
944 1.000 2.000 3.000];
945 · "$scilab_string = $matrix->as_scilab( format => "%s", name => "",
947 semi => 0 );"
948 This function is just an alias for "as_matlab()", since both Scilab
950 and Matlab have the same matrix format.
951 · "$minimum = Math::MatrixReal::min($number1,$number2);"
953 Returns the minimum of the two numbers ""number1"" and ""number2"".
955 · "$maximum = Math::MatrixReal::max($number1,$number2);"
957 Returns the maximum of the two numbers ""number1"" and ""number2"".
959 · "$minimal_cost_matrix = $cost_matrix->kleene();"
961 Copies the matrix "$cost_matrix" (which has to be quadratic!) to a
963 new matrix of the same size (i.e., "clones" the input matrix) and
964 applies Kleene's algorithm to it.
965 See Math::Kleene(3) for more details about this algorithm!
967 The method returns an object reference to the new matrix.
969 Matrix "$cost_matrix" is not changed by this method in any way.
971 · "($norm_matrix,$norm_vector) = $matrix->normalize($vector);"
973 This method is used to improve the numerical stability when solving
975 linear equation systems.
976 Suppose you have a matrix "A" and a vector "b" and you want to find
978 out a vector "x" so that "A * x = b", i.e., the vector "x" which
979 solves the equation system represented by the matrix "A" and the
980 vector "b".
981 Applying this method to the pair (A,b) yields a pair (A',b') where
983 each row has been divided by (the absolute value of) the greatest
984 coefficient appearing in that row. So this coefficient becomes
985 equal to "1" (or "-1") in the new pair (A',b') (all others become
986 smaller than one and greater than minus one).
987 Note that this operation does not change the equation system itself
989 because the same division is carried out on either side of the
990 equation sign!
991 The method requires a quadratic (!) matrix "$matrix" and a vector
993 "$vector" for input (the vector must be a column vector with the
994 same number of rows as the input matrix) and returns a list of two
995 items which are object references to a new matrix and a new vector,
996 in this order.
997 The output matrix and vector are clones of the input matrix and
999 vector to which the operation explained above has been applied.
1000 The input matrix and vector are not changed by this in any way.
1002 Example of how this method can affect the result of the methods to
1004 solve equation systems (explained immediately below following this
1005 method):
1006 Consider the following little program:
1008 #!perl -w
1010 use Math::MatrixReal qw(new_from_string);
1012 $A = Math::MatrixReal->new_from_string(<<"MATRIX");
1014 [ 1 2 3 ]
1015 [ 5 7 11 ]
1016 [ 23 19 13 ]
1017 MATRIX
1018 $b = Math::MatrixReal->new_from_string(<<"MATRIX");
1020 [ 0 ]
1021 [ 1 ]
1022 [ 29 ]
1023 MATRIX
1024 $LR = $A->decompose_LR();
1026 if (($dim,$x,$B) = $LR->solve_LR($b))
1027 {
1028 $test = $A * $x;
1029 print "x = \n$x";
1030 print "A * x = \n$test";
1031 }
1032 ($A_,$b_) = $A->normalize($b);
1034 $LR = $A_->decompose_LR();
1036 if (($dim,$x,$B) = $LR->solve_LR($b_))
1037 {
1038 $test = $A * $x;
1039 print "x = \n$x";
1040 print "A * x = \n$test";
1041 }
1042 This will print:
1044 x =
1046 [ 1.000000000000E+00 ]
1047 [ 1.000000000000E+00 ]
1048 [ -1.000000000000E+00 ]
1049 A * x =
1050 [ 4.440892098501E-16 ]
1051 [ 1.000000000000E+00 ]
1052 [ 2.900000000000E+01 ]
1053 x =
1054 [ 1.000000000000E+00 ]
1055 [ 1.000000000000E+00 ]
1056 [ -1.000000000000E+00 ]
1057 A * x =
1058 [ 0.000000000000E+00 ]
1059 [ 1.000000000000E+00 ]
1060 [ 2.900000000000E+01 ]
1061 You can see that in the second example (where "normalize()" has
1063 been used), the result is "better", i.e., more accurate!
1064 · "$LR_matrix = $matrix->decompose_LR();"
1066 This method is needed to solve linear equation systems.
1068 Suppose you have a matrix "A" and a vector "b" and you want to find
1070 out a vector "x" so that "A * x = b", i.e., the vector "x" which
1071 solves the equation system represented by the matrix "A" and the
1072 vector "b".
1073 You might also have a matrix "A" and a whole bunch of different
1075 vectors "b1".."bk" for which you need to find vectors "x1".."xk" so
1076 that "A * xi = bi", for "i=1..k".
1077 Using Gaussian transformations (multiplying a row or column with a
1079 factor, swapping two rows or two columns and adding a multiple of
1080 one row or column to another), it is possible to decompose any
1081 matrix "A" into two triangular matrices, called "L" and "R" (for
1082 "Left" and "Right").
1083 "L" has one's on the main diagonal (the elements (1,1), (2,2),
1085 (3,3) and so so), non-zero values to the left and below of the main
1086 diagonal and all zero's in the upper right half of the matrix.
1087 "R" has non-zero values on the main diagonal as well as to the
1089 right and above of the main diagonal and all zero's in the lower
1090 left half of the matrix, as follows:
1091 [ 1 0 0 0 0 ] [ x x x x x ]
1093 [ x 1 0 0 0 ] [ 0 x x x x ]
1094 L = [ x x 1 0 0 ] R = [ 0 0 x x x ]
1095 [ x x x 1 0 ] [ 0 0 0 x x ]
1096 [ x x x x 1 ] [ 0 0 0 0 x ]
1097 Note that ""L * R"" is equivalent to matrix "A" in the sense that
1099 "L * R * x = b <==> A * x = b" for all vectors "x", leaving out
1100 of account permutations of the rows and columns (these are taken
1101 care of "magically" by this module!) and numerical errors.
1102 Trick:
1104 Because we know that "L" has one's on its main diagonal, we can
1106 store both matrices together in the same array without information
1107 loss! I.e.,
1108 [ R R R R R ]
1110 [ L R R R R ]
1111 LR = [ L L R R R ]
1112 [ L L L R R ]
1113 [ L L L L R ]
1114 Beware, though, that "LR" and ""L * R"" are not the same!!!
1116 Note also that for the same reason, you cannot apply the method
1118 "normalize()" to an "LR" decomposition matrix. Trying to do so will
1119 yield meaningless rubbish!
1120 (You need to apply "normalize()" to each pair (Ai,bi) BEFORE decom‐
1122 posing the matrix "Ai'"!)
1123 Now what does all this help us in solving linear equation systems?
1125 It helps us because a triangular matrix is the next best thing that
1127 can happen to us besides a diagonal matrix (a matrix that has non-
1128 zero values only on its main diagonal - in which case the solution
1129 is trivial, simply divide ""b[i]"" by ""A[i,i]"" to get ""x[i]""!).
1130 To find the solution to our problem ""A * x = b"", we divide this
1132 problem in parts: instead of solving "A * x = b" directly, we first
1133 decompose "A" into "L" and "R" and then solve ""L * y = b"" and
1134 finally ""R * x = y"" (motto: divide and rule!).
1135 From the illustration above it is clear that solving ""L * y = b""
1137 and ""R * x = y"" is straightforward: we immediately know that
1138 "y[1] = b[1]". We then deduce swiftly that
1139 y[2] = b[2] - L[2,1] * y[1]
1141 (and we know ""y[1]"" by now!), that
1143 y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
1145 and so on.
1147 Having effortlessly calculated the vector "y", we now proceed to
1149 calculate the vector "x" in a similar fashion: we see immediately
1150 that "x[n] = y[n] / R[n,n]". It follows that
1151 x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
1153 and
1155 x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
1157 / R[n-2,n-2]
1158 and so on.
1160 You can see that - especially when you have many vectors "b1".."bk"
1162 for which you are searching solutions to "A * xi = bi" - this
1163 scheme is much more efficient than a straightforward, "brute force"
1164 approach.
1165 This method requires a quadratic matrix as its input matrix.
1167 If you don't have that many equations, fill up with zero's (i.e.,
1169 do nothing to fill the superfluous rows if it's a "fresh" matrix,
1170 i.e., a matrix that has been created with "new()" or "shadow()").
1171 The method returns an object reference to a new matrix containing
1173 the matrices "L" and "R".
1174 The input matrix is not changed by this method in any way.
1176 Note that you can "copy()" or "clone()" the result of this method
1178 without losing its "magical" properties (for instance concerning
1179 the hidden permutations of its rows and columns).
1180 However, as soon as you are applying any method that alters the
1182 contents of the matrix, its "magical" properties are stripped off,
1183 and the matrix immediately reverts to an "ordinary" matrix (with
1184 the values it just happens to contain at that moment, be they mean‐
1185 ingful as an ordinary matrix or not!).
1186 · "($dimension,$x_vector,$base_matrix) =
1188 $LR_matrix""->""solve_LR($b_vector);"
1189 Use this method to actually solve an equation system.
1191 Matrix "$LR_matrix" must be a (quadratic) matrix returned by the
1193 method "decompose_LR()", the LR decomposition matrix of the matrix
1194 "A" of your equation system "A * x = b".
1195 The input vector "$b_vector" is the vector "b" in your equation
1197 system "A * x = b", which must be a column vector and have the same
1198 number of rows as the input matrix "$LR_matrix".
1199 The method returns a list of three items if a solution exists or an
1201 empty list otherwise (!).
1202 Therefore, you should always use this method like this:
1204 if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
1206 {
1207 # do something with the solution...
1208 }
1209 else
1210 {
1211 # do something with the fact that there is no solution...
1212 }
1213 The three items returned are: the dimension "$dimension" of the
1215 solution space (which is zero if only one solution exists, one if
1216 the solution is a straight line, two if the solution is a plane,
1217 and so on), the solution vector "$x_vector" (which is the vector
1218 "x" of your equation system "A * x = b") and a matrix
1219 "$base_matrix" representing a base of the solution space (a set of
1220 vectors which put up the solution space like the spokes of an
1221 umbrella).
1222 Only the first "$dimension" columns of this base matrix actually
1224 contain entries, the remaining columns are all zero.
1225 Now what is all this stuff with that "base" good for?
1227 The output vector "x" is ALWAYS a solution of your equation system
1229 "A * x = b".
1230 But also any vector "$vector"
1232 $vector = $x_vector->clone();
1234 $machine_infinity = 1E+99; # or something like that
1236 for ( $i = 1; $i <= $dimension; $i++ )
1238 {
1239 $vector += rand($machine_infinity) * $base_matrix->column($i);
1240 }
1241 is a solution to your problem "A * x = b", i.e., if "$A_matrix"
1243 contains your matrix "A", then
1244 print abs( $A_matrix * $vector - $b_vector ), "\n";
1246 should print a number around 1E-16 or so!
1248 By the way, note that you can actually calculate those vectors
1250 "$vector" a little more efficient as follows:
1251 $rand_vector = $x_vector->shadow();
1253 $machine_infinity = 1E+99; # or something like that
1255 for ( $i = 1; $i <= $dimension; $i++ )
1257 {
1258 $rand_vector->assign($i,1, rand($machine_infinity) );
1259 }
1260 $vector = $x_vector + ( $base_matrix * $rand_vector );
1262 Note that the input matrix and vector are not changed by this
1264 method in any way.
1265 · "$inverse_matrix = $LR_matrix->invert_LR();"
1267 Use this method to calculate the inverse of a given matrix
1269 "$LR_matrix", which must be a (quadratic) matrix returned by the
1270 method "decompose_LR()".
1271 The method returns an object reference to a new matrix of the same
1273 size as the input matrix containing the inverse of the matrix that
1274 you initially fed into "decompose_LR()" IF THE INVERSE EXISTS, or
1275 an empty list otherwise.
1276 Therefore, you should always use this method in the following way:
1278 if ( $inverse_matrix = $LR->invert_LR() )
1280 {
1281 # do something with the inverse matrix...
1282 }
1283 else
1284 {
1285 # do something with the fact that there is no inverse matrix...
1286 }
1287 Note that by definition (disregarding numerical errors), the prod‐
1289 uct of the initial matrix and its inverse (or vice-versa) is always
1290 a matrix containing one's on the main diagonal (elements (1,1),
1291 (2,2), (3,3) and so on) and zero's elsewhere.
1292 The input matrix is not changed by this method in any way.
1294 · "$condition = $matrix->condition($inverse_matrix);"
1296 In fact this method is just a shortcut for
1298 abs($matrix) * abs($inverse_matrix)
1300 Both input matrices must be quadratic and have the same size, and
1302 the result is meaningful only if one of them is the inverse of the
1303 other (for instance, as returned by the method "invert_LR()").
1304 The number returned is a measure of the "condition" of the given
1306 matrix "$matrix", i.e., a measure of the numerical stability of the
1307 matrix.
1308 This number is always positive, and the smaller its value, the bet‐
1310 ter the condition of the matrix (the better the stability of all
1311 subsequent computations carried out using this matrix).
1312 Numerical stability means for example that if
1314 abs( $vec_correct - $vec_with_error ) < $epsilon
1316 holds, there must be a "$delta" which doesn't depend on the vector
1318 "$vec_correct" (nor "$vec_with_error", by the way) so that
1319 abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
1321 also holds.
1323 · "$determinant = $LR_matrix->det_LR();"
1325 Calculates the determinant of a matrix, whose LR decomposition
1327 matrix "$LR_matrix" must be given (which must be a (quadratic)
1328 matrix returned by the method "decompose_LR()").
1329 In fact the determinant is a by-product of the LR decomposition: It
1331 is (in principle, that is, except for the sign) simply the product
1332 of the elements on the main diagonal (elements (1,1), (2,2), (3,3)
1333 and so on) of the LR decomposition matrix.
1334 (The sign is taken care of "magically" by this module)
1336 · "$order = $LR_matrix->order_LR();"
1338 Calculates the order (called "Rang" in German) of a matrix, whose
1340 LR decomposition matrix "$LR_matrix" must be given (which must be a
1341 (quadratic) matrix returned by the method "decompose_LR()").
1342 This number is a measure of the number of linear independent row
1344 and column vectors (= number of linear independent equations in the
1345 case of a matrix representing an equation system) of the matrix
1346 that was initially fed into "decompose_LR()".
1347 If "n" is the number of rows and columns of the (quadratic!)
1349 matrix, then "n - order" is the dimension of the solution space of
1350 the associated equation system.
1351 · "$rank = $LR_matrix->rank_LR();"
1353 This is an alias for the "order_LR()" function. The "order" is usu‐
1355 ally called the "rank" in the United States.
1356 · "$scalar_product = $vector1->scalar_product($vector2);"
1358 Returns the scalar product of vector "$vector1" and vector "$vec‐
1360 tor2".
1361 Both vectors must be column vectors (i.e., a matrix having several
1363 rows but only one column).
1364 This is a (more efficient!) shortcut for
1366 $temp = ~$vector1 * $vector2;
1368 $scalar_product = $temp->element(1,1);
1369 or the sum "i=1..n" of the products "vector1[i] * vector2[i]".
1371 Provided none of the two input vectors is the null vector, then the
1373 two vectors are orthogonal, i.e., have an angle of 90 degrees
1374 between them, exactly when their scalar product is zero, and
1375 vice-versa.
1376 · "$vector_product = $vector1->vector_product($vector2);"
1378 Returns the vector product of vector "$vector1" and vector "$vec‐
1380 tor2".
1381 Both vectors must be column vectors (i.e., a matrix having several
1383 rows but only one column).
1384 Currently, the vector product is only defined for 3 dimensions
1386 (i.e., vectors with 3 rows); all other vectors trigger an error
1387 message.
1388 In 3 dimensions, the vector product of two vectors "x" and "y" is
1390 defined as
1391 ⎪ x[1] y[1] e[1] ⎪
1393 determinant ⎪ x[2] y[2] e[2] ⎪
1394 ⎪ x[3] y[3] e[3] ⎪
1395 where the ""x[i]"" and ""y[i]"" are the components of the two vec‐
1397 tors "x" and "y", respectively, and the ""e[i]"" are unity vectors
1398 (i.e., vectors with a length equal to one) with a one in row "i"
1399 and zero's elsewhere (this means that you have numbers and vectors
1400 as elements in this matrix!).
1401 This determinant evaluates to the rather simple formula
1403 z[1] = x[2] * y[3] - x[3] * y[2]
1405 z[2] = x[3] * y[1] - x[1] * y[3]
1406 z[3] = x[1] * y[2] - x[2] * y[1]
1407 A characteristic property of the vector product is that the result‐
1409 ing vector is orthogonal to both of the input vectors (if neither
1410 of both is the null vector, otherwise this is trivial), i.e., the
1411 scalar product of each of the input vectors with the resulting vec‐
1412 tor is always zero.
1413 · "$length = $vector->length();"
1415 This is actually a shortcut for
1417 $length = sqrt( $vector->scalar_product($vector) );
1419 and returns the length of a given (column!) vector "$vector".
1421 Note that the "length" calculated by this method is in fact the
1423 "two"-norm of a vector "$vector"!
1424 The general definition for norms of vectors is the following:
1426 sub vector_norm
1428 {
1429 croak "Usage: \$norm = \$vector->vector_norm(\$n);"
1430 if (@_ != 2);
1431 my($vector,$n) = @_;
1433 my($rows,$cols) = ($vector->[1],$vector->[2]);
1434 my($k,$comp,$sum);
1435 croak "Math::MatrixReal::vector_norm(): vector is not a column vector"
1437 unless ($cols == 1);
1438 croak "Math::MatrixReal::vector_norm(): norm index must be > 0"
1440 unless ($n > 0);
1441 croak "Math::MatrixReal::vector_norm(): norm index must be integer"
1443 unless ($n == int($n));
1444 $sum = 0;
1446 for ( $k = 0; $k < $rows; $k++ )
1447 {
1448 $comp = abs( $vector->[0][$k][0] );
1449 $sum += $comp ** $n;
1450 }
1451 return( $sum ** (1 / $n) );
1452 }
1453 Note that the case "n = 1" is the "one"-norm for matrices applied
1455 to a vector, the case "n = 2" is the euclidian norm or length of a
1456 vector, and if "n" goes to infinity, you have the "infinity"- or
1457 "maximum"-norm for matrices applied to a vector!
1458 · "$xn_vector = $matrix->""solve_GSM($x0_vector,$b_vector,$epsilon);"
1460 · "$xn_vector = $matrix->""solve_SSM($x0_vector,$b_vector,$epsilon);"
1462 · "$xn_vector = $matrix->""solve_RM($x0_vector,$b_vec‐
1464 tor,$weight,$epsilon);"
1465 In some cases it might not be practical or desirable to solve an
1467 equation system ""A * x = b"" using an analytical algorithm like
1468 the "decompose_LR()" and "solve_LR()" method pair.
1469 In fact in some cases, due to the numerical properties (the "condi‐
1471 tion") of the matrix "A", the numerical error of the obtained
1472 result can be greater than by using an approximative (iterative)
1473 algorithm like one of the three implemented here.
1474 All three methods, GSM ("Global Step Method" or "Gesamtschrittver‐
1476 fahren"), SSM ("Single Step Method" or "Einzelschrittverfahren")
1477 and RM ("Relaxation Method" or "Relaxationsverfahren"), are fix-
1478 point iterations, that is, can be described by an iteration func‐
1479 tion ""x(t+1) = Phi( x(t) )"" which has the property:
1480 Phi(x) = x <==> A * x = b
1482 We can define "Phi(x)" as follows:
1484 Phi(x) := ( En - A ) * x + b
1486 where "En" is a matrix of the same size as "A" ("n" rows and col‐
1488 umns) with one's on its main diagonal and zero's elsewhere.
1489 This function has the required property.
1491 Proof:
1493 A * x = b
1495 <==> -( A * x ) = -b
1497 <==> -( A * x ) + x = -b + x
1499 <==> -( A * x ) + x + b = x
1501 <==> x - ( A * x ) + b = x
1503 <==> ( En - A ) * x + b = x
1505 This last step is true because
1507 x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
1509 is the same as
1511 ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
1513 qed
1515 Note that actually solving the equation system ""A * x = b"" means
1517 to calculate
1518 a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] = b[i]
1520 <==> a[i,i] x[i] =
1522 b[i]
1523 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
1524 + a[i,i] x[i]
1525 <==> x[i] =
1527 ( b[i]
1528 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
1529 + a[i,i] x[i]
1530 ) / a[i,i]
1531 <==> x[i] =
1533 ( b[i] -
1534 ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
1535 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
1536 ) / a[i,i]
1537 There is one major restriction, though: a fix-point iteration is
1539 guaranteed to converge only if the first derivative of the itera‐
1540 tion function has an absolute value less than one in an area around
1541 the point "x(*)" for which ""Phi( x(*) ) = x(*)"" is to be true,
1542 and if the start vector "x(0)" lies within that area!
1543 This is best verified grafically, which unfortunately is impossible
1545 to do in this textual documentation!
1546 See literature on Numerical Analysis for details!
1548 In our case, this restriction translates to the following three
1550 conditions:
1551 There must exist a norm so that the norm of the matrix of the iter‐
1553 ation function, "( En - A )", has a value less than one, the matrix
1554 "A" may not have any zero value on its main diagonal and the ini‐
1555 tial vector "x(0)" must be "good enough", i.e., "close enough" to
1556 the solution "x(*)".
1557 (Remember school math: the first derivative of a straight line
1559 given by ""y = a * x + b"" is "a"!)
1560 The three methods expect a (quadratic!) matrix "$matrix" as their
1562 first argument, a start vector "$x0_vector", a vector "$b_vector"
1563 (which is the vector "b" in your equation system ""A * x = b""), in
1564 the case of the "Relaxation Method" ("RM"), a real number "$weight"
1565 best between zero and two, and finally an error limit (real number)
1566 "$epsilon".
1567 (Note that the weight "$weight" used by the "Relaxation Method"
1569 ("RM") is NOT checked to lie within any reasonable range!)
1570 The three methods first test the first two conditions of the three
1572 conditions listed above and return an empty list if these condi‐
1573 tions are not fulfilled.
1574 Therefore, you should always test their return value using some
1576 code like:
1577 if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
1579 {
1580 # do something with the solution...
1581 }
1582 else
1583 {
1584 # do something with the fact that there is no solution...
1585 }
1586 Otherwise, they iterate until "abs( Phi(x) - x ) < epsilon".
1588 (Beware that theoretically, infinite loops might result if the
1590 starting vector is too far "off" the solution! In practice, this
1591 shouldn't be a problem. Anyway, you can always press <ctrl-C> if
1592 you think that the iteration takes too long!)
1593 The difference between the three methods is the following:
1595 In the "Global Step Method" ("GSM"), the new vector ""x(t+1)""
1597 (called "y" here) is calculated from the vector "x(t)" (called "x"
1598 here) according to the formula:
1599 y[i] =
1601 ( b[i]
1602 - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
1603 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
1604 ) / a[i,i]
1605 In the "Single Step Method" ("SSM"), the components of the vector
1607 ""x(t+1)"" which have already been calculated are used to calculate
1608 the remaining components, i.e.
1609 y[i] =
1611 ( b[i]
1612 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
1613 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
1614 ) / a[i,i]
1615 In the "Relaxation method" ("RM"), the components of the vector
1617 ""x(t+1)"" are calculated by "mixing" old and new value (like cold
1618 and hot water), and the weight "$weight" determines the "aperture"
1619 of both the "hot water tap" as well as of the "cold water tap",
1620 according to the formula:
1621 y[i] =
1623 ( b[i]
1624 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
1625 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
1626 ) / a[i,i]
1627 y[i] = weight * y[i] + (1 - weight) * x[i]
1628 Note that the weight "$weight" should be greater than zero and less
1630 than two (!).
1631 The three methods are supposed to be of different efficiency.
1633 Experiment!
1634 Remember that in most cases, it is probably advantageous to first
1636 "normalize()" your equation system prior to solving it!
1637
1639 SYNOPSIS
1640 · Unary operators:
1642 ""-"", ""~"", ""abs"", "test", ""!"", '""'
1644 · Binary operators:
1646 "".""
1648 Binary (arithmetic) operators:
1650 ""+"", ""-"", ""*"", ""**"", ""+="", ""-="", ""*="",
1652 ""/="",""**=""
1653 · Binary (relational) operators:
1655 ""=="", ""!="", ""<"", ""<="", "">"", "">=""
1657 ""eq"", ""ne"", ""lt"", ""le"", ""gt"", ""ge""
1659 Note that the latter (""eq"", ""ne"", ... ) are just synonyms of
1661 the former (""=="", ""!="", ... ), defined for convenience only.
1662 DESCRIPTION
1664 '.' Concatenation
1666 Returns the two matrices concatenated side by side.
1668 Example: $c = $a . $b;
1670 For example, if
1672 $a=[ 1 2 ] $b=[ 5 6 ]
1674 [ 3 4 ] [ 7 8 ]
1675 then
1676 $c=[ 1 2 5 6 ]
1678 [ 3 4 7 8 ]
1679 Note that only matrices with the same number of rows may be
1681 concatenated.
1682 '-' Unary minus
1684 Returns the negative of the given matrix, i.e., the matrix
1686 with all elements multiplied with the factor "-1".
1687 Example:
1689 $matrix = -$matrix;
1691 '~' Transposition
1693 Returns the transposed of the given matrix.
1695 Examples:
1697 $temp = ~$vector * $vector;
1699 $length = sqrt( $temp->element(1,1) );
1700 if (~$matrix == $matrix) { # matrix is symmetric ... }
1702 abs Norm
1704 Returns the "one"-Norm of the given matrix.
1706 Example:
1708 $error = abs( $A * $x - $b );
1710 test Boolean test
1712 Tests wether there is at least one non-zero element in the
1714 matrix.
1715 Example:
1717 if ($xn_vector) { # result of iteration is not zero ... }
1719 '!' Negated boolean test
1721 Tests wether the matrix contains only zero's.
1723 Examples:
1725 if (! $b_vector) { # heterogenous equation system ... }
1727 else { # homogenous equation system ... }
1728 unless ($x_vector) { # not the null-vector! }
1730 '""""'
1732 "Stringify" operator
1733 Converts the given matrix into a string.
1735 Uses scientific representation to keep precision loss to a
1737 minimum in case you want to read this string back in again
1738 later with "new_from_string()".
1739 Uses a 13-digit mantissa and a 20-character field for each
1741 element so that lines will wrap nicely on an 80-column screen.
1742 Examples:
1744 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
1746 [ 1 0 ]
1747 [ 0 -1 ]
1748 MATRIX
1749 print "$matrix";
1750 [ 1.000000000000E+00 0.000000000000E+00 ]
1752 [ 0.000000000000E+00 -1.000000000000E+00 ]
1753 $string = "$matrix";
1755 $test = Math::MatrixReal->new_from_string($string);
1756 if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
1757 '+' Addition
1759 Returns the sum of the two given matrices.
1761 Examples:
1763 $matrix_S = $matrix_A + $matrix_B;
1765 $matrix_A += $matrix_B;
1767 '-' Subtraction
1769 Returns the difference of the two given matrices.
1771 Examples:
1773 $matrix_D = $matrix_A - $matrix_B;
1775 $matrix_A -= $matrix_B;
1777 Note that this is the same as:
1779 $matrix_S = $matrix_A + -$matrix_B;
1781 $matrix_A += -$matrix_B;
1783 (The latter are less efficient, though)
1785 '*' Multiplication
1787 Returns the matrix product of the two given matrices or the
1789 product of the given matrix and scalar factor.
1790 Examples:
1792 $matrix_P = $matrix_A * $matrix_B;
1794 $matrix_A *= $matrix_B;
1796 $vector_b = $matrix_A * $vector_x;
1798 $matrix_B = -1 * $matrix_A;
1800 $matrix_B = $matrix_A * -1;
1802 $matrix_A *= -1;
1804 '/' Division
1806 Currently a shortcut for doing $a * $b ** -1 is $a / $b, which
1808 works for square matrices. One can also use 1/$a .
1809 '**' Exponentiation
1811 Returns the matrix raised to an integer power. If 0 is passed,
1813 the identity matrix is returned. If a negative integer is
1814 passed, it computes the inverse (if it exists) and then raised
1815 the inverse to the absolute value of the integer. The matrix
1816 must be quadratic.
1817 Examples:
1819 $matrix2 = $matrix ** 2;
1821 $matrix **= 2;
1823 $inv2 = $matrix ** -2;
1825 $ident = $matrix ** 0;
1827 '==' Equality
1829 Tests two matrices for equality.
1831 Example:
1833 if ( $A * $x == $b ) { print "EUREKA!\n"; }
1835 Note that in most cases, due to numerical errors (due to the
1837 finite precision of computer arithmetics), it is a bad idea to
1838 compare two matrices or vectors this way.
1839 Better use the norm of the difference of the two matrices you
1841 want to compare and compare that norm with a small number,
1842 like this:
1843 if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
1845 '!=' Inequality
1847 Tests two matrices for inequality.
1849 Example:
1851 while ($x0_vector != $xn_vector) { # proceed with iteration ... }
1853 (Stops when the iteration becomes stationary)
1855 Note that (just like with the '==' operator), it is usually a
1857 bad idea to compare matrices or vectors this way. Compare the
1858 norm of the difference of the two matrices with a small number
1859 instead.
1860 '<' Less than
1862 Examples:
1864 if ( $matrix1 < $matrix2 ) { # ... }
1866 if ( $vector < $epsilon ) { # ... }
1868 if ( 1E-12 < $vector ) { # ... }
1870 if ( $A * $x - $b < 1E-12 ) { # ... }
1872 These are just shortcuts for saying:
1874 if ( abs($matrix1) < abs($matrix2) ) { # ... }
1876 if ( abs($vector) < abs($epsilon) ) { # ... }
1878 if ( abs(1E-12) < abs($vector) ) { # ... }
1880 if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
1882 Uses the "one"-norm for matrices and Perl's built-in "abs()"
1884 for scalars.
1885 '<=' Less than or equal
1887 As with the '<' operator, this is just a shortcut for the same
1889 expression with "abs()" around all arguments.
1890 Example:
1892 if ( $A * $x - $b <= 1E-12 ) { # ... }
1894 which in fact is the same as:
1896 if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
1898 Uses the "one"-norm for matrices and Perl's built-in "abs()"
1900 for scalars.
1901 '>' Greater than
1903 As with the '<' and '<=' operator, this
1905 if ( $xn - $x0 > 1E-12 ) { # ... }
1907 is just a shortcut for:
1909 if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
1911 Uses the "one"-norm for matrices and Perl's built-in "abs()"
1913 for scalars.
1914 '>=' Greater than or equal
1916 As with the '<', '<=' and '>' operator, the following
1918 if ( $LR >= $A ) { # ... }
1920 is simply a shortcut for:
1922 if ( abs($LR) >= abs($A) ) { # ... }
1924 Uses the "one"-norm for matrices and Perl's built-in "abs()"
1926 for scalars.
1927
1929 Math::VectorReal, Math::PARI, Math::MatrixBool, DFA::Kleene,
1930 Math::Kleene, Set::IntegerRange, Set::IntegerFast .
1931
1933 This man page documents Math::MatrixReal version 2.02.
1934 The latest version can always be found at
1936 http://leto.net/code/Math-MatrixReal/
1937
1939 Steffen Beyer <sb@engelschall.com>, Rodolphe Ortalo <ortalo@laas.fr>,
1940 Jonathan Leto <jonathan@leto.net>.
1941 Currently maintained by Jonathan Leto, send all bugs/patches to me.
1943
1945 Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for
1946 Algebra and Linear Algebra at the university (RWTH Aachen, Germany),
1947 and to Prof. Esser and his assistant, Mr. Jarausch, for their fascinat‐
1948 ing lectures in Numerical Analysis!
1949
1951 Copyright (c) 1996-2002 by Steffen Beyer, Rodolphe Ortalo, Jonathan
1952 Leto. All rights reserved.
1953
1955 This package is free software; you can redistribute it and/or modify it
1956 under the same terms as Perl itself.
1957perl v5.8.8 2008-02-22 Math::MatrixReal(3)