1Math::MatrixReal(3) User Contributed Perl Documentation Math::MatrixReal(3)
2
6 Math::MatrixReal - Matrix of Reals
7 Implements the data type "matrix of real numbers" (and consequently
9 also "vector of real numbers").
10
12 my $a = Math::MatrixReal->new_random(5, 5);
13 my $b = $a->new_random(10, 30, { symmetric=>1, bounded_by=>[-1,1] });
15 my $c = $b * $a ** 3;
17 my $d = $b->new_from_rows( [ [ 5, 3 ,4], [3, 4, 5], [ 2, 4, 1 ] ] );
19 print $a;
21 my $row = ($a * $b)->row(3);
23 my $col = (5*$c)->col(2);
25 my $transpose = ~$c;
27 my $transpose = $c->transpose;
29 my $inverse = $a->inverse;
31 my $inverse = 1/$a;
33 my $inverse = $a ** -1;
35 my $determinant= $a->det;
37 · $matrix->display_precision($integer)
39 Sets the default precision when matrices are printed or
41 stringified. $matrix->display_precision(0) will only show the
42 integer part of all the entries of $matrix and
43 $matrix->display_precision() will return to the default scientific
44 display notation. This method does not effect the precision of the
45 calculations.
46
48 Constructor Methods And Such
49 · use Math::MatrixReal;
50 Makes the methods and overloaded operators of this module available
52 to your program.
53 · $new_matrix = new Math::MatrixReal($rows,$columns);
55 The matrix object constructor method. A new matrix of size $rows by
57 $columns will be created, with the value 0.0 for all elements.
58 Note that this method is implicitly called by many of the other
60 methods in this module.
61 · $new_matrix = $some_matrix->new($rows,$columns);
63 Another way of calling the matrix object constructor method.
65 Matrix $some_matrix is not changed by this in any way.
67 · $new_matrix = $matrix->new_from_cols( [
69 $column_vector|$array_ref|$string, ... ] )
70 Creates a new matrix given a reference to an array of any of the
72 following:
73 · column vectors ( n by 1 Math::MatrixReal matrices )
75 · references to arrays
77 · strings properly formatted to create a column with
79 Math::MatrixReal's new_from_string command
80 You may mix and match these as you wish. However, all must be of
82 the same dimension--no padding happens automatically. Example:
83 my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
85 print $matrix;
86 will print
88 [ 1.000000000000E+00 3.000000000000E+00 ]
90 [ 2.000000000000E+00 4.000000000000E+00 ]
91 · new_from_rows( [ $row_vector|$array_ref|$string, ... ] )
93 Creates a new matrix given a reference to an array of any of the
95 following:
96 · row vectors ( 1 by n Math::MatrixReal matrices )
98 · references to arrays
100 · strings properly formatted to create a row with
102 Math::MatrixReal's new_from_string command
103 You may mix and match these as you wish. However, all must be of
105 the same dimension--no padding happens automatically. Example:
106 my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
108 print $matrix;
109 will print
111 [ 1.000000000000E+00 2.000000000000E+00 ]
113 [ 3.000000000000E+00 4.000000000000E+00 ]
114 · $new_matrix = Math::MatrixReal->new_random($rows, $cols, %options
116 );
117 This method allows you to create a random matrix with various
119 properties controlled by the %options matrix, which is optional.
120 The default values of the %options matrix are { integer => 0,
121 symmetric => 0, tridiagonal => 0, diagonal => 0, bounded_by =>
122 [0,10] } .
123 Example:
125 $matrix = Math::MatrixReal->new_random(4, { diagonal => 1, integer => 1 } );
127 print $matrix;
128 will print a 4x4 random diagonal matrix with integer entries
130 between zero and ten, something like
131 [ 5.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
133 [ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
134 [ 0.000000000000E+00 0.000000000000E+00 1.000000000000E+00 0.000000000000E+00 ]
135 [ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 8.000000000000E+00 ]
136 · $new_matrix = Math::MatrixReal->new_diag( $array_ref );
138 This method allows you to create a diagonal matrix by only
140 specifying the diagonal elements. Example:
141 $matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
143 print $matrix;
144 will print
146 [ 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
148 [ 0.000000000000E+00 2.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
149 [ 0.000000000000E+00 0.000000000000E+00 3.000000000000E+00 0.000000000000E+00 ]
150 [ 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 4.000000000000E+00 ]
151 · $new_matrix = Math::MatrixReal->new_tridiag( $lower, $diag, $upper
153 );
154 This method allows you to create a tridiagonal matrix by only
156 specifying the lower diagonal, diagonal and upper diagonal,
157 respectively.
158 $matrix = Math::MatrixReal->new_tridiag( [ 6, 4, 2 ], [1,2,3,4], [1, 8, 9] );
160 print $matrix;
161 will print
163 [ 1.000000000000E+00 1.000000000000E+00 0.000000000000E+00 0.000000000000E+00 ]
165 [ 6.000000000000E+00 2.000000000000E+00 8.000000000000E+00 0.000000000000E+00 ]
166 [ 0.000000000000E+00 4.000000000000E+00 3.000000000000E+00 9.000000000000E+00 ]
167 [ 0.000000000000E+00 0.000000000000E+00 2.000000000000E+00 4.000000000000E+00 ]
168 · $new_matrix = Math::MatrixReal->new_from_string($string);
170 This method allows you to read in a matrix from a string (for
172 instance, from the keyboard, from a file or from your code).
173 The syntax is simple: each row must start with ""[ "" and end with
175 "" ]\n"" (""\n"" being the newline character and "" "" a space or
176 tab) and contain one or more numbers, all separated from each other
177 by spaces or tabs.
178 Additional spaces or tabs can be added at will, but no comments.
180 Examples:
182 $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
184 $matrix = Math::MatrixReal->new_from_string($string);
185 print "$matrix";
186 By the way, this prints
188 [ 1.000000000000E+00 2.000000000000E+00 3.000000000000E+00 ]
190 [ 2.000000000000E+00 2.000000000000E+00 -1.000000000000E+00 ]
191 [ 1.000000000000E+00 1.000000000000E+00 1.000000000000E+00 ]
192 But you can also do this in a much more comfortable way using the
194 shell-like "here-document" syntax:
195 $matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
197 [ 1 0 0 0 0 0 1 ]
198 [ 0 1 0 0 0 0 0 ]
199 [ 0 0 1 0 0 0 0 ]
200 [ 0 0 0 1 0 0 0 ]
201 [ 0 0 0 0 1 0 0 ]
202 [ 0 0 0 0 0 1 0 ]
203 [ 1 0 0 0 0 0 -1 ]
204 MATRIX
205 You can even use variables in the matrix:
207 $c1 = 2 / 3;
209 $c2 = -2 / 5;
210 $c3 = 26 / 9;
211 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
213 [ 3 2 0 ]
215 [ 0 3 2 ]
216 [ $c1 $c2 $c3 ]
217 MATRIX
219 (Remember that you may use spaces and tabs to format the matrix to
221 your taste)
222 Note that this method uses exactly the same representation for a
224 matrix as the "stringify" operator "": this means that you can
225 convert any matrix into a string with "$string = "$matrix";" and
226 read it back in later (for instance from a file!).
227 Note however that you may suffer a precision loss in this process
229 because only 13 digits are supported in the mantissa when printed!!
230 If the string you supply (or someone else supplies) does not obey
232 the syntax mentioned above, an exception is raised, which can be
233 caught by "eval" as follows:
234 print "Please enter your matrix (in one line): ";
236 $string = <STDIN>;
237 $string =~ s/\\n/\n/g;
238 eval { $matrix = Math::MatrixReal->new_from_string($string); };
239 if ($@)
240 {
241 print "$@";
242 # ...
243 # (error handling)
244 }
245 else
246 {
247 # continue...
248 }
249 or as follows:
251 eval { $matrix = Math::MatrixReal->new_from_string(<<"MATRIX"); };
253 [ 3 2 0 ]
254 [ 0 3 2 ]
255 [ $c1 $c2 $c3 ]
256 MATRIX
257 if ($@)
258 # ...
259 Actually, the method shown above for reading a matrix from the
261 keyboard is a little awkward, since you have to enter a lot of
262 "\n"'s for the newlines.
263 A better way is shown in this piece of code:
265 while (1)
267 {
268 print "\nPlease enter your matrix ";
269 print "(multiple lines, <ctrl-D> = done):\n";
270 eval { $new_matrix =
271 Math::MatrixReal->new_from_string(join('',<STDIN>)); };
272 if ($@)
273 {
274 $@ =~ s/\s+at\b.*?$//;
275 print "${@}Please try again.\n";
276 }
277 else { last; }
278 }
279 Possible error messages of the "new_from_string()" method are:
281 Math::MatrixReal::new_from_string(): syntax error in input string
283 Math::MatrixReal::new_from_string(): empty input string
284 If the input string has rows with varying numbers of columns, the
286 following warning will be printed to STDERR:
287 Math::MatrixReal::new_from_string(): missing elements will be set to zero!
289 If everything is okay, the method returns an object reference to
291 the (newly allocated) matrix containing the elements you specified.
292 · $new_matrix = $some_matrix->shadow();
294 Returns an object reference to a NEW but EMPTY matrix (filled with
296 zero's) of the SAME SIZE as matrix "$some_matrix".
297 Matrix "$some_matrix" is not changed by this in any way.
299 · $matrix1->copy($matrix2);
301 Copies the contents of matrix "$matrix2" to an ALREADY EXISTING
303 matrix "$matrix1" (which must have the same size as matrix
304 "$matrix2"!).
305 Matrix "$matrix2" is not changed by this in any way.
307 · $twin_matrix = $some_matrix->clone();
309 Returns an object reference to a NEW matrix of the SAME SIZE as
311 matrix "$some_matrix". The contents of matrix "$some_matrix" have
312 ALREADY BEEN COPIED to the new matrix "$twin_matrix". This is the
313 method that the operator "=" is overloaded to when you type "$a =
314 $b", when $a and $b are matrices.
315 Matrix "$some_matrix" is not changed by this in any way.
317 Matrix Row, Column and Element operations
319 · $row_vector = $matrix->row($row);
320 This is a projection method which returns an object reference to a
322 NEW matrix (which in fact is a (row) vector since it has only one
323 row) to which row number "$row" of matrix "$matrix" has already
324 been copied.
325 Matrix "$matrix" is not changed by this in any way.
327 · $column_vector = $matrix->column($column);
329 This is a projection method which returns an object reference to a
331 NEW matrix (which in fact is a (column) vector since it has only
332 one column) to which column number "$column" of matrix "$matrix"
333 has already been copied.
334 Matrix "$matrix" is not changed by this in any way.
336 · $matrix->assign($row,$column,$value);
338 Explicitly assigns a value "$value" to a single element of the
340 matrix "$matrix", located in row "$row" and column "$column",
341 thereby replacing the value previously stored there.
342 · $value = $matrix->element($row,$column);
344 Returns the value of a specific element of the matrix "$matrix",
346 located in row "$row" and column "$column".
347 · $new_matrix = $matrix->each( \&function );
349 Creates a new matrix by evaluating a code reference on each element
351 of the given matrix. The function is passed the element, the row
352 index and the column index, in that order. The value the function
353 returns ( or the value of the last executed statement ) is the
354 value given to the corresponding element in $new_matrix.
355 Example:
357 # add 1 to every element in the matrix
359 $matrix = $matrix->each ( sub { (shift) + 1 } );
360 Example:
362 my $cofactor = $matrix->each( sub { my(undef,$i,$j) = @_;
364 ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
365 : -1*$matrix->minor($i,$j)->det();
366 } );
367 This code needs some explanation. For each element of $matrix, it
369 throws away the actual value and stores the row and column indexes
370 in $i and $j. Then it sets element [$i,$j] in $cofactor to the
371 determinant of "$matrix->minor($i,$j)" if it is an "even" element,
372 or "-1*$matrix->minor($i,$j)" if it is an "odd" element.
373 · $new_matrix = $matrix->each_diag( \&function );
375 Creates a new matrix by evaluating a code reference on each
377 diagonal element of the given matrix. The function is passed the
378 element, the row index and the column index, in that order. The
379 value the function returns ( or the value of the last executed
380 statement ) is the value given to the corresponding element in
381 $new_matrix.
382 · $matrix->swap_col( $col1, $col2 );
384 This method takes two one-based column numbers and swaps the values
386 of each element in each column. "$matrix->swap_col(2,3)" would
387 replace column 2 in $matrix with column 3, and replace column 3
388 with column 2.
389 · $matrix->swap_row( $row1, $row2 );
391 This method takes two one-based row numbers and swaps the values of
393 each element in each row. "$matrix->swap_row(2,3)" would replace
394 row 2 in $matrix with row 3, and replace row 3 with row 2.
395 · $matrix->assign_row( $row_number , $new_row_vector );
397 This method takes a one-based row number and assigns row
399 $row_number of $matrix with $new_row_vector and returns the
400 resulting matrix. "$matrix->assign_row(5, $x)" would replace row 2
401 in $matrix with the row vector $x.
402 Matrix Operations
404 · $det = $matrix->det();
405 Returns the determinant of the matrix, without going through the
407 rigamarole of computing a LR decomposition. This method should be
408 much faster than LR decomposition if the matrix is diagonal or
409 triangular. Otherwise, it is just a wrapper for
410 "$matrix->decompose_LR->det_LR". If the determinant is zero, there
411 is no inverse and vice-versa. Only quadratic matrices have
412 determinants.
413 · "$inverse = $matrix->inverse();"
415 Returns the inverse of a matrix, without going through the
417 rigamarole of computing a LR decomposition. If no inverse exists,
418 undef is returned and an error is printed via "carp()". This is
419 nothing but a wrapper for "$matrix->decompose_LR->invert_LR".
420 · "($rows,$columns) = $matrix->dim();"
422 Returns a list of two items, representing the number of rows and
424 columns the given matrix "$matrix" contains.
425 · "$norm_one = $matrix->norm_one();"
427 Returns the "one"-norm of the given matrix "$matrix".
429 The "one"-norm is defined as follows:
431 For each column, the sum of the absolute values of the elements in
433 the different rows of that column is calculated. Finally, the
434 maximum of these sums is returned.
435 Note that the "one"-norm and the "maximum"-norm are mathematically
437 equivalent, although for the same matrix they usually yield a
438 different value.
439 Therefore, you should only compare values that have been calculated
441 using the same norm!
442 Throughout this package, the "one"-norm is (arbitrarily) used for
444 all comparisons, for the sake of uniformity and comparability,
445 except for the iterative methods "solve_GSM()", "solve_SSM()" and
446 "solve_RM()" which use either norm depending on the matrix itself.
447 · "$norm_max = $matrix->norm_max();"
449 Returns the "maximum"-norm of the given matrix $matrix.
451 The "maximum"-norm is defined as follows:
453 For each row, the sum of the absolute values of the elements in the
455 different columns of that row is calculated. Finally, the maximum
456 of these sums is returned.
457 Note that the "maximum"-norm and the "one"-norm are mathematically
459 equivalent, although for the same matrix they usually yield a
460 different value.
461 Therefore, you should only compare values that have been calculated
463 using the same norm!
464 Throughout this package, the "one"-norm is (arbitrarily) used for
466 all comparisons, for the sake of uniformity and comparability,
467 except for the iterative methods "solve_GSM()", "solve_SSM()" and
468 "solve_RM()" which use either norm depending on the matrix itself.
469 · "$norm_sum = $matrix->norm_sum();"
471 This is a very simple norm which is defined as the sum of the
473 absolute values of every element.
474 · $p_norm = $matrix->norm_p($n);>
476 This function returns the "p-norm" of a vector. The argument $n
478 must be a number greater than or equal to 1 or the string "Inf".
479 The p-norm is defined as (sum(x_i^p))^(1/p). In words, it raised
480 each element to the p-th power, adds them up, and then takes the
481 p-th root of that number. If the string "Inf" is passed, the
482 "infinity-norm" is computed, which is really the limit of the
483 p-norm as p goes to infinity. It is defined as the maximum element
484 of the vector. Also, note that the familiar Euclidean distance
485 between two vectors is just a special case of a p-norm, when p is
486 equal to 2.
487 Example:
489 $a = Math::MatrixReal->new_from_cols([[1,2,3]]);
490 $p1 = $a->norm_p(1);
491 $p2 = $a->norm_p(2);
492 $p3 = $a->norm_p(3);
493 $pinf = $a->norm_p("Inf");
494 print "(1,2,3,Inf) norm:\n$p1\n$p2\n$p3\n$pinf\n";
496 $i1 = $a->new_from_rows([[1,0]]);
498 $i2 = $a->new_from_rows([[0,1]]);
499 # this should be sqrt(2) since it is the same as the
501 # hypotenuse of a 1 by 1 right triangle
502 $dist = ($i1-$i2)->norm_p(2);
504 print "Distance is $dist, which should be " . sqrt(2) . "\n";
505 Output:
507 (1,2,3,Inf) norm:
509 6
510 3.74165738677394139
511 3.30192724889462668
512 3
513 Distance is 1.41421356237309505, which should be 1.41421356237309505
515 · $frob_norm = "$matrix->norm_frobenius();"
517 This norm is similar to that of a p-norm where p is 2, except it
519 acts on a matrix, not a vector. Each element of the matrix is
520 squared, this is added up, and then a square root is taken.
521 · "$matrix->spectral_radius();"
523 Returns the maximum value of the absolute value of all eigenvalues.
525 Currently this computes all eigenvalues, then sifts through them to
526 find the largest in absolute value. Needless to say, this is very
527 inefficient, and in the future an algorithm that computes only the
528 largest eigenvalue may be implemented.
529 · "$matrix1->transpose($matrix2);"
531 Calculates the transposed matrix of matrix $matrix2 and stores the
533 result in matrix "$matrix1" (which must already exist and have the
534 same size as matrix "$matrix2"!).
535 This operation can also be carried out "in-place", i.e., input and
537 output matrix may be identical.
538 Transposition is a symmetry operation: imagine you rotate the
540 matrix along the axis of its main diagonal (going through elements
541 (1,1), (2,2), (3,3) and so on) by 180 degrees.
542 Another way of looking at it is to say that rows and columns are
544 swapped. In fact the contents of element "(i,j)" are swapped with
545 those of element "(j,i)".
546 Note that (especially for vectors) it makes a big difference if you
548 have a row vector, like this:
549 [ -1 0 1 ]
551 or a column vector, like this:
553 [ -1 ]
555 [ 0 ]
556 [ 1 ]
557 the one vector being the transposed of the other!
559 This is especially true for the matrix product of two vectors:
561 [ -1 ]
563 [ -1 0 1 ] * [ 0 ] = [ 2 ] , whereas
564 [ 1 ]
565 * [ -1 0 1 ]
567 [ -1 ] [ 1 0 -1 ]
568 [ 0 ] * [ -1 0 1 ] = [ -1 ] [ 1 0 -1 ] = [ 0 0 0 ]
569 [ 1 ] [ 0 ] [ 0 0 0 ] [ -1 0 1 ]
570 [ 1 ] [ -1 0 1 ]
571 So be careful about what you really mean!
573 Hint: throughout this module, whenever a vector is explicitly
575 required for input, a COLUMN vector is expected!
576 · "$trace = $matrix->trace();"
578 This returns the trace of the matrix, which is defined as the sum
580 of the diagonal elements. The matrix must be quadratic.
581 · "$minor = $matrix->minor($row,$col);"
583 Returns the minor matrix corresponding to $row and $col. $matrix
585 must be quadratic. If $matrix is n rows by n cols, the minor of
586 $row and $col will be an (n-1) by (n-1) matrix. The minor is
587 defined as crossing out the row and the col specified and returning
588 the remaining rows and columns as a matrix. This method is used by
589 "cofactor()".
590 · "$cofactor = $matrix->cofactor();"
592 The cofactor matrix is constructed as follows:
594 For each element, cross out the row and column that it sits in.
596 Now, take the determinant of the matrix that is left in the other
597 rows and columns. Multiply the determinant by (-1)^(i+j), where i
598 is the row index, and j is the column index. Replace the given
599 element with this value.
600 The cofactor matrix can be used to find the inverse of the matrix.
602 One formula for the inverse of a matrix is the cofactor matrix
603 transposed divided by the original determinant of the matrix.
604 The following two inverses should be exactly the same:
606 my $inverse1 = $matrix->inverse;
608 my $inverse2 = ~($matrix->cofactor)->each( sub { (shift)/$matrix->det() } );
609 Caveat: Although the cofactor matrix is simple algorithm to compute
611 the inverse of a matrix, and can be used with pencil and paper for
612 small matrices, it is comically slower than the native "inverse()"
613 function. Here is a small benchmark:
614 # $matrix1 is 15x15
616 $det = $matrix1->det;
617 timethese( 10,
618 {'inverse' => sub { $matrix1->inverse(); },
619 'cofactor' => sub { (~$matrix1->cofactor)->each ( sub { (shift)/$det; } ) }
620 } );
621 Benchmark: timing 10 iterations of LR, cofactor, inverse...
624 inverse: 1 wallclock secs ( 0.56 usr + 0.00 sys = 0.56 CPU) @ 17.86/s (n=10)
625 cofactor: 36 wallclock secs (36.62 usr + 0.01 sys = 36.63 CPU) @ 0.27/s (n=10)
626 · "$adjoint = $matrix->adjoint();"
628 The adjoint is just the transpose of the cofactor matrix. This
630 method is just an alias for " ~($matrix->cofactor)".
631 Arithmetic Operations
633 · "$matrix1->add($matrix2,$matrix3);"
634 Calculates the sum of matrix "$matrix2" and matrix "$matrix3" and
636 stores the result in matrix "$matrix1" (which must already exist
637 and have the same size as matrix "$matrix2" and matrix
638 "$matrix3"!).
639 This operation can also be carried out "in-place", i.e., the output
641 and one (or both) of the input matrices may be identical.
642 · "$matrix1->subtract($matrix2,$matrix3);"
644 Calculates the difference of matrix "$matrix2" minus matrix
646 "$matrix3" and stores the result in matrix "$matrix1" (which must
647 already exist and have the same size as matrix "$matrix2" and
648 matrix "$matrix3"!).
649 This operation can also be carried out "in-place", i.e., the output
651 and one (or both) of the input matrices may be identical.
652 Note that this operation is the same as
654 "$matrix1->add($matrix2,-$matrix3);", although the latter is a
655 little less efficient.
656 · "$matrix1->multiply_scalar($matrix2,$scalar);"
658 Calculates the product of matrix "$matrix2" and the number
660 "$scalar" (i.e., multiplies each element of matrix "$matrix2" with
661 the factor "$scalar") and stores the result in matrix "$matrix1"
662 (which must already exist and have the same size as matrix
663 "$matrix2"!).
664 This operation can also be carried out "in-place", i.e., input and
666 output matrix may be identical.
667 · "$product_matrix = $matrix1->multiply($matrix2);"
669 Calculates the product of matrix "$matrix1" and matrix "$matrix2"
671 and returns an object reference to a new matrix "$product_matrix"
672 in which the result of this operation has been stored.
673 Note that the dimensions of the two matrices "$matrix1" and
675 "$matrix2" (i.e., their numbers of rows and columns) must harmonize
676 in the following way (example):
677 [ 2 2 ]
679 [ 2 2 ]
680 [ 2 2 ]
681 [ 1 1 1 ] [ * * ]
683 [ 1 1 1 ] [ * * ]
684 [ 1 1 1 ] [ * * ]
685 [ 1 1 1 ] [ * * ]
686 I.e., the number of columns of matrix "$matrix1" has to be the same
688 as the number of rows of matrix "$matrix2".
689 The number of rows and columns of the resulting matrix
691 "$product_matrix" is determined by the number of rows of matrix
692 "$matrix1" and the number of columns of matrix "$matrix2",
693 respectively.
694 · "$matrix1->negate($matrix2);"
696 Calculates the negative of matrix "$matrix2" (i.e., multiplies all
698 elements with "-1") and stores the result in matrix "$matrix1"
699 (which must already exist and have the same size as matrix
700 "$matrix2"!).
701 This operation can also be carried out "in-place", i.e., input and
703 output matrix may be identical.
704 · "$matrix_to_power = $matrix1->exponent($integer);"
706 Raises the matrix to the $integer power. Obviously, $integer must
708 be an integer. If it is zero, the identity matrix is returned. If a
709 negative integer is given, the inverse will be computed (if it
710 exists) and then raised the the absolute value of $integer. The
711 matrix must be quadratic.
712 Boolean Matrix Operations
714 · "$matrix->is_quadratic();"
715 Returns a boolean value indicating if the given matrix is quadratic
717 (also know as "square" or "n by n"). A matrix is quadratic if it
718 has the same number of rows as it does columns.
719 · "$matrix->is_square();"
721 This is an alias for "is_quadratic()".
723 · "$matrix->is_symmetric();"
725 Returns a boolean value indicating if the given matrix is
727 symmetric. By definition, a matrix is symmetric if and only if
728 (M[i,j]=M[j,i]). This is equivalent to "($matrix == ~$matrix)" but
729 without memory allocation. Only quadratic matrices can be
730 symmetric.
731 Notes: A symmetric matrix always has real eigenvalues/eigenvectors.
733 A matrix plus its transpose is always symmetric.
734 · "$matrix->is_skew_symmetric();"
736 Returns a boolean value indicating if the given matrix is skew
738 symmetric. By definition, a matrix is symmetric if and only if
739 (M[i,j]=-M[j,i]). This is equivalent to "($matrix == -(~$matrix))"
740 but without memory allocation. Only quadratic matrices can be skew
741 symmetric.
742 · "$matrix->is_diagonal();"
744 Returns a boolean value indicating if the given matrix is diagonal,
746 i.e. all of the nonzero elements are on the main diagonal. Only
747 quadratic matrices can be diagonal.
748 · "$matrix->is_tridiagonal();"
750 Returns a boolean value indicating if the given matrix is
752 tridiagonal, i.e. all of the nonzero elements are on the main
753 diagonal or the diagonals above and below the main diagonal. Only
754 quadratic matrices can be tridiagonal.
755 · "$matrix->is_upper_triangular();"
757 Returns a boolean value indicating if the given matrix is upper
759 triangular, i.e. all of the nonzero elements not on the main
760 diagonal are above it. Only quadratic matrices can be upper
761 triangular. Note: diagonal matrices are both upper and lower
762 triangular.
763 · "$matrix->is_lower_triangular();"
765 Returns a boolean value indicating if the given matrix is lower
767 triangular, i.e. all of the nonzero elements not on the main
768 diagonal are below it. Only quadratic matrices can be lower
769 triangular. Note: diagonal matrices are both upper and lower
770 triangular.
771 · "$matrix->is_orthogonal();"
773 Returns a boolean value indicating if the given matrix is
775 orthogonal. An orthogonal matrix is has the property that the
776 transpose equals the inverse of the matrix. Instead of computing
777 each and comparing them, this method multiplies the matrix by it's
778 transpose, and returns true if this turns out to be the identity
779 matrix, false otherwise. Only quadratic matrices can orthogonal.
780 · "$matrix->is_binary();"
782 Returns a boolean value indicating if the given matrix is binary.
784 A matrix is binary if it contains only zeroes or ones.
785 · "$matrix->is_gramian();"
787 Returns a boolean value indicating if the give matrix is Gramian.
789 A matrix $A is Gramian if and only if there exists a square matrix
790 $B such that "$A = ~$B*$B". This is equivalent to checking if $A is
791 symmetric and has all nonnegative eigenvalues, which is what
792 Math::MatrixReal uses to check for this property.
793 · "$matrix->is_LR();"
795 Returns a boolean value indicating if the matrix is an LR
797 decomposition matrix.
798 · "$matrix->is_positive();"
800 Returns a boolean value indicating if the matrix contains only
802 positive entries. Note that a zero entry is not positive and will
803 cause "is_positive()" to return false.
804 · "$matrix->is_negative();"
806 Returns a boolean value indicating if the matrix contains only
808 negative entries. Note that a zero entry is not negative and will
809 cause "is_negative()" to return false.
810 · "$matrix->is_periodic($k);"
812 Returns a boolean value indicating if the matrix is periodic with
814 period $k. This is true if "$matrix ** ($k+1) == $matrix". When
815 "$k == 1", this reduces down to the "is_idempotent()" function.
816 · "$matrix->is_idempotent();"
818 Returns a boolean value indicating if the matrix is idempotent,
820 which is defined as the square of the matrix being equal to the
821 original matrix, i.e "$matrix ** 2 == $matrix".
822 · "$matrix->is_row_vector();"
824 Returns a boolean value indicating if the matrix is a row vector.
826 A row vector is a matrix which is 1xn. Note that the 1x1 matrix is
827 both a row and column vector.
828 · "$matrix->is_col_vector();"
830 Returns a boolean value indicating if the matrix is a col vector.
832 A col vector is a matrix which is nx1. Note that the 1x1 matrix is
833 both a row and column vector.
834 Eigensystems
836 · "($l, $V) = $matrix->sym_diagonalize();"
837 This method performs the diagonalization of the quadratic symmetric
839 matrix M stored in $matrix. On output, l is a column vector
840 containing all the eigenvalues of M and V is an orthogonal matrix
841 which columns are the corresponding normalized eigenvectors. The
842 primary property of an eigenvalue l and an eigenvector x is of course
843 that: M * x = l * x.
844 The method uses a Householder reduction to tridiagonal form followed
846 by a QL algoritm with implicit shifts on this tridiagonal. (The
847 tridiagonal matrix is kept internally in a compact form in this
848 routine to save memory.) In fact, this routine wraps the
849 householder() and tri_diagonalize() methods described below when
850 their intermediate results are not desired. The overall algorithmic
851 complexity of this technique is O(N^3). According to several books,
852 the coefficient hidden by the 'O' is one of the best possible for
853 general (symmetric) matrixes.
854 · "($T, $Q) = $matrix->householder();"
856 This method performs the Householder algorithm which reduces the n by
858 n real symmetric matrix M contained in $matrix to tridiagonal form.
859 On output, T is a symmetric tridiagonal matrix (only diagonal and
860 off-diagonal elements are non-zero) and Q is an orthogonal matrix
861 performing the tranformation between M and T ("$M == $Q * $T * ~$Q").
862 · "($l, $V) = $T->tri_diagonalize([$Q]);"
864 This method diagonalizes the symmetric tridiagonal matrix T. On
866 output, $l and $V are similar to the output values described for
867 sym_diagonalize().
868 The optional argument $Q corresponds to an orthogonal transformation
870 matrix Q that should be used additionally during V (eigenvectors)
871 computation. It should be supplied if the desired eigenvectors
872 correspond to a more general symmetric matrix M previously reduced by
873 the householder() method, not a mere tridiagonal. If T is really a
874 tridiagonal matrix, Q can be omitted (it will be internally created
875 in fact as an identity matrix). The method uses a QL algorithm (with
876 implicit shifts).
877 · "$l = $matrix->sym_eigenvalues();"
879 This method computes the eigenvalues of the quadratic symmetric
881 matrix M stored in $matrix. On output, l is a column vector
882 containing all the eigenvalues of M. Eigenvectors are not computed
883 (on the contrary of "sym_diagonalize()") and this method is more
884 efficient (even though it uses a similar algorithm with two phases).
885 However, understand that the algorithmic complexity of this technique
886 is still also O(N^3). But the coefficient hidden by the 'O' is better
887 by a factor of..., well, see your benchmark, it's wiser.
888 This routine wraps the householder_tridiagonal() and
890 tri_eigenvalues() methods described below when the intermediate
891 tridiagonal matrix is not needed.
892 · "$T = $matrix->householder_tridiagonal();"
894 This method performs the Householder algorithm which reduces the n by
896 n real symmetric matrix M contained in $matrix to tridiagonal form.
897 On output, T is the obtained symmetric tridiagonal matrix (only
898 diagonal and off-diagonal elements are non-zero). The operation is
899 similar to the householder() method, but potentially a little more
900 efficient as the transformation matrix is not computed.
901 · $l = $T->tri_eigenvalues();
903 This method computesthe eigenvalues of the symmetric tridiagonal
905 matrix T. On output, $l is a vector containing the eigenvalues
906 (similar to "sym_eigenvalues()"). This method is much more efficient
907 than tri_diagonalize() when eigenvectors are not needed.
908 Miscellaneous
910 · $matrix->zero();
911 Assigns a zero to every element of the matrix "$matrix", i.e.,
913 erases all values previously stored there, thereby effectively
914 transforming the matrix into a "zero"-matrix or "null"-matrix, the
915 neutral element of the addition operation in a Ring.
916 (For instance the (quadratic) matrices with "n" rows and columns
918 and matrix addition and multiplication form a Ring. Most prominent
919 characteristic of a Ring is that multiplication is not commutative,
920 i.e., in general, ""matrix1 * matrix2"" is not the same as
921 ""matrix2 * matrix1""!)
922 · $matrix->one();
924 Assigns one's to the elements on the main diagonal (elements (1,1),
926 (2,2), (3,3) and so on) of matrix "$matrix" and zero's to all
927 others, thereby erasing all values previously stored there and
928 transforming the matrix into a "one"-matrix, the neutral element of
929 the multiplication operation in a Ring.
930 (If the matrix is quadratic (which this method doesn't require,
932 though), then multiplying this matrix with itself yields this same
933 matrix again, and multiplying it with some other matrix leaves that
934 other matrix unchanged!)
935 · "$latex_string = $matrix->as_latex( align=> "c", format => "%s",
937 name => "" );"
938 This function returns the matrix as a LaTeX string. It takes a hash
940 as an argument which is used to control the style of the output.
941 The hash element "align" may be "c","l" or "r", corresponding to
942 center, left and right, respectively. The "format" element is a
943 format string that is given to "sprintf" to control the style of
944 number format, such a floating point or scientific notation. The
945 "name" element can be used so that a LaTeX string of "$name = " is
946 prepended to the string.
947 Example:
949 my $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
951 print $a->as_latex( ( format => "%.2f", align => "l",name => "A" ) );
952 Output:
954 $A = $ $
955 \left( \begin{array}{ll}
956 1.23&1.00 \\
957 5.68&2.00 \\
958 9.10&3.00
959 \end{array} \right)
960 $
961 · "$yacas_string = $matrix->as_yacas( format => "%s", name => "",
963 semi => 0 );"
964 This function returns the matrix as a string that can be read by
966 Yacas. It takes a hash as an an argument which controls the style
967 of the output. The "format" element is a format string that is
968 given to "sprintf" to control the style of number format, such a
969 floating point or scientific notation. The "name" element can be
970 used so that "$name = " is prepended to the string. The <semi>
971 element can be set to 1 to that a semicolon is appended (so Matlab
972 does not print out the matrix.)
973 Example:
975 $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
977 print $a->as_yacas( ( format => "%.2f", align => "l",name => "A" ) );
978 Output:
980 A := {{1.23,1.00},{5.68,2.00},{9.10,3.00}}
982 · "$matlab_string = $matrix->as_matlab( format => "%s", name => "",
984 semi => 0 );"
985 This function returns the matrix as a string that can be read by
987 Matlab. It takes a hash as an an argument which controls the style
988 of the output. The "format" element is a format string that is
989 given to "sprintf" to control the style of number format, such a
990 floating point or scientific notation. The "name" element can be
991 used so that "$name = " is prepended to the string. The <semi>
992 element can be set to 1 to that a semicolon is appended (so Matlab
993 does not print out the matrix.)
994 Example:
996 my $a = Math::MatrixReal->new_from_rows([[ 1.234, 5.678, 9.1011],[1,2,3]] );
998 print $a->as_matlab( ( format => "%.3f", name => "A",semi => 1 ) );
999 Output:
1001 A = [ 1.234 5.678 9.101;
1002 1.000 2.000 3.000];
1003 · "$scilab_string = $matrix->as_scilab( format => "%s", name => "",
1005 semi => 0 );"
1006 This function is just an alias for "as_matlab()", since both Scilab
1008 and Matlab have the same matrix format.
1009 · "$minimum = Math::MatrixReal::min($number1,$number2);" "$minimum =
1011 Math::MatrixReal::min($matrix);" "$minimum = $matrix-"min;>
1012 Returns the minimum of the two numbers ""number1"" and ""number2""
1014 if called with two arguments, or returns the value of the smallest
1015 element of a matrix if called with one argument or as an object
1016 method.
1017 · "$maximum = Math::MatrixReal::max($number1,$number2);" "$maximum =
1019 Math::MatrixReal::max($number1,$number2);" "$maximum =
1020 Math::MatrixReal::max($matrix);" "$maximum = $matrix-"max;>
1021 Returns the maximum of the two numbers ""number1"" and ""number2""
1023 if called with two arguments, or returns the value of the largest
1024 element of a matrix if called with one arguemnt or as on object
1025 method.
1026 · "$minimal_cost_matrix = $cost_matrix->kleene();"
1028 Copies the matrix "$cost_matrix" (which has to be quadratic!) to a
1030 new matrix of the same size (i.e., "clones" the input matrix) and
1031 applies Kleene's algorithm to it.
1032 See Math::Kleene(3) for more details about this algorithm!
1034 The method returns an object reference to the new matrix.
1036 Matrix "$cost_matrix" is not changed by this method in any way.
1038 · "($norm_matrix,$norm_vector) = $matrix->normalize($vector);"
1040 This method is used to improve the numerical stability when solving
1042 linear equation systems.
1043 Suppose you have a matrix "A" and a vector "b" and you want to find
1045 out a vector "x" so that "A * x = b", i.e., the vector "x" which
1046 solves the equation system represented by the matrix "A" and the
1047 vector "b".
1048 Applying this method to the pair (A,b) yields a pair (A',b') where
1050 each row has been divided by (the absolute value of) the greatest
1051 coefficient appearing in that row. So this coefficient becomes
1052 equal to "1" (or "-1") in the new pair (A',b') (all others become
1053 smaller than one and greater than minus one).
1054 Note that this operation does not change the equation system itself
1056 because the same division is carried out on either side of the
1057 equation sign!
1058 The method requires a quadratic (!) matrix "$matrix" and a vector
1060 "$vector" for input (the vector must be a column vector with the
1061 same number of rows as the input matrix) and returns a list of two
1062 items which are object references to a new matrix and a new vector,
1063 in this order.
1064 The output matrix and vector are clones of the input matrix and
1066 vector to which the operation explained above has been applied.
1067 The input matrix and vector are not changed by this in any way.
1069 Example of how this method can affect the result of the methods to
1071 solve equation systems (explained immediately below following this
1072 method):
1073 Consider the following little program:
1075 #!perl -w
1077 use Math::MatrixReal qw(new_from_string);
1079 $A = Math::MatrixReal->new_from_string(<<"MATRIX");
1081 [ 1 2 3 ]
1082 [ 5 7 11 ]
1083 [ 23 19 13 ]
1084 MATRIX
1085 $b = Math::MatrixReal->new_from_string(<<"MATRIX");
1087 [ 0 ]
1088 [ 1 ]
1089 [ 29 ]
1090 MATRIX
1091 $LR = $A->decompose_LR();
1093 if (($dim,$x,$B) = $LR->solve_LR($b))
1094 {
1095 $test = $A * $x;
1096 print "x = \n$x";
1097 print "A * x = \n$test";
1098 }
1099 ($A_,$b_) = $A->normalize($b);
1101 $LR = $A_->decompose_LR();
1103 if (($dim,$x,$B) = $LR->solve_LR($b_))
1104 {
1105 $test = $A * $x;
1106 print "x = \n$x";
1107 print "A * x = \n$test";
1108 }
1109 This will print:
1111 x =
1113 [ 1.000000000000E+00 ]
1114 [ 1.000000000000E+00 ]
1115 [ -1.000000000000E+00 ]
1116 A * x =
1117 [ 4.440892098501E-16 ]
1118 [ 1.000000000000E+00 ]
1119 [ 2.900000000000E+01 ]
1120 x =
1121 [ 1.000000000000E+00 ]
1122 [ 1.000000000000E+00 ]
1123 [ -1.000000000000E+00 ]
1124 A * x =
1125 [ 0.000000000000E+00 ]
1126 [ 1.000000000000E+00 ]
1127 [ 2.900000000000E+01 ]
1128 You can see that in the second example (where "normalize()" has
1130 been used), the result is "better", i.e., more accurate!
1131 · "$LR_matrix = $matrix->decompose_LR();"
1133 This method is needed to solve linear equation systems.
1135 Suppose you have a matrix "A" and a vector "b" and you want to find
1137 out a vector "x" so that "A * x = b", i.e., the vector "x" which
1138 solves the equation system represented by the matrix "A" and the
1139 vector "b".
1140 You might also have a matrix "A" and a whole bunch of different
1142 vectors "b1".."bk" for which you need to find vectors "x1".."xk" so
1143 that "A * xi = bi", for "i=1..k".
1144 Using Gaussian transformations (multiplying a row or column with a
1146 factor, swapping two rows or two columns and adding a multiple of
1147 one row or column to another), it is possible to decompose any
1148 matrix "A" into two triangular matrices, called "L" and "R" (for
1149 "Left" and "Right").
1150 "L" has one's on the main diagonal (the elements (1,1), (2,2),
1152 (3,3) and so so), non-zero values to the left and below of the main
1153 diagonal and all zero's in the upper right half of the matrix.
1154 "R" has non-zero values on the main diagonal as well as to the
1156 right and above of the main diagonal and all zero's in the lower
1157 left half of the matrix, as follows:
1158 [ 1 0 0 0 0 ] [ x x x x x ]
1160 [ x 1 0 0 0 ] [ 0 x x x x ]
1161 L = [ x x 1 0 0 ] R = [ 0 0 x x x ]
1162 [ x x x 1 0 ] [ 0 0 0 x x ]
1163 [ x x x x 1 ] [ 0 0 0 0 x ]
1164 Note that ""L * R"" is equivalent to matrix "A" in the sense that
1166 "L * R * x = b <==> A * x = b" for all vectors "x", leaving out
1167 of account permutations of the rows and columns (these are taken
1168 care of "magically" by this module!) and numerical errors.
1169 Trick:
1171 Because we know that "L" has one's on its main diagonal, we can
1173 store both matrices together in the same array without information
1174 loss! I.e.,
1175 [ R R R R R ]
1177 [ L R R R R ]
1178 LR = [ L L R R R ]
1179 [ L L L R R ]
1180 [ L L L L R ]
1181 Beware, though, that "LR" and ""L * R"" are not the same!!!
1183 Note also that for the same reason, you cannot apply the method
1185 "normalize()" to an "LR" decomposition matrix. Trying to do so will
1186 yield meaningless rubbish!
1187 (You need to apply "normalize()" to each pair (Ai,bi) BEFORE
1189 decomposing the matrix "Ai'"!)
1190 Now what does all this help us in solving linear equation systems?
1192 It helps us because a triangular matrix is the next best thing that
1194 can happen to us besides a diagonal matrix (a matrix that has non-
1195 zero values only on its main diagonal - in which case the solution
1196 is trivial, simply divide ""b[i]"" by ""A[i,i]"" to get ""x[i]""!).
1197 To find the solution to our problem ""A * x = b"", we divide this
1199 problem in parts: instead of solving "A * x = b" directly, we first
1200 decompose "A" into "L" and "R" and then solve ""L * y = b"" and
1201 finally ""R * x = y"" (motto: divide and rule!).
1202 From the illustration above it is clear that solving ""L * y = b""
1204 and ""R * x = y"" is straightforward: we immediately know that
1205 "y[1] = b[1]". We then deduce swiftly that
1206 y[2] = b[2] - L[2,1] * y[1]
1208 (and we know ""y[1]"" by now!), that
1210 y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]
1212 and so on.
1214 Having effortlessly calculated the vector "y", we now proceed to
1216 calculate the vector "x" in a similar fashion: we see immediately
1217 that "x[n] = y[n] / R[n,n]". It follows that
1218 x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]
1220 and
1222 x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
1224 / R[n-2,n-2]
1225 and so on.
1227 You can see that - especially when you have many vectors "b1".."bk"
1229 for which you are searching solutions to "A * xi = bi" - this
1230 scheme is much more efficient than a straightforward, "brute force"
1231 approach.
1232 This method requires a quadratic matrix as its input matrix.
1234 If you don't have that many equations, fill up with zero's (i.e.,
1236 do nothing to fill the superfluous rows if it's a "fresh" matrix,
1237 i.e., a matrix that has been created with "new()" or "shadow()").
1238 The method returns an object reference to a new matrix containing
1240 the matrices "L" and "R".
1241 The input matrix is not changed by this method in any way.
1243 Note that you can "copy()" or "clone()" the result of this method
1245 without losing its "magical" properties (for instance concerning
1246 the hidden permutations of its rows and columns).
1247 However, as soon as you are applying any method that alters the
1249 contents of the matrix, its "magical" properties are stripped off,
1250 and the matrix immediately reverts to an "ordinary" matrix (with
1251 the values it just happens to contain at that moment, be they
1252 meaningful as an ordinary matrix or not!).
1253 · "($dimension,$x_vector,$base_matrix) =
1255 $LR_matrix""->""solve_LR($b_vector);"
1256 Use this method to actually solve an equation system.
1258 Matrix "$LR_matrix" must be a (quadratic) matrix returned by the
1260 method "decompose_LR()", the LR decomposition matrix of the matrix
1261 "A" of your equation system "A * x = b".
1262 The input vector "$b_vector" is the vector "b" in your equation
1264 system "A * x = b", which must be a column vector and have the same
1265 number of rows as the input matrix "$LR_matrix".
1266 The method returns a list of three items if a solution exists or an
1268 empty list otherwise (!).
1269 Therefore, you should always use this method like this:
1271 if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
1273 {
1274 # do something with the solution...
1275 }
1276 else
1277 {
1278 # do something with the fact that there is no solution...
1279 }
1280 The three items returned are: the dimension "$dimension" of the
1282 solution space (which is zero if only one solution exists, one if
1283 the solution is a straight line, two if the solution is a plane,
1284 and so on), the solution vector "$x_vector" (which is the vector
1285 "x" of your equation system "A * x = b") and a matrix
1286 "$base_matrix" representing a base of the solution space (a set of
1287 vectors which put up the solution space like the spokes of an
1288 umbrella).
1289 Only the first "$dimension" columns of this base matrix actually
1291 contain entries, the remaining columns are all zero.
1292 Now what is all this stuff with that "base" good for?
1294 The output vector "x" is ALWAYS a solution of your equation system
1296 "A * x = b".
1297 But also any vector "$vector"
1299 $vector = $x_vector->clone();
1301 $machine_infinity = 1E+99; # or something like that
1303 for ( $i = 1; $i <= $dimension; $i++ )
1305 {
1306 $vector += rand($machine_infinity) * $base_matrix->column($i);
1307 }
1308 is a solution to your problem "A * x = b", i.e., if "$A_matrix"
1310 contains your matrix "A", then
1311 print abs( $A_matrix * $vector - $b_vector ), "\n";
1313 should print a number around 1E-16 or so!
1315 By the way, note that you can actually calculate those vectors
1317 "$vector" a little more efficient as follows:
1318 $rand_vector = $x_vector->shadow();
1320 $machine_infinity = 1E+99; # or something like that
1322 for ( $i = 1; $i <= $dimension; $i++ )
1324 {
1325 $rand_vector->assign($i,1, rand($machine_infinity) );
1326 }
1327 $vector = $x_vector + ( $base_matrix * $rand_vector );
1329 Note that the input matrix and vector are not changed by this
1331 method in any way.
1332 · "$inverse_matrix = $LR_matrix->invert_LR();"
1334 Use this method to calculate the inverse of a given matrix
1336 "$LR_matrix", which must be a (quadratic) matrix returned by the
1337 method "decompose_LR()".
1338 The method returns an object reference to a new matrix of the same
1340 size as the input matrix containing the inverse of the matrix that
1341 you initially fed into "decompose_LR()" IF THE INVERSE EXISTS, or
1342 an empty list otherwise.
1343 Therefore, you should always use this method in the following way:
1345 if ( $inverse_matrix = $LR->invert_LR() )
1347 {
1348 # do something with the inverse matrix...
1349 }
1350 else
1351 {
1352 # do something with the fact that there is no inverse matrix...
1353 }
1354 Note that by definition (disregarding numerical errors), the
1356 product of the initial matrix and its inverse (or vice-versa) is
1357 always a matrix containing one's on the main diagonal (elements
1358 (1,1), (2,2), (3,3) and so on) and zero's elsewhere.
1359 The input matrix is not changed by this method in any way.
1361 · "$condition = $matrix->condition($inverse_matrix);"
1363 In fact this method is just a shortcut for
1365 abs($matrix) * abs($inverse_matrix)
1367 Both input matrices must be quadratic and have the same size, and
1369 the result is meaningful only if one of them is the inverse of the
1370 other (for instance, as returned by the method "invert_LR()").
1371 The number returned is a measure of the "condition" of the given
1373 matrix "$matrix", i.e., a measure of the numerical stability of the
1374 matrix.
1375 This number is always positive, and the smaller its value, the
1377 better the condition of the matrix (the better the stability of all
1378 subsequent computations carried out using this matrix).
1379 Numerical stability means for example that if
1381 abs( $vec_correct - $vec_with_error ) < $epsilon
1383 holds, there must be a "$delta" which doesn't depend on the vector
1385 "$vec_correct" (nor "$vec_with_error", by the way) so that
1386 abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta
1388 also holds.
1390 · "$determinant = $LR_matrix->det_LR();"
1392 Calculates the determinant of a matrix, whose LR decomposition
1394 matrix "$LR_matrix" must be given (which must be a (quadratic)
1395 matrix returned by the method "decompose_LR()").
1396 In fact the determinant is a by-product of the LR decomposition: It
1398 is (in principle, that is, except for the sign) simply the product
1399 of the elements on the main diagonal (elements (1,1), (2,2), (3,3)
1400 and so on) of the LR decomposition matrix.
1401 (The sign is taken care of "magically" by this module)
1403 · "$order = $LR_matrix->order_LR();"
1405 Calculates the order (called "Rang" in German) of a matrix, whose
1407 LR decomposition matrix "$LR_matrix" must be given (which must be a
1408 (quadratic) matrix returned by the method "decompose_LR()").
1409 This number is a measure of the number of linear independent row
1411 and column vectors (= number of linear independent equations in the
1412 case of a matrix representing an equation system) of the matrix
1413 that was initially fed into "decompose_LR()".
1414 If "n" is the number of rows and columns of the (quadratic!)
1416 matrix, then "n - order" is the dimension of the solution space of
1417 the associated equation system.
1418 · "$rank = $LR_matrix->rank_LR();"
1420 This is an alias for the "order_LR()" function. The "order" is
1422 usually called the "rank" in the United States.
1423 · "$scalar_product = $vector1->scalar_product($vector2);"
1425 Returns the scalar product of vector "$vector1" and vector
1427 "$vector2".
1428 Both vectors must be column vectors (i.e., a matrix having several
1430 rows but only one column).
1431 This is a (more efficient!) shortcut for
1433 $temp = ~$vector1 * $vector2;
1435 $scalar_product = $temp->element(1,1);
1436 or the sum "i=1..n" of the products "vector1[i] * vector2[i]".
1438 Provided none of the two input vectors is the null vector, then the
1440 two vectors are orthogonal, i.e., have an angle of 90 degrees
1441 between them, exactly when their scalar product is zero, and vice-
1442 versa.
1443 · "$vector_product = $vector1->vector_product($vector2);"
1445 Returns the vector product of vector "$vector1" and vector
1447 "$vector2".
1448 Both vectors must be column vectors (i.e., a matrix having several
1450 rows but only one column).
1451 Currently, the vector product is only defined for 3 dimensions
1453 (i.e., vectors with 3 rows); all other vectors trigger an error
1454 message.
1455 In 3 dimensions, the vector product of two vectors "x" and "y" is
1457 defined as
1458 | x[1] y[1] e[1] |
1460 determinant | x[2] y[2] e[2] |
1461 | x[3] y[3] e[3] |
1462 where the ""x[i]"" and ""y[i]"" are the components of the two
1464 vectors "x" and "y", respectively, and the ""e[i]"" are unity
1465 vectors (i.e., vectors with a length equal to one) with a one in
1466 row "i" and zero's elsewhere (this means that you have numbers and
1467 vectors as elements in this matrix!).
1468 This determinant evaluates to the rather simple formula
1470 z[1] = x[2] * y[3] - x[3] * y[2]
1472 z[2] = x[3] * y[1] - x[1] * y[3]
1473 z[3] = x[1] * y[2] - x[2] * y[1]
1474 A characteristic property of the vector product is that the
1476 resulting vector is orthogonal to both of the input vectors (if
1477 neither of both is the null vector, otherwise this is trivial),
1478 i.e., the scalar product of each of the input vectors with the
1479 resulting vector is always zero.
1480 · "$length = $vector->length();"
1482 This is actually a shortcut for
1484 $length = sqrt( $vector->scalar_product($vector) );
1486 and returns the length of a given column or row vector "$vector".
1488 Note that the "length" calculated by this method is in fact the
1490 "two"-norm (also know as the Euclidean norm) of a vector "$vector"!
1491 The general definition for norms of vectors is the following:
1493 sub vector_norm
1495 {
1496 croak "Usage: \$norm = \$vector->vector_norm(\$n);"
1497 if (@_ != 2);
1498 my($vector,$n) = @_;
1500 my($rows,$cols) = ($vector->[1],$vector->[2]);
1501 my($k,$comp,$sum);
1502 croak "Math::MatrixReal::vector_norm(): vector is not a column vector"
1504 unless ($cols == 1);
1505 croak "Math::MatrixReal::vector_norm(): norm index must be > 0"
1507 unless ($n > 0);
1508 croak "Math::MatrixReal::vector_norm(): norm index must be integer"
1510 unless ($n == int($n));
1511 $sum = 0;
1513 for ( $k = 0; $k < $rows; $k++ )
1514 {
1515 $comp = abs( $vector->[0][$k][0] );
1516 $sum += $comp ** $n;
1517 }
1518 return( $sum ** (1 / $n) );
1519 }
1520 Note that the case "n = 1" is the "one"-norm for matrices applied
1522 to a vector, the case "n = 2" is the euclidian norm or length of a
1523 vector, and if "n" goes to infinity, you have the "infinity"- or
1524 "maximum"-norm for matrices applied to a vector!
1525 · "$xn_vector = $matrix->""solve_GSM($x0_vector,$b_vector,$epsilon);"
1527 · "$xn_vector = $matrix->""solve_SSM($x0_vector,$b_vector,$epsilon);"
1529 · "$xn_vector =
1531 $matrix->""solve_RM($x0_vector,$b_vector,$weight,$epsilon);"
1532 In some cases it might not be practical or desirable to solve an
1534 equation system ""A * x = b"" using an analytical algorithm like
1535 the "decompose_LR()" and "solve_LR()" method pair.
1536 In fact in some cases, due to the numerical properties (the
1538 "condition") of the matrix "A", the numerical error of the obtained
1539 result can be greater than by using an approximative (iterative)
1540 algorithm like one of the three implemented here.
1541 All three methods, GSM ("Global Step Method" or
1543 "Gesamtschrittverfahren"), SSM ("Single Step Method" or
1544 "Einzelschrittverfahren") and RM ("Relaxation Method" or
1545 "Relaxationsverfahren"), are fix-point iterations, that is, can be
1546 described by an iteration function ""x(t+1) = Phi( x(t) )"" which
1547 has the property:
1548 Phi(x) = x <==> A * x = b
1550 We can define "Phi(x)" as follows:
1552 Phi(x) := ( En - A ) * x + b
1554 where "En" is a matrix of the same size as "A" ("n" rows and
1556 columns) with one's on its main diagonal and zero's elsewhere.
1557 This function has the required property.
1559 Proof:
1561 A * x = b
1563 <==> -( A * x ) = -b
1565 <==> -( A * x ) + x = -b + x
1567 <==> -( A * x ) + x + b = x
1569 <==> x - ( A * x ) + b = x
1571 <==> ( En - A ) * x + b = x
1573 This last step is true because
1575 x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
1577 is the same as
1579 ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
1581 qed
1583 Note that actually solving the equation system ""A * x = b"" means
1585 to calculate
1586 a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] = b[i]
1588 <==> a[i,i] x[i] =
1590 b[i]
1591 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
1592 + a[i,i] x[i]
1593 <==> x[i] =
1595 ( b[i]
1596 - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
1597 + a[i,i] x[i]
1598 ) / a[i,i]
1599 <==> x[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 There is one major restriction, though: a fix-point iteration is
1607 guaranteed to converge only if the first derivative of the
1608 iteration function has an absolute value less than one in an area
1609 around the point "x(*)" for which ""Phi( x(*) ) = x(*)"" is to be
1610 true, and if the start vector "x(0)" lies within that area!
1611 This is best verified graphically, which unfortunately is
1613 impossible to do in this textual documentation!
1614 See literature on Numerical Analysis for details!
1616 In our case, this restriction translates to the following three
1618 conditions:
1619 There must exist a norm so that the norm of the matrix of the
1621 iteration function, "( En - A )", has a value less than one, the
1622 matrix "A" may not have any zero value on its main diagonal and the
1623 initial vector "x(0)" must be "good enough", i.e., "close enough"
1624 to the solution "x(*)".
1625 (Remember school math: the first derivative of a straight line
1627 given by ""y = a * x + b"" is "a"!)
1628 The three methods expect a (quadratic!) matrix "$matrix" as their
1630 first argument, a start vector "$x0_vector", a vector "$b_vector"
1631 (which is the vector "b" in your equation system ""A * x = b""), in
1632 the case of the "Relaxation Method" ("RM"), a real number "$weight"
1633 best between zero and two, and finally an error limit (real number)
1634 "$epsilon".
1635 (Note that the weight "$weight" used by the "Relaxation Method"
1637 ("RM") is NOT checked to lie within any reasonable range!)
1638 The three methods first test the first two conditions of the three
1640 conditions listed above and return an empty list if these
1641 conditions are not fulfilled.
1642 Therefore, you should always test their return value using some
1644 code like:
1645 if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
1647 {
1648 # do something with the solution...
1649 }
1650 else
1651 {
1652 # do something with the fact that there is no solution...
1653 }
1654 Otherwise, they iterate until "abs( Phi(x) - x ) < epsilon".
1656 (Beware that theoretically, infinite loops might result if the
1658 starting vector is too far "off" the solution! In practice, this
1659 shouldn't be a problem. Anyway, you can always press <ctrl-C> if
1660 you think that the iteration takes too long!)
1661 The difference between the three methods is the following:
1663 In the "Global Step Method" ("GSM"), the new vector ""x(t+1)""
1665 (called "y" here) is calculated from the vector "x(t)" (called "x"
1666 here) according to the formula:
1667 y[i] =
1669 ( b[i]
1670 - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
1671 a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
1672 ) / a[i,i]
1673 In the "Single Step Method" ("SSM"), the components of the vector
1675 ""x(t+1)"" which have already been calculated are used to calculate
1676 the remaining components, i.e.
1677 y[i] =
1679 ( b[i]
1680 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
1681 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
1682 ) / a[i,i]
1683 In the "Relaxation method" ("RM"), the components of the vector
1685 ""x(t+1)"" are calculated by "mixing" old and new value (like cold
1686 and hot water), and the weight "$weight" determines the "aperture"
1687 of both the "hot water tap" as well as of the "cold water tap",
1688 according to the formula:
1689 y[i] =
1691 ( b[i]
1692 - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] + # note the "y[]"!
1693 a[i,i+1] x[i+1] + ... + a[i,n] x[n] ) # note the "x[]"!
1694 ) / a[i,i]
1695 y[i] = weight * y[i] + (1 - weight) * x[i]
1696 Note that the weight "$weight" should be greater than zero and less
1698 than two (!).
1699 The three methods are supposed to be of different efficiency.
1701 Experiment!
1702 Remember that in most cases, it is probably advantageous to first
1704 "normalize()" your equation system prior to solving it!
1705
1707 SYNOPSIS
1708 · Unary operators:
1709 ""-"", ""~"", ""abs"", "test", ""!"", '""'
1711 · Binary operators:
1713 "".""
1715 Binary (arithmetic) operators:
1717 ""+"", ""-"", ""*"", ""**"", ""+="", ""-="", ""*="", ""/="",""**=""
1719 · Binary (relational) operators:
1721 ""=="", ""!="", ""<"", ""<="", "">"", "">=""
1723 ""eq"", ""ne"", ""lt"", ""le"", ""gt"", ""ge""
1725 Note that the latter (""eq"", ""ne"", ... ) are just synonyms of the
1727 former (""=="", ""!="", ... ), defined for convenience only.
1728 DESCRIPTION
1730 '.' Concatenation
1731 Returns the two matrices concatenated side by side.
1733 Example: $c = $a . $b;
1735 For example, if
1737 $a=[ 1 2 ] $b=[ 5 6 ]
1739 [ 3 4 ] [ 7 8 ]
1740 then
1741 $c=[ 1 2 5 6 ]
1743 [ 3 4 7 8 ]
1744 Note that only matrices with the same number of rows may be
1746 concatenated.
1747 '-' Unary minus
1749 Returns the negative of the given matrix, i.e., the matrix with
1751 all elements multiplied with the factor "-1".
1752 Example:
1754 $matrix = -$matrix;
1756 '~' Transposition
1758 Returns the transposed of the given matrix.
1760 Examples:
1762 $temp = ~$vector * $vector;
1764 $length = sqrt( $temp->element(1,1) );
1765 if (~$matrix == $matrix) { # matrix is symmetric ... }
1767 abs Norm
1769 Returns the "one"-Norm of the given matrix.
1771 Example:
1773 $error = abs( $A * $x - $b );
1775 test Boolean test
1777 Tests wether there is at least one non-zero element in the matrix.
1779 Example:
1781 if ($xn_vector) { # result of iteration is not zero ... }
1783 '!' Negated boolean test
1785 Tests wether the matrix contains only zero's.
1787 Examples:
1789 if (! $b_vector) { # heterogenous equation system ... }
1791 else { # homogenous equation system ... }
1792 unless ($x_vector) { # not the null-vector! }
1794 '""""'
1796 "Stringify" operator
1797 Converts the given matrix into a string.
1799 Uses scientific representation to keep precision loss to a minimum
1801 in case you want to read this string back in again later with
1802 "new_from_string()".
1803 By default a 13-digit mantissa and a 20-character field for each
1805 element is used so that lines will wrap nicely on an 80-column
1806 screen.
1807 Examples:
1809 $matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
1811 [ 1 0 ]
1812 [ 0 -1 ]
1813 MATRIX
1814 print "$matrix";
1815 [ 1.000000000000E+00 0.000000000000E+00 ]
1817 [ 0.000000000000E+00 -1.000000000000E+00 ]
1818 $string = "$matrix";
1820 $test = Math::MatrixReal->new_from_string($string);
1821 if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }
1822 '+' Addition
1824 Returns the sum of the two given matrices.
1826 Examples:
1828 $matrix_S = $matrix_A + $matrix_B;
1830 $matrix_A += $matrix_B;
1832 '-' Subtraction
1834 Returns the difference of the two given matrices.
1836 Examples:
1838 $matrix_D = $matrix_A - $matrix_B;
1840 $matrix_A -= $matrix_B;
1842 Note that this is the same as:
1844 $matrix_S = $matrix_A + -$matrix_B;
1846 $matrix_A += -$matrix_B;
1848 (The latter are less efficient, though)
1850 '*' Multiplication
1852 Returns the matrix product of the two given matrices or the
1854 product of the given matrix and scalar factor.
1855 Examples:
1857 $matrix_P = $matrix_A * $matrix_B;
1859 $matrix_A *= $matrix_B;
1861 $vector_b = $matrix_A * $vector_x;
1863 $matrix_B = -1 * $matrix_A;
1865 $matrix_B = $matrix_A * -1;
1867 $matrix_A *= -1;
1869 '/' Division
1871 Currently a shortcut for doing $a * $b ** -1 is $a / $b, which
1873 works for square matrices. One can also use 1/$a .
1874 '**' Exponentiation
1876 Returns the matrix raised to an integer power. If 0 is passed, the
1878 identity matrix is returned. If a negative integer is passed, it
1879 computes the inverse (if it exists) and then raised the inverse to
1880 the absolute value of the integer. The matrix must be quadratic.
1881 Examples:
1883 $matrix2 = $matrix ** 2;
1885 $matrix **= 2;
1887 $inv2 = $matrix ** -2;
1889 $ident = $matrix ** 0;
1891 '==' Equality
1893 Tests two matrices for equality.
1895 Example:
1897 if ( $A * $x == $b ) { print "EUREKA!\n"; }
1899 Note that in most cases, due to numerical errors (due to the
1901 finite precision of computer arithmetics), it is a bad idea to
1902 compare two matrices or vectors this way.
1903 Better use the norm of the difference of the two matrices you want
1905 to compare and compare that norm with a small number, like this:
1906 if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }
1908 '!=' Inequality
1910 Tests two matrices for inequality.
1912 Example:
1914 while ($x0_vector != $xn_vector) { # proceed with iteration ... }
1916 (Stops when the iteration becomes stationary)
1918 Note that (just like with the '==' operator), it is usually a bad
1920 idea to compare matrices or vectors this way. Compare the norm of
1921 the difference of the two matrices with a small number instead.
1922 '<' Less than
1924 Examples:
1926 if ( $matrix1 < $matrix2 ) { # ... }
1928 if ( $vector < $epsilon ) { # ... }
1930 if ( 1E-12 < $vector ) { # ... }
1932 if ( $A * $x - $b < 1E-12 ) { # ... }
1934 These are just shortcuts for saying:
1936 if ( abs($matrix1) < abs($matrix2) ) { # ... }
1938 if ( abs($vector) < abs($epsilon) ) { # ... }
1940 if ( abs(1E-12) < abs($vector) ) { # ... }
1942 if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }
1944 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
1946 scalars.
1947 '<=' Less than or equal
1949 As with the '<' operator, this is just a shortcut for the same
1951 expression with "abs()" around all arguments.
1952 Example:
1954 if ( $A * $x - $b <= 1E-12 ) { # ... }
1956 which in fact is the same as:
1958 if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }
1960 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
1962 scalars.
1963 '>' Greater than
1965 As with the '<' and '<=' operator, this
1967 if ( $xn - $x0 > 1E-12 ) { # ... }
1969 is just a shortcut for:
1971 if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }
1973 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
1975 scalars.
1976 '>=' Greater than or equal
1978 As with the '<', '<=' and '>' operator, the following
1980 if ( $LR >= $A ) { # ... }
1982 is simply a shortcut for:
1984 if ( abs($LR) >= abs($A) ) { # ... }
1986 Uses the "one"-norm for matrices and Perl's built-in "abs()" for
1988 scalars.
1989
1991 Math::VectorReal, Math::PARI, Math::MatrixBool, Math::Vec, DFA::Kleene,
1992 Math::Kleene, Set::IntegerRange, Set::IntegerFast .
1993
1995 This man page documents Math::MatrixReal version 2.05.
1996 The latest released version can always be found at
1998 http://leto.net/code/Math-MatrixReal/ and there is also a subversion
1999 repository available at http://leto.net/svn/Math-MatrixReal/ .
2000
2002 Steffen Beyer <sb@engelschall.com>, Rodolphe Ortalo <ortalo@laas.fr>,
2003 Jonathan Leto <jonathan@leto.net>.
2004 Currently maintained by Jonathan Leto, send all bugs/patches to me.
2006
2008 Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for
2009 Algebra and Linear Algebra at the university (RWTH Aachen, Germany),
2010 and to Prof. Esser and his assistant, Mr. Jarausch, for their
2011 fascinating lectures in Numerical Analysis!
2012
2014 Copyright (c) 1996-2008 by Steffen Beyer, Rodolphe Ortalo, Jonathan
2015 Leto. All rights reserved.
2016
2018 This package is free software; you can redistribute it and/or modify it
2019 under the same terms as Perl itself.
2020perl v5.12.0 2010-05-03 Math::MatrixReal(3)