Math::MatrixReal(3pm)

📦
version 2.13 · 3 versions
, in debian bookworm, debian trixie, ubuntu jammy, ubuntu noble, ubuntu resolute, ubuntu stonking

💖 Math::MatrixReal 💖

okay so bestie, welcome to the Math::MatrixReal(3pm) page, brought to you by User Contributed Perl Documentation. this is the whole thing, start to finish, and we're doing it in one sitting. grab a snack. 🍓

🌸 NAME 🌸

Math::MatrixReal — Matrix of Reals

It implements the data type "matrix of real numbers" (and consequently also "vector of real numbers").

✨ SYNOPSIS ✨

my $a = Math::MatrixReal->new_random(5, 5);

my $b = $a->new_random(10, 30, { symmetric=>1, bounded_by=>[-1,1] });

my $c = $b * $a ** 3;

my $d = $b->new_from_rows( [ [ 5, 3 ,4], [3, 4, 5], [ 2, 4, 1 ] ] );

print $a;

my $row        = ($a * $b)->row(3);

my $col        = (5*$c)->col(2);

my $transpose  = ~$c;

my $transpose  = $c->transpose;

my $inverse    = $a->inverse;

my $inverse    = 1/$a;

my $inverse    = $a ** -1;

my $determinant= $a->det;
  • $matrix->display_precision($integer) — sets the default precision for when matrices get printed or stringified. $matrix->display_precision(0) will only show the integer part of all the entries of $matrix, and $matrix->display_precision() will return to the default scientific display notation. heads up tho: this method does not effect the precision of the calculations. it's a ✨display✨ thing only, bestie.

🧚 FUNCTIONS 🧚

💅 Constructor Methods And Such 💅

  • use Math::MatrixReal; — makes the methods and overloaded operators of this module available to your program. the classic. 🫶
  • $new_matrix = new Math::MatrixReal($rows,$columns); — the matrix object constructor method, the main character itself. a new matrix of size $rows by $columns gets created, with the value 0.0 for all elements. ngl, this method gets implicitly called by a whole bunch of the other methods in this module, so it's everywhere. 👑
  • $new_matrix = $some_matrix->new($rows,$columns); — another way of calling the matrix object constructor method. and no cap: matrix $some_matrix is not changed by this in any way. 💅
  • $new_matrix = $matrix->new_from_cols( [ $column_vector|$array_ref|$string, ... ] ) — creates a new matrix given a reference to an array of any of the following:
    • column vectors ( n by 1 Math::MatrixReal matrices )
    • references to arrays
    • strings properly formatted to create a column with Math::MatrixReal's new_from_string command

You may mix and match these as you wish. However, all must be of the same dimension--no padding happens automatically. Example:

my $matrix = Math::MatrixReal->new_from_cols( [ [1,2], [3,4] ] );
print $matrix;

will print

[  1.000000000000E+00  3.000000000000E+00 ]
[  2.000000000000E+00  4.000000000000E+00 ]
  • new_from_rows( [ $row_vector|$array_ref|$string, ... ] ) — creates a new matrix given a reference to an array of any of the following:
    • row vectors ( 1 by n Math::MatrixReal matrices )
    • references to arrays
    • strings properly formatted to create a row with Math::MatrixReal's new_from_string command

You may mix and match these as you wish. However, all must be of the same dimension--no padding happens automatically. Example:

my $matrix = Math::MatrixReal->new_from_rows( [ [1,2], [3,4] ] );
print $matrix;

will print

[  1.000000000000E+00  2.000000000000E+00 ]
[  3.000000000000E+00  4.000000000000E+00 ]
  • $new_matrix = Math::MatrixReal->new_random($rows, $cols, %options ); — this method allows you to create a random matrix with various properties controlled by the %options matrix, which is optional. the default values of the %options matrix are { integer => 0, symmetric => 0, tridiagonal => 0, diagonal => 0, bounded_by => [0,10] } .

Example:

$matrix = Math::MatrixReal->new_random(4, { diagonal => 1, integer => 1 }  );
print $matrix;

will print a 4x4 random diagonal matrix with integer entries between zero and ten, something like

[  5.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  2.000000000000E+00  0.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  0.000000000000E+00  1.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  8.000000000000E+00 ]
  • $new_matrix = Math::MatrixReal->new_diag( $array_ref ); — this method allows you to create a diagonal matrix by only specifying the diagonal elements. Example:
$matrix = Math::MatrixReal->new_diag( [ 1,2,3,4 ] );
print $matrix;

will print

[  1.000000000000E+00  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  2.000000000000E+00  0.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  0.000000000000E+00  3.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  0.000000000000E+00  0.000000000000E+00  4.000000000000E+00 ]
  • $new_matrix = Math::MatrixReal->new_tridiag( $lower, $diag, $upper ); — this method allows you to create a tridiagonal matrix by only specifying the lower diagonal, diagonal and upper diagonal, respectively. 🦋
$matrix = Math::MatrixReal->new_tridiag( [ 6, 4, 2 ], [1,2,3,4], [1, 8, 9] );
print $matrix;

will print

[  1.000000000000E+00  1.000000000000E+00  0.000000000000E+00  0.000000000000E+00 ]
[  6.000000000000E+00  2.000000000000E+00  8.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00  4.000000000000E+00  3.000000000000E+00  9.000000000000E+00 ]
[  0.000000000000E+00  0.000000000000E+00  2.000000000000E+00  4.000000000000E+00 ]
  • $new_matrix = Math::MatrixReal->new_from_string($string); — this method allows you to read in a matrix from a string (for instance, from the keyboard, from a file or from your code). iconic honestly.

The syntax is simple: each row must start with [ and end with ]\n (\n being the newline character and a space or tab) and contain one or more numbers, all separated from each other by spaces or tabs.

Additional spaces or tabs can be added at will, but no comments. 💀

Examples:

$string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
$matrix = Math::MatrixReal->new_from_string($string);
print "$matrix";

By the way, this prints

[  1.000000000000E+00  2.000000000000E+00  3.000000000000E+00 ]
[  2.000000000000E+00  2.000000000000E+00 -1.000000000000E+00 ]
[  1.000000000000E+00  1.000000000000E+00  1.000000000000E+00 ]

But you can also do this in a much more comfortable way using the shell-like "here-document" syntax:

$matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
[  1  0  0  0  0  0  1  ]
[  0  1  0  0  0  0  0  ]
[  0  0  1  0  0  0  0  ]
[  0  0  0  1  0  0  0  ]
[  0  0  0  0  1  0  0  ]
[  0  0  0  0  0  1  0  ]
[  1  0  0  0  0  0 -1  ]
MATRIX

You can even use variables in the matrix:

$c1 =   2  /  3;
$c2 =  -2  /  5;
$c3 =  26  /  9;

$matrix = Math::MatrixReal->new_from_string(<<"MATRIX");

    [   3    2    0   ]
    [   0    3    2   ]
    [  $c1  $c2  $c3  ]

MATRIX

(Remember that you may use spaces and tabs to format the matrix to your taste)

Note that this method uses exactly the same representation for a matrix as the "stringify" operator: this means that you can convert any matrix into a string with $string = "$matrix"; and read it back in later (for instance from a file!).

Note however that you may suffer a precision loss in this process because only 13 digits are supported in the mantissa when printed!! 😭

If the string you supply (or someone else supplies) does not obey the syntax mentioned above, an exception is raised, which can be caught by eval as follows:

print "Please enter your matrix (in one line): ";
$string = <STDIN>;
$string =~ s/\\n/\n/g;
eval { $matrix = Math::MatrixReal->new_from_string($string); };
if ($@)
{
    print "$@";
    # ...
    # (error handling)
}
else
{
    # continue...
}

or as follows:

eval { $matrix = Math::MatrixReal->new_from_string(<<"MATRIX"); };
[   3    2    0   ]
[   0    3    2   ]
[  $c1  $c2  $c3  ]
MATRIX
if ($@)
# ...

Actually, the method shown above for reading a matrix from the keyboard is a little awkward, since you have to enter a lot of "\n"'s for the newlines. 🙈

A better way is shown in this piece of code:

while (1)
{
    print "\nPlease enter your matrix ";
    print "(multiple lines, <ctrl-D> = done):\n";
    eval { $new_matrix =
        Math::MatrixReal->new_from_string(join('',<STDIN>)); };
    if ($@)
    {
        $@ =~ s/\s+at\b.*?$//;
        print "${@}Please try again.\n";
    }
    else { last; }
}

Possible error messages of the new_from_string() method are:

Math::MatrixReal::new_from_string(): syntax error in input string
Math::MatrixReal::new_from_string(): empty input string

If the input string has rows with varying numbers of columns, the following warning will be printed to STDERR:

Math::MatrixReal::new_from_string(): missing elements will be set to zero!

If everything is okay, the method returns an object reference to the (newly allocated) matrix containing the elements you specified. slay. ✨

  • $new_matrix = $some_matrix->shadow(); — returns an object reference to a NEW but EMPTY matrix (filled with zero's) of the SAME SIZE as matrix $some_matrix. and again, for the people in the back: matrix $some_matrix is not changed by this in any way. 💅
  • $matrix1->copy($matrix2); — copies the contents of matrix $matrix2 to an ALREADY EXISTING matrix $matrix1 (which must have the same size as matrix $matrix2!). same size, bestie, it's not optional. Matrix $matrix2 is not changed by this in any way.
  • $twin_matrix = $some_matrix->clone(); — returns an object reference to a NEW matrix of the SAME SIZE as matrix $some_matrix. The contents of matrix $some_matrix have ALREADY BEEN COPIED to the new matrix $twin_matrix. This is the method that the operator = is overloaded to when you type $a = $b, when $a and $b are matrices. Matrix $some_matrix is not changed by this in any way. twins!! 👯
  • $matrix = Math::MatrixReal->reshape($rows, $cols, $array_ref);

Return a matrix with the dimensions you asked for ($rows x $cols) where the elements get pulled straight out of that array reference $array_ref. And bestie, listen — the elements are read in column-major order (literally the same way Fortran arrays are stored, iconic behavior).

$matrix = Math::MatrixReal->reshape(4, 3, [1..12]);

Which gives you this absolute baddie of a matrix:

[ 1    5    9 ]
[ 2    6   10 ]
[ 3    7   11 ]
[ 4    8   12 ]

✨ Matrix Row, Column and Element operations ✨

  • $value = $matrix->element($row,$column); — hands you back the value of one specific element of $matrix, the one chilling in row $row and column $column.
NOTE: Unlike Perl, matrices are indexed with base-one indexes. So the very first element of the matrix lives in the first line, first column: perl $elem = $matrix->element(1, 1); # first element of the matrix.
  • $matrix->assign($row,$column,$value); — explicitly slaps the value $value onto a single element of $matrix, the one at row $row and column $column, which means whatever was sitting there before is gone, replaced, no cap.
  • $row_vector = $matrix->row($row); — a projection method 🪞 that returns an object reference to a NEW matrix (which is technically a (row) vector because it's only got one row) that already has row number $row of $matrix copied into it.

Matrix $matrix is not changed by this in any way.

  • $column_vector = $matrix->column($column); — another projection method 🪞, this time returning an object reference to a NEW matrix (a (column) vector, since it's got exactly one column) that already has column number $column of $matrix copied into it.

Matrix $matrix is not changed by this in any way.

  • @all_elements = $matrix->as_list; — grabs the contents of a Math::MatrixReal object and gives them to you as a Perl list. Obsessed.

Example:

my $matrix = Math::MatrixReal->new_from_rows([ [1, 2], [3, 4] ]);
my @list = $matrix->as_list; # 1, 2, 3, 4

This method is suitable for use with OpenGL. Say, for example, you need to rotate a model around the X-axis by 90 degrees clockwise. That could be achieved via:

use Math::Trig;
use OpenGL;
...;
my $axis = [1, 0, 0];
my $angle = 90;
...
my ($x, $y, $z) = @$axis;
my $f = $angle;
my $cos_f = cos(deg2rad($f));
my $sin_f = sin(deg2rad($f));
my $rotation = Math::MatrixReal->new_from_rows([
   [$cos_f+(1-$cos_f)*$x**2,    (1-$cos_f)*$x*$y-$sin_f*$z, (1-$cos_f)*$x*$z+$sin_f*$y, 0 ],
   [(1-$cos_f)*$y*$z+$sin_f*$z, $cos_f+(1-$cos_f)*$y**2 ,   (1-$cos_f)*$y*$z-$sin_f*$x, 0 ],
   [(1-$cos_f)*$z*$x-$sin_f*$y, (1-$cos_f)*$z*$y+$sin_f*$x, $cos_f+(1-$cos_f)*$z**2    ,0 ],
   [0,                          0,                          0,                          1 ],
]);
...;
my $model_initial = Math::MatrixReal->new_diag( [1, 1, 1, 1] ); # identity matrix
my $model = $model_initial * $rotation;
$model = ~$model; # OpenGL operates on transposed matrices
my $model_oga = OpenGL::Array->new_list(GL_FLOAT, $model->as_list);
$shader->SetMatrix(model => $model_oga); # instance of OpenGL::Shader

See OpenGL, OpenGL::Shader, OpenGL::Array, and rotation matrix over at https://en.wikipedia.org/wiki/Rotation_matrix for the lore 🦋

  • $new_matrix = $matrix->each( \&function ); — builds a brand new matrix by running a code reference over every single element of the given matrix. The function gets handed the element, the row index and the column index, in that order. Whatever the function returns (or the value of the last executed statement) becomes the value of the corresponding element in $new_matrix.

Example:

# add 1 to every element in the matrix
$matrix = $matrix->each ( sub { (shift) + 1 } );

Example:

my $cofactor = $matrix->each( sub { my(undef,$i,$j) = @_;
    ($i+$j) % 2 == 0 ? $matrix->minor($i,$j)->det()
    : -1*$matrix->minor($i,$j)->det();
    } );

This code needs some explaining, bestie. For each element of $matrix, it throws away the actual value and stashes the row and column indexes in $i and $j. Then it sets element [$i,$j] in $cofactor to the determinant of $matrix->minor($i,$j) if it's an "even" element, or -1*$matrix->minor($i,$j) if it's an "odd" element. Math ate and left no crumbs 💅

  • $new_matrix = $matrix->each_diag( \&function ); — same energy, except it builds a new matrix by evaluating a code reference on each diagonal element of the given matrix. The function gets the element, the row index and the column index, in that order. Whatever it returns (or the value of the last executed statement) is the value given to the corresponding element in $new_matrix.
  • $matrix->swap_col( $col1, $col2 ); — takes two one-based column numbers and swaps the values of every element in each column. So $matrix->swap_col(2,3) would replace column 2 in $matrix with column 3, and replace column 3 with column 2. Trade offer accepted 🤝
  • $matrix->swap_row( $row1, $row2 ); — takes two one-based row numbers and swaps the values of every element in each row. $matrix->swap_row(2,3) would replace row 2 in $matrix with row 3, and replace row 3 with row 2.
  • $matrix->assign_row( $row_number , $new_row_vector ); — takes a one-based row number and assigns row $row_number of $matrix with $new_row_vector, returning the resulting matrix. So $matrix->assign_row(5, $x) would replace row 5 in $matrix with the row vector $x.
  • $matrix->maximum(); and $matrix->minimum(); — these two work similarly, one for computing the maximum element or elements from a matrix, the other for the minimum element or elements from a matrix. They work in a similar way to the Octave/MatLab max/min functions.

When computing the maximum or minimum from a vector (vertical or horizontal), only one element is returned. When computing the maximum or minimum from a matrix, the maximum/minimum element for each column is returned in an array reference.

When called in list context, the function returns a pair, where the first element is the maximum/minimum element (or elements) and the second is the position of that value in the vector (first occurrence), or the row where it occurs, for matrices.

Consider the matrix and vector below for the following examples:

           [ 1 9 4 ]
      $A = [ 3 5 2 ]       $B = [ 8 7 9 5 3 ]
           [ 8 7 6 ]

When used in scalar context:

$max = $A->maximum();    # $max = [ 8, 9, 6 ]
$min = $B->minimum();    # $min = 3

When used in list context:

($min, $pos) = $A->minimum(); # $min = [ 1 5 2 ]
                              # $pos = [ 1 2 2 ]
($max, $pos) = $B->maximum(); # $max = 9
                              # $pos = 3

🔥 Matrix Operations 🔥

  • $det = $matrix->det(); — returns the determinant of the matrix, without going through the whole rigamarole of computing an LR decomposition. This method should be much faster than LR decomposition if the matrix is diagonal or triangular. Otherwise, it's just a wrapper for $matrix->decompose_LR->det_LR. If the determinant is zero, there is no inverse, and vice-versa. Only quadratic matrices have determinants.
  • $inverse = $matrix->inverse(); — returns the inverse of a matrix, without going through the rigamarole of computing an LR decomposition. If no inverse exists, undef is returned and an error is printed via carp(). This is nothing but a wrapper for $matrix->decompose_LR->invert_LR.
  • ($rows,$columns) = $matrix->dim(); — returns a list of two items: the number of rows and the number of columns that the given matrix $matrix contains.
  • $norm_one = $matrix->norm_one(); — returns the "one"-norm of the given matrix $matrix.

The "one"-norm is defined as follows:

For each column, the sum of the absolute values of the elements in the different rows of that column is calculated. Finally, the maximum of these sums is returned.

Note that the "one"-norm and the "maximum"-norm are mathematically equivalent, although for the same matrix they usually yield a different value.

Therefore, you should only compare values that have been calculated using the same norm!

Throughout this package, the "one"-norm is (arbitrarily) used for all comparisons, for the sake of uniformity and comparability, except for the iterative methods solve_GSM(), solve_SSM() and solve_RM() which use either norm depending on the matrix itself.

  • $norm_max = $matrix->norm_max(); — returns the "maximum"-norm of the given matrix $matrix.

The "maximum"-norm is defined as follows:

For each row, the sum of the absolute values of the elements in the different columns of that row is calculated. Finally, the maximum of these sums is returned.

Note that the "maximum"-norm and the "one"-norm are mathematically equivalent, although for the same matrix they usually yield a different value.

Therefore, you should only compare values that have been calculated using the same norm!

Throughout this package, the "one"-norm is (arbitrarily) used for all comparisons, for the sake of uniformity and comparability, except for the iterative methods solve_GSM(), solve_SSM() and solve_RM() which use either norm depending on the matrix itself.

  • $norm_sum = $matrix->norm_sum(); — this is a very simple norm which is defined as the sum of the absolute values of every element. Lowkey the easiest one here.
  • $p_norm = $matrix->norm_p($n);> — this function returns the "p-norm" of a vector. The argument $n must be a number greater than or equal to 1, or the string "Inf". The p-norm is defined as (sum(x_i^p))^(1/p). In words: it raises each element to the p-th power, adds them all up, and then takes the p-th root of that number. If the string "Inf" is passed, the "infinity-norm" is computed, which is really the limit of the p-norm as p goes to infinity. It's defined as the maximum element of the vector. Also, note that the familiar Euclidean distance between two vectors is just a special case of a p-norm, when p is equal to 2. Slay 🌈

Example vibes incoming, bestie 💖:

$a = Math::MatrixReal->new_from_cols([[1,2,3]]);
$p1   = $a->norm_p(1);
    $p2   = $a->norm_p(2);
    $p3   = $a->norm_p(3);
$pinf = $a->norm_p("Inf");

print "(1,2,3,Inf) norm:\n$p1\n$p2\n$p3\n$pinf\n";

$i1 = $a->new_from_rows([[1,0]]);
$i2 = $a->new_from_rows([[0,1]]);

# this should be sqrt(2) since it is the same as the
# hypotenuse of a 1 by 1 right triangle

$dist  = ($i1-$i2)->norm_p(2);
print "Distance is $dist, which should be " . sqrt(2) . "\n";

And the output, bestie — it literally ate 💅:

(1,2,3,Inf) norm:
6
3.74165738677394139
3.30192724889462668
3

Distance is 1.41421356237309505, which should be 1.41421356237309505
  • $frob_norm = "$matrix->norm_frobenius();" — okay so this norm is giving p-norm-with-p=2 energy, except it does its thing on a matrix, not a vector, okay? Every single element of the matrix gets squared, all of that gets added up, and then you take a square root. That's the whole vibe. 🌈
  • "$matrix->spectral_radius();" — this hands you the biggest absolute value out of all the eigenvalues. Fr fr though: right now it computes every eigenvalue and then sifts through them hunting for the largest in absolute value. Needless to say this is very inefficient 😭, and someday an algorithm that computes only the largest eigenvalue may be implemented. Not today though.
  • "$matrix1->transpose($matrix2);" — this calculates the transposed matrix of matrix $matrix2 and drops the result into matrix $matrix1. And heads up bestie: $matrix1 has to already exist AND have the same size as matrix $matrix2! 👑

You can also do this "in-place", like, input and output matrix can literally be the identical one. No cap.

Transposition is a symmetry operation — picture rotating the matrix 180 degrees along the axis of its main diagonal (the one going through elements (1,1), (2,2), (3,3) and so on).

Another way to clock it: rows and columns get swapped. In fact, the contents of element (i,j) get swapped with those of element (j,i).

Note that (especially for vectors) it makes a huge difference whether you've got a row vector, like this:

[ -1 0 1 ]

or a column vector, like this:

[ -1 ]
[  0 ]
[  1 ]

the one vector being the transposed of the other, periodt!

This is especially true for the matrix product of two vectors:

               [ -1 ]
[ -1 0 1 ] * [  0 ]  =  [ 2 ] ,  whereas
               [  1 ]

                             *     [ -1  0  1 ]
[ -1 ]                                            [  1  0 -1 ]
[  0 ] * [ -1 0 1 ]  =  [ -1 ]   [  1  0 -1 ]  =  [  0  0  0 ]
[  1 ]                  [  0 ]   [  0  0  0 ]     [ -1  0  1 ]
                        [  1 ]   [ -1  0  1 ]

So be careful about what you really mean, okay?? 😭

Hint: throughout this module, whenever a vector is explicitly required for input, a COLUMN vector is expected! 🎀
  • "$trace = $matrix->trace();" — returns the trace of the matrix, which is defined as the sum of the diagonal elements. The matrix must be quadratic. That's it, that's the whole thing ✨.
  • "$minor = $matrix->minor($row,$col);" — returns the minor matrix corresponding to $row and $col. $matrix must be quadratic. If $matrix is n rows by n cols, the minor of $row and $col will be an (n-1) by (n-1) matrix. The minor is defined as crossing out the row and the col you specified and returning the remaining rows and columns as a matrix. This method is used by "cofactor()". 💜
  • "$cofactor = $matrix->cofactor();" — the cofactor matrix gets constructed like this:

For each element, cross out the row and column that it sits in. Now, take the determinant of the matrix that is left in the other rows and columns. Multiply the determinant by (-1)^(i+j), where i is the row index, and j is the column index. Replace the given element with this value.

The cofactor matrix can be used to find the inverse of the matrix. One formula for the inverse of a matrix is the cofactor matrix transposed divided by the original determinant of the matrix.

The following two inverses should be exactly the same:

my $inverse1 = $matrix->inverse;
my $inverse2 = ~($matrix->cofactor)->each( sub { (shift)/$matrix->det() } );
Caveat bestie 🚨: although the cofactor matrix is a simple algorithm to compute the inverse of a matrix, and can be used with pencil and paper for small matrices, it is comically slower than the native "inverse()" function. Here is a small benchmark:
# $matrix1 is 15x15
$det = $matrix1->det;
timethese( 10,
    {'inverse' => sub { $matrix1->inverse(); },
      'cofactor' => sub { (~$matrix1->cofactor)->each ( sub { (shift)/$det; } ) }
    } );

Benchmark: timing 10 iterations of LR, cofactor, inverse...
    inverse:  1 wallclock secs ( 0.56 usr +  0.00 sys =  0.56 CPU) @ 17.86/s (n=10)
cofactor: 36 wallclock secs (36.62 usr +  0.01 sys = 36.63 CPU) @  0.27/s (n=10)
  • "$adjoint = $matrix->adjoint();" — the adjoint is just the transpose of the cofactor matrix. This method is literally just an alias for " ~($matrix->cofactor)". Lowkey the easiest glow-up in this whole module 🌸
  • "$part_of_matrix = $matrix->submatrix(x1,y1,x2,Y2);" — submatrix lets you select only part of an existing matrix so you can produce a brand new one. This method takes four arguments to define a selection area:
    • firstly: Coordinate of top left corner to select (x1,y1)
    • secondly: Coordinate of bottom right corner to select (x2,y2)

Example:

my $matrix = Math::MatrixReal->new_from_string(<<'MATRIX');
[  0  0  0  0  0  0  0  ]
[  0  0  0  0  0  0  0  ]
[  0  0  0  0  0  0  0  ]
[  0  0  0  0  0  0  0  ]
[  0  0  0  0  1  0  1  ]
[  0  0  0  0  0  1  0  ]
[  0  0  0  0  1  0  1  ]
MATRIX

my $submatrix = $matrix->submatrix(5,5,7,7);
$submatrix->display_precision(0);
print $submatrix;

Output:

[  1  0  1  ]
[  0  1  0  ]
[  1  0  1  ]

✨ Arithmetic Operations ✨

  • "$matrix1->add($matrix2,$matrix3);" — calculates the sum of matrix "$matrix2" and matrix "$matrix3" and stores the result in matrix "$matrix1". And yeah bestie, $matrix1 must already exist and have the same size as matrix "$matrix2" and matrix "$matrix3"! 👑

This operation can also be carried out "in-place", i.e., the output and one (or both) of the input matrices may be identical. Iconic.

  • "$matrix1->subtract($matrix2,$matrix3);" — calculates the difference of matrix "$matrix2" minus matrix "$matrix3" and stores the result in matrix "$matrix1" (which, yes, must already exist and have the same size as matrix "$matrix2" and matrix "$matrix3"!). 🦋

This operation can also be carried out "in-place", i.e., the output and one (or both) of the input matrices may be identical.

Note that this operation is the same as "$matrix1->add($matrix2,-$matrix3);", although the latter is a little less efficient. Ngl, just use subtract.

  • "$matrix1->multiply_scalar($matrix2,$scalar);" — calculates the product of matrix "$matrix2" and the number "$scalar" (i.e., multiplies each element of matrix "$matrix2" with the factor "$scalar") and stores the result in matrix "$matrix1" (which must already exist and have the same size as matrix "$matrix2"!). 🍓

This operation can also be carried out "in-place", i.e., input and output matrix may be identical.

  • "$product_matrix = $matrix1->multiply($matrix2);" — calculates the product of matrix "$matrix1" and matrix "$matrix2" and returns an object reference to a new matrix "$product_matrix" in which the result of this operation has been stored. ✨

Note that the dimensions of the two matrices "$matrix1" and "$matrix2" (i.e., their numbers of rows and columns) must harmonize in the following way (example):

                          [ 2 2 ]
                          [ 2 2 ]
                          [ 2 2 ]

              [ 1 1 1 ]   [ * * ]
              [ 1 1 1 ]   [ * * ]
              [ 1 1 1 ]   [ * * ]
              [ 1 1 1 ]   [ * * ]

I.e., the number of columns of matrix "$matrix1" has to be the same as the number of rows of matrix "$matrix2". Vibe check: they gotta match, no exceptions 💅

The number of rows and columns of the resulting matrix "$product_matrix" is determined by the number of rows of matrix "$matrix1" and the number of columns of matrix "$matrix2", respectively.

  • "$matrix1->negate($matrix2);" — calculates the negative of matrix "$matrix2" (i.e., multiplies all elements with "-1") and stores the result in matrix "$matrix1" (which must already exist and have the same size as matrix "$matrix2"!). 💀

This operation can also be carried out "in-place", i.e., input and output matrix may be identical.

  • "$matrix_to_power = $matrix1->exponent($integer);" — raises the matrix to the $integer power. Obviously, $integer must be an integer, bestie, don't try anything cute 🙈. If it is zero, the identity matrix is returned. If a negative integer is given, the inverse will be computed (if it exists) and then raised the the absolute value of $integer. The matrix must be quadratic.

🌈 Boolean Matrix Operations 🌈

  • "$matrix->is_quadratic();" — returns a boolean value indicating if the given matrix is quadratic (also know as "square" or "n by n"). A matrix is quadratic if it has the same number of rows as it does columns. Symmetry in its purest form ✨
  • "$matrix->is_square();" — this is an alias for "is_quadratic()". Same energy, different name 💅
  • "$matrix->is_symmetric();" — returns a boolean value indicating if the given matrix is symmetric. By definition, a matrix is symmetric if and only if (M[i,j]=M[j,i]). This is equivalent to "($matrix == ~$matrix)" but without memory allocation. Only quadratic matrices can be symmetric. 👑

Notes: A symmetric matrix always has real eigenvalues/eigenvectors. A matrix plus its transpose is always symmetric.

  • "$matrix->is_skew_symmetric();" — returns a boolean value indicating if the given matrix is skew symmetric. By definition, a matrix is symmetric if and only if (M[i,j]=-M[j,i]). This is equivalent to "($matrix == -(~$matrix))" but without memory allocation. Only quadratic matrices can be skew symmetric. 🦋
  • "$matrix->is_diagonal();"

Returns a boolean value telling you if the matrix is diagonal — like, every single nonzero element is parked on the main diagonal, nowhere else. Heads up though: only quadratic matrices can be diagonal. ✨

  • $matrix->is_tridiagonal(); — boolean check for tridiagonal energy: all the nonzero elements live on the main diagonal or on the diagonals directly above and below it. Same rule applies, bestie — only quadratic matrices can be tridiagonal. 🌸
  • $matrix->is_upper_triangular(); — boolean saying if the matrix is upper triangular, i.e. every nonzero element that is not on the main diagonal sits above it. Only quadratic matrices can be upper triangular. Also note: diagonal matrices are both upper and lower triangular. 💖
  • $matrix->is_lower_triangular(); — same vibe, flipped: lower triangular means every nonzero element that is not on the main diagonal sits below it. Only quadratic matrices can be lower triangular. And yes, again — diagonal matrices are both upper and lower triangular. 🦋
  • $matrix->is_orthogonal(); — boolean for whether the matrix is orthogonal. An orthogonal matrix has the property that the transpose equals the inverse of the matrix. Instead of computing each one and comparing them, this method multiplies the matrix by its transpose, and returns true if that turns out to be the identity matrix, false otherwise. Only quadratic matrices can be orthogonal. 👑
  • $matrix->is_binary(); — boolean telling you if the matrix is binary, aka it contains only zeroes or ones. 🔥
  • $matrix->is_gramian(); — boolean for whether the given matrix is Gramian. A matrix $A is Gramian if and only if there exists a square matrix $B such that "$A = ~$B*$B". That's equivalent to checking whether $A is symmetric and has all nonnegative eigenvalues — which is literally what Math::MatrixReal uses to check for this property. 🧚
  • $matrix->is_LR(); — boolean telling you if the matrix is an LR decomposition matrix. 💜
  • $matrix->is_positive(); — boolean for whether the matrix contains only positive entries. Careful though: a zero entry is not positive, and one little zero will make is_positive() return false. 💀
  • $matrix->is_negative(); — boolean for whether the matrix contains only negative entries. Same warning, same drama: a zero entry is not negative and it'll make is_negative() return false. 💀
  • $matrix->is_periodic($k); — boolean for whether the matrix is periodic with period $k. This is true if "$matrix ** ($k+1) == $matrix". And when "$k == 1", this reduces down to the is_idempotent() function. 🎀
  • $matrix->is_idempotent(); — boolean for whether the matrix is idempotent, which is defined as the square of the matrix being equal to the original matrix, i.e. "$matrix ** 2 == $matrix". ✨
  • $matrix->is_row_vector(); — boolean for whether the matrix is a row vector. A row vector is a matrix which is 1xn. Note that the 1x1 matrix is both a row and column vector. 🌸
  • $matrix->is_col_vector(); — boolean for whether the matrix is a col vector. A col vector is a matrix which is nx1. And same little plot twist: the 1x1 matrix is both a row and column vector. 🌈

✨ Eigensystems ✨

  • ($l, $V) = $matrix->sym_diagonalize(); — this one diagonalizes the quadratic symmetric matrix M stored in $matrix. On output, l is a column vector holding all the eigenvalues of M, and V is an orthogonal matrix whose columns are the corresponding normalized eigenvectors. The whole point of an eigenvalue l and an eigenvector x is of course that:
M * x = l * x

The method uses a Householder reduction to tridiagonal form, followed by a QL algorithm with implicit shifts on that tridiagonal. (The tridiagonal matrix is kept internally in a compact form in this routine, to save memory.) Fun fact: this routine actually wraps the householder() and tri_diagonalize() methods described below, for when you don't want their intermediate results. The overall algorithmic complexity of this technique is O(N^3). According to several books, the coefficient hidden by the 'O' is one of the best possible for general (symmetric) matrixes. 👑

  • ($T, $Q) = $matrix->householder(); — runs the Householder algorithm, which reduces the n by n real symmetric matrix M contained in $matrix to tridiagonal form. On output, T is a symmetric tridiagonal matrix (only diagonal and off-diagonal elements are non-zero) and Q is an orthogonal matrix performing the transformation between M and T:
$M == $Q * $T * ~$Q

🦋

  • ($l, $V) = $T->tri_diagonalize([$Q]); — diagonalizes the symmetric tridiagonal matrix T. On output, $l and $V are just like the output values described for sym_diagonalize().

The optional argument $Q corresponds to an orthogonal transformation matrix Q that should be used additionally during V (eigenvectors) computation. Supply it if the eigenvectors you want correspond to a more general symmetric matrix M that was previously reduced by the householder() method, not a mere tridiagonal. If T really is a tridiagonal matrix, Q can be omitted (it'll be internally created in fact as an identity matrix). The method uses a QL algorithm (with implicit shifts). 💅

  • $l = $matrix->sym_eigenvalues(); — computes the eigenvalues of the quadratic symmetric matrix M stored in $matrix. On output, l is a column vector containing all the eigenvalues of M. Eigenvectors are not computed (on the contrary of sym_diagonalize()) and this method is more efficient, even though it uses a similar algorithm with two phases. But understand: the algorithmic complexity of this technique is still also O(N^3). The coefficient hidden by the 'O' is better by a factor of... well, see your benchmark, it's wiser. 😌

This routine wraps the householder_tridiagonal() and tri_eigenvalues() methods described below, for when the intermediate tridiagonal matrix is not needed.

  • $T = $matrix->householder_tridiagonal(); — performs the Householder algorithm, which reduces the n by n real symmetric matrix M contained in $matrix to tridiagonal form. On output, T is the obtained symmetric tridiagonal matrix (only diagonal and off-diagonal elements are non-zero). The operation is similar to the householder() method, but potentially a little more efficient since the transformation matrix is not computed. 🔥
  • $l = $T->tri_eigenvalues(); — computes the eigenvalues of the symmetric tridiagonal matrix T. On output, $l is a vector containing the eigenvalues (similar to sym_eigenvalues()). This method is much more efficient than tri_diagonalize() when eigenvectors aren't needed. 🎀

✨ Miscellaneous ✨

  • $matrix->zero(); — assigns a zero to every element of the matrix $matrix, i.e., it erases all values previously stored there, thereby effectively transforming the matrix into a "zero"-matrix or "null"-matrix — the neutral element of the addition operation in a Ring. 💀

(For instance, the (quadratic) matrices with n rows and columns and matrix addition and multiplication form a Ring. Most prominent characteristic of a Ring is that multiplication is not commutative, i.e., in general, "matrix1 * matrix2" is not the same as "matrix2 * matrix1"!)

  • $matrix->one(); — assigns one's to the elements on the main diagonal (elements (1,1), (2,2), (3,3) and so on) of matrix $matrix, and zero's to all the others, thereby erasing all values previously stored there and transforming the matrix into a "one"-matrix — the neutral element of the multiplication operation in a Ring. 👑

(If the matrix is quadratic (which this method doesn't require, though), then multiplying this matrix with itself yields this same matrix again, and multiplying it with some other matrix leaves that other matrix unchanged!)

  • $latex_string = $matrix->as_latex( align=> "c", format => "%s", name => "" ); — returns the matrix as a LaTeX string. It takes a hash as an argument, which is used to control the style of the output. The hash element align may be "c", "l" or "r", corresponding to center, left and right, respectively. The format element is a format string that is given to sprintf to control the style of number format, such as floating point or scientific notation. The name element can be used so that a LaTeX string of $name = is prepended to the string. 💖

Example:

my $a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
print $a->as_latex( ( format => "%.2f", align => "l",name => "A" ) );

Output:

$A = $ $
\left( \begin{array}{ll}
1.23&1.00 \\
5.68&2.00 \\
9.10&3.00
\end{array} \right)
$
  • $yacas_string = $matrix->as_yacas( format => "%s", name => "", semi => 0 ); — returns the matrix as a string that can be read by Yacas. It takes a hash as an argument which controls the style of the output. The format element is a format string that is given to sprintf to control the style of number format, such as floating point or scientific notation. The name element can be used so that $name = is prepended to the string. The semi element can be set to 1 so that a semicolon is appended (so Matlab does not print out the matrix.) ✨

Example:

$a = Math::MatrixReal->new_from_cols([[ 1.234, 5.678, 9.1011],[1,2,3]] );
print $a->as_yacas( ( format => "%.2f", align => "l",name => "A" ) );

Output:

A := {{1.23,1.00},{5.68,2.00},{9.10,3.00}}
  • $matlab_string = $matrix->as_matlab( format => "%s", name => "", semi => 0 );

This one returns your matrix as a string that Matlab can actually read — like, straight-up ingest it. 💖 It takes a hash as an argument to control the styling of the output:

  • the format element is a format string that gets handed to sprintf to control the number format style, so floating point or scientific notation, your call
  • the name element lets you prepend $name = onto the string
  • the semi element can be set to 1 so a semicolon gets appended (which means Matlab won't print the matrix out)

Example:

my $a = Math::MatrixReal->new_from_rows([[ 1.234, 5.678, 9.1011],[1,2,3]] );
print $a->as_matlab( ( format => "%.3f", name => "A",semi => 1 ) );

Output:

A = [ 1.234 5.678 9.101;
 1.000 2.000 3.000];
  • $scilab_string = $matrix->as_scilab( format => "%s", name => "", semi => 0 ); — bestie, this is literally just an alias for as_matlab(), because Scilab and Matlab share the exact same matrix format. Same energy, different name. ✨
  • $minimum = Math::MatrixReal::min($number1,$number2); / $minimum = Math::MatrixReal::min($matrix); / $minimum = $matrix->min; — give it two numbers and you get the minimum of "number1" and "number2". Give it one argument, or call it as an object method, and it hands you the value of the smallest element in the matrix. 🌸
  • $maximum = Math::MatrixReal::max($number1,$number2); / $maximum = Math::MatrixReal::max($number1,$number2); / $maximum = Math::MatrixReal::max($matrix); / $maximum = $matrix->max; — same deal, opposite direction: two numbers in, maximum of "number1" and "number2" out. One argument, or object method, and you get the value of the largest element of the matrix. 👑
  • $minimal_cost_matrix = $cost_matrix->kleene(); — this one copies the matrix $cost_matrix (which has to be quadratic, no exceptions!) into a brand new matrix of the same size — i.e. it "clones" the input matrix — and then applies Kleene's algorithm to it. 🧚

Peep Math::Kleene(3) if you want the full tea on how this algorithm works!

The method returns an object reference to the new matrix.

Matrix $cost_matrix is not changed by this method in any way. Zero drama, promise. 💅

  • ($norm_matrix,$norm_vector) = $matrix->normalize($vector); — this method exists to level up the numerical stability when you're solving linear equation systems. 🔥

So imagine you've got a matrix A and a vector b, and you're hunting for a vector x such that A * x = b — that's the vector x that solves the equation system represented by matrix A and vector b.

Run this method on the pair (A,b) and you get back a pair (A',b') where every row has been divided by (the absolute value of) the greatest coefficient showing up in that row. Which means that coefficient becomes equal to 1 (or -1) in the new pair (A',b') — and all the others shrink down to smaller than one and greater than minus one.

And here's the thing: this operation does not change the equation system itself, because the same division happens on either side of the equation sign! It's giving balance. ⚖️

The method needs a quadratic (!) matrix $matrix and a vector $vector for input (the vector has to be a column vector with the same number of rows as the input matrix), and it returns a list of two items — object references to a new matrix and a new vector, in that order.

The output matrix and vector are clones of the input matrix and vector with the operation above already applied.

The input matrix and vector are not changed by this in any way. 💜

Example of how this method can affect the result of the methods that solve equation systems (explained right below, immediately following this method):

Consider the following little program:

#!perl -w

use Math::MatrixReal qw(new_from_string);

$A = Math::MatrixReal->new_from_string(<<"MATRIX");
[  1   2   3  ]
[  5   7  11  ]
[ 23  19  13  ]
MATRIX

$b = Math::MatrixReal->new_from_string(<<"MATRIX");
[   0   ]
[   1   ]
[  29   ]
MATRIX

$LR = $A->decompose_LR();
if (($dim,$x,$B) = $LR->solve_LR($b))
{
    $test = $A * $x;
    print "x = \n$x";
    print "A * x = \n$test";
}

($A_,$b_) = $A->normalize($b);

$LR = $A_->decompose_LR();
if (($dim,$x,$B) = $LR->solve_LR($b_))
{
    $test = $A * $x;
    print "x = \n$x";
    print "A * x = \n$test";
}

This will print:

x =
[  1.000000000000E+00 ]
[  1.000000000000E+00 ]
[ -1.000000000000E+00 ]
A * x =
[  4.440892098501E-16 ]
[  1.000000000000E+00 ]
[  2.900000000000E+01 ]
x =
[  1.000000000000E+00 ]
[  1.000000000000E+00 ]
[ -1.000000000000E+00 ]
A * x =
[  0.000000000000E+00 ]
[  1.000000000000E+00 ]
[  2.900000000000E+01 ]

See how in the second example (the one where normalize() was used) the result is "better", i.e. more accurate?! That's the glow-up. 😭✨

  • $LR_matrix = $matrix->decompose_LR(); — this method is the one you need to actually solve linear equation systems.

Say you've got a matrix A and a vector b and you wanna find a vector x so that A * x = b — i.e. the vector x that solves the equation system represented by matrix A and vector b.

Or maybe you've got a matrix A and a whole bunch of different vectors b1..bk and you need vectors x1..xk so that A * xi = bi, for i=1..k.

Using Gaussian transformations (multiplying a row or column by a factor, swapping two rows or two columns, and adding a multiple of one row or column to another), it's possible to decompose any matrix A into two triangular matrices, called L and R (for "Left" and "Right").

L has one's on the main diagonal (the elements (1,1), (2,2), (3,3) and so so), non-zero values to the left and below the main diagonal, and all zero's in the upper right half of the matrix.

R has non-zero values on the main diagonal as well as to the right and above of the main diagonal, and all zero's in the lower left half of the matrix, like this:

          [ 1 0 0 0 0 ]      [ x x x x x ]
          [ x 1 0 0 0 ]      [ 0 x x x x ]
      L = [ x x 1 0 0 ]  R = [ 0 0 x x x ]
          [ x x x 1 0 ]      [ 0 0 0 x x ]
          [ x x x x 1 ]      [ 0 0 0 0 x ]

Note that "L * R" is equivalent to matrix A in the sense that L * R * x = b <==> A * x = b for all vectors x — leaving out of account permutations of the rows and columns (those get taken care of "magically" by this module!) and numerical errors.

Trick:

Because we know that L has one's on its main diagonal, we can store both matrices together in the same array without losing any information! I.e.,

                 [ R R R R R ]
                 [ L R R R R ]
            LR = [ L L R R R ]
                 [ L L L R R ]
                 [ L L L L R ]
Beware though, bestie — LR and "L * R" are not the same thing!!! 💀

Note also that, for that same reason, you cannot apply the method normalize() to an LR decomposition matrix. Trying to do so will yield meaningless rubbish!

(You need to apply normalize() to each pair (Ai,bi) BEFORE decomposing the matrix Ai'!)

Now, what does all this actually do for us in solving linear equation systems?

It helps us because a triangular matrix is the next best thing that can happen to us besides a diagonal matrix (a matrix that has non-zero values only on its main diagonal — in which case the solution is trivial, you simply divide "b[i]" by "A[i,i]" to get "x[i]"!).

To find the solution to our problem "A * x = b", we divide the problem into parts: instead of solving A * x = b directly, we first decompose A into L and R, and then solve "L * y = b" and finally "R * x = y" (motto: divide and rule!).

From the illustration above it's clear that solving "L * y = b" and "R * x = y" is straightforward: we immediately know that y[1] = b[1]. We then deduce swiftly that

y[2] = b[2] - L[2,1] * y[1]

(and we know "y[1]" by now!), that

y[3] = b[3] - L[3,1] * y[1] - L[3,2] * y[2]

and so on.

Having effortlessly calculated the vector y, we now proceed to calculate the vector x in a similar fashion: we see immediately that x[n] = y[n] / R[n,n]. It follows that

x[n-1] = ( y[n-1] - R[n-1,n] * x[n] ) / R[n-1,n-1]

and

x[n-2] = ( y[n-2] - R[n-2,n-1] * x[n-1] - R[n-2,n] * x[n] )
         / R[n-2,n-2]

and so on.

You can see that — especially when you've got many vectors b1..bk that you're searching solutions to A * xi = bi for — this scheme is much more efficient than a straightforward, "brute force" approach. It's giving efficiency icon. 🦋

This method requires a quadratic matrix as its input matrix.

If you don't have that many equations, fill up with zero's (i.e., do nothing to fill the superfluous rows if it's a "fresh" matrix — a matrix that's been created with new() or shadow()).

The method returns an object reference to a new matrix containing the matrices L and R.

The input matrix is not changed by this method in any way.

Note that you can copy() or clone() the result of this method without losing its "magical" properties (for instance, the ones concerning the hidden permutations of its rows and columns).

However — as soon as you apply any method that alters the contents of the matrix, its "magical" properties get stripped off, and the matrix immediately reverts to an "ordinary" matrix (with the values it just happens to contain at that moment, whether those are meaningful as an ordinary matrix or not!). 🙈

  • ($dimension,$x_vector,$base_matrix) = $LR_matrix->solve_LR($b_vector); — use this method to actually solve an equation system.

Matrix $LR_matrix must be a (quadratic) matrix returned by the method decompose_LR() — the LR decomposition matrix of the matrix A of your equation system A * x = b.

The input vector $b_vector is the vector b in your equation system A * x = b, which must be a column vector and have the same number of rows as the input matrix $LR_matrix.

The method hands you back a list of three things if a solution exists — or an empty list if it doesn't (!). 😭

So like, always call it like this, bestie, no exceptions:
  if ( ($dim,$x_vec,$base) = $LR->solve_LR($b_vec) )
  {
      # do something with the solution...
  }
  else
  {
      # do something with the fact that there is no solution...
  }

The three things you get back are:

  • the dimension $dimension of the solution space (it's zero if there's only one solution, one if the solution is a straight line, two if the solution is a plane, and so on forever, the way that... yeah),
  • the solution vector $x_vector (which is the vector x in your equation system A * x = b), and
  • a matrix $base_matrix that's a base of the solution space (a set of vectors that hold up the solution space like the spokes of an umbrella ☂️).

Only the first $dimension columns of this base matrix actually contain entries, the remaining columns are all zero. Periodt. 💅

Now what is all this stuff with that "base" good for?

The output vector x is ALWAYS a solution of your equation system A * x = b.

But also any vector $vector

  $vector = $x_vector->clone();

  $machine_infinity = 1E+99; # or something like that

  for ( $i = 1; $i <= $dimension; $i++ )
  {
      $vector += rand($machine_infinity) * $base_matrix->column($i);
  }

is a solution to your problem A * x = b, i.e., if $A_matrix contains your matrix A, then

  print abs( $A_matrix * $vector - $b_vector ), "\n";

should print a number around 1E-16 or so! 🤯 iconic.

By the way, note that you can actually calculate those vectors $vector a little more efficient as follows:

  $rand_vector = $x_vector->shadow();

  $machine_infinity = 1E+99; # or something like that

  for ( $i = 1; $i <= $dimension; $i++ )
  {
      $rand_vector->assign($i,1, rand($machine_infinity) );
  }

  $vector = $x_vector + ( $base_matrix * $rand_vector );

Note that the input matrix and vector are not changed by this method in any way. 🫶

  • $inverse_matrix = $LR_matrix->invert_LR(); — use this method to calculate the inverse of a given matrix $LR_matrix, which must be a (quadratic) matrix returned by the method decompose_LR().

The method returns an object reference to a new matrix of the same size as the input matrix containing the inverse of the matrix that you initially fed into decompose_LR() IF THE INVERSE EXISTS, or an empty list otherwise.

So like, you should always use this method in the following way:

  if ( $inverse_matrix = $LR->invert_LR() )
  {
      # do something with the inverse matrix...
  }
  else
  {
      # do something with the fact that there is no inverse matrix...
  }

Note that by definition (disregarding numerical errors), the product of the initial matrix and its inverse (or vice-versa) is always a matrix containing one's on the main diagonal (elements (1,1), (2,2), (3,3) and so on) and zero's elsewhere. it's giving identity matrix, iykyk ✨

The input matrix is not changed by this method in any way.

  • $condition = $matrix->condition($inverse_matrix); — in fact this method is just a shortcut for
  abs($matrix) * abs($inverse_matrix)

Both input matrices must be quadratic and have the same size, and the result is meaningful only if one of them is the inverse of the other (for instance, as returned by the method invert_LR()).

The number returned is a measure of the "condition" of the given matrix $matrix, i.e., a measure of the numerical stability of the matrix.

This number is always positive, and the smaller its value, the better the condition of the matrix (the better the stability of all subsequent computations carried out using this matrix). lowkey the vibe check of your matrix 🌈

Numerical stability means for example that if

  abs( $vec_correct - $vec_with_error ) < $epsilon

holds, there must be a $delta which doesn't depend on the vector $vec_correct (nor $vec_with_error, by the way) so that

  abs( $matrix * $vec_correct - $matrix * $vec_with_error ) < $delta

also holds.

  • $determinant = $LR_matrix->det_LR(); — calculates the determinant of a matrix, whose LR decomposition matrix $LR_matrix must be given (which must be a (quadratic) matrix returned by the method decompose_LR()).

In fact the determinant is a by-product of the LR decomposition: It is (in principle, that is, except for the sign) simply the product of the elements on the main diagonal (elements (1,1), (2,2), (3,3) and so on) of the LR decomposition matrix.

(The sign is taken care of "magically" by this module 🧚✨)

  • $order = $LR_matrix->order_LR(); — calculates the order (called "Rang" in German) of a matrix, whose LR decomposition matrix $LR_matrix must be given (which must be a (quadratic) matrix returned by the method decompose_LR()).

This number is a measure of the number of linear independent row and column vectors (= number of linear independent equations in the case of a matrix representing an equation system) of the matrix that was initially fed into decompose_LR().

If n is the number of rows and columns of the (quadratic!) matrix, then n - order is the dimension of the solution space of the associated equation system. 👑

  • $rank = $LR_matrix->rank_LR(); — this is an alias for the order_LR() function. The "order" is usually called the "rank" in the United States. 🇺🇸
  • $scalar_product = $vector1->scalar_product($vector2); — returns the scalar product of vector $vector1 and vector $vector2.

Both vectors must be column vectors (i.e., a matrix having several rows but only one column).

This is a (more efficient!) shortcut for

  $temp           = ~$vector1 * $vector2;
  $scalar_product =  $temp->element(1,1);

or the sum i=1..n of the products vector1[i] * vector2[i].

Provided none of the two input vectors is the null vector, then the two vectors are orthogonal, i.e., have an angle of 90 degrees between them, exactly when their scalar product is zero, and vice-versa. 📐

  • $vector_product = $vector1->vector_product($vector2); — returns the vector product of vector $vector1 and vector $vector2.

Both vectors must be column vectors (i.e., a matrix having several rows but only one column).

Currently, the vector product is only defined for 3 dimensions (i.e., vectors with 3 rows); all other vectors trigger an error message. 💀

In 3 dimensions, the vector product of two vectors x and y is defined as

                |  x[1]  y[1]  e[1]  |
    determinant |  x[2]  y[2]  e[2]  |
                |  x[3]  y[3]  e[3]  |

where the x[i] and y[i] are the components of the two vectors x and y, respectively, and the e[i] are unity vectors (i.e., vectors with a length equal to one) with a one in row i and zero's elsewhere (this means that you have numbers and vectors as elements in this matrix!).

This determinant evaluates to the rather simple formula

  z[1] = x[2] * y[3] - x[3] * y[2]
  z[2] = x[3] * y[1] - x[1] * y[3]
  z[3] = x[1] * y[2] - x[2] * y[1]

A characteristic property of the vector product is that the resulting vector is orthogonal to both of the input vectors (if neither of both is the null vector, otherwise this is trivial), i.e., the scalar product of each of the input vectors with the resulting vector is always zero.

  • $length = $vector->length(); — this is actually a shortcut for
  $length = sqrt( $vector->scalar_product($vector) );

and returns the length of a given column or row vector $vector.

Note that the "length" calculated by this method is in fact the "two"-norm (also know as the Euclidean norm) of a vector $vector! 🌸

The general definition for norms of vectors is the following:

  sub vector_norm
  {
      croak "Usage: \$norm = \$vector->vector_norm(\$n);"
        if (@_ != 2);

      my($vector,$n) = @_;
      my($rows,$cols) = ($vector->[1],$vector->[2]);
      my($k,$comp,$sum);

      croak "Math::MatrixReal::vector_norm(): vector is not a column vector"
        unless ($cols == 1);

      croak "Math::MatrixReal::vector_norm(): norm index must be > 0"
        unless ($n > 0);

      croak "Math::MatrixReal::vector_norm(): norm index must be integer"
        unless ($n == int($n));

      $sum = 0;
      for ( $k = 0; $k < $rows; $k++ )
      {
          $comp = abs( $vector->[0][$k][0] );
          $sum += $comp ** $n;
      }
      return( $sum ** (1 / $n) );
  }

Note that the case n = 1 is the "one"-norm for matrices applied to a vector, the case n = 2 is the euclidian norm or length of a vector, and if n goes to infinity, you have the "infinity"- or "maximum"-norm for matrices applied to a vector! 🚀

  • $xn_vector = $matrix->solve_GSM($x0_vector,$b_vector,$epsilon);
  • $xn_vector = $matrix->solve_SSM($x0_vector,$b_vector,$epsilon);
  • $xn_vector = $matrix->solve_RM($x0_vector,$b_vector,$weight,$epsilon);

In some cases it might not be practical or desirable to solve an equation system A * x = b using an analytical algorithm like the decompose_LR() and solve_LR() method pair.

In fact in some cases, due to the numerical properties (the "condition") of the matrix A, the numerical error of the obtained result can be greater than by using an approximative (iterative) algorithm like one of the three implemented here. 🙈

All three methods, GSM ("Global Step Method" or "Gesamtschrittverfahren"), SSM ("Single Step Method" or "Einzelschrittverfahren") and RM ("Relaxation Method" or "Relaxationsverfahren"), are fix-point iterations, that is, can be described by an iteration function x(t+1) = Phi( x(t) ) which has the property:

  Phi(x)  =  x    <==>    A * x  =  b

We can define Phi(x) as follows:

  Phi(x)  :=  ( En - A ) * x  +  b

where En is a matrix of the same size as A (n rows and columns) with one's on its main diagonal and zero's elsewhere.

This function has the required property.

Proof:

           A * x        =   b
             <==>  -( A * x )      =  -b

             <==>  -( A * x ) + x  =  -b + x

             <==>  -( A * x ) + x + b  =  x

             <==>  x - ( A * x ) + b  =  x

             <==>  ( En - A ) * x + b  =  x

and this last glow-up is legit because this:

             x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]

is literally the same thing as

             ( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]

qed bestie 💅

okay so real talk — actually solving the equation system A * x = b means you're out here calculating:

                   a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n]  =  b[i]

             <==>  a[i,i] x[i]  =
                   b[i]
                   - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
                   + a[i,i] x[i]

             <==>  x[i]  =
                   ( b[i]
                       - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] )
                       + a[i,i] x[i]
                   ) / a[i,i]

             <==>  x[i]  =
                   ( b[i] -
                       ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
                         a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
                   ) / a[i,i]
⚠️ but there's one MAJOR catch, and it's giving plot twist: a fix-point iteration is only guaranteed to converge if the first derivative of the iteration function has an absolute value less than one in an area around the point x(*) where Phi( x(*) ) = x(*) is supposed to be true — and only if the start vector x(0) sits inside that area! 🔥

And the way you'd normally check that is graphically, which — sadly — is impossible to do in this textual documentation 😭

Go read some literature on Numerical Analysis for the deets, bestie 📚

So in our case, that restriction boils down to these three conditions, no cap:

  • There has to exist a norm such that the norm of the matrix of the iteration function, ( En - A ), is less than one.
  • The matrix A may not have any zero value on its main diagonal.
  • The initial vector x(0) must be "good enough", aka "close enough" to the solution x(*).

(flashback to school math rn: the first derivative of a straight line given by y = a * x + b is a! 💀)

All three methods want a (quadratic!) matrix $matrix as their first argument, a start vector $x0_vector, a vector $b_vector (that's the b in your A * x = b equation system), and — in the case of the "Relaxation Method" (RM) — a real number $weight best kept between zero and two, plus finally an error limit (real number) $epsilon.

heads up: the weight $weight used by the "Relaxation Method" (RM) is NOT checked to lie within any reasonable range! 🙈

The three methods first test the first two of those three conditions and return an empty list if they aren't fulfilled.

So like, always test their return value, okay? Some code like this:

  if ( $xn_vector = $A_matrix->solve_GSM($x0_vector,$b_vector,1E-12) )
  {
      # do something with the solution...
  }
  else
  {
      # do something with the fact that there is no solution...
  }

Otherwise they just iterate until abs( Phi(x) - x ) < epsilon.

⚠️ beware: theoretically, infinite loops can happen if the starting vector is too far "off" the solution! In practice it shouldn't be an issue, and anyway you can always smash <ctrl-C> if the iteration is taking way too long 😤

So what's the difference between the three? glad you asked ✨

💫 The Global Step Method (GSM)

In GSM, the new vector x(t+1) (called y here) is calculated from the vector x(t) (called x here) according to this formula:

  y[i] =
  ( b[i]
      - ( a[i,1] x[1] + ... + a[i,i-1] x[i-1] +
          a[i,i+1] x[i+1] + ... + a[i,n] x[n] )
  ) / a[i,i]

💫 The Single Step Method (SSM)

In SSM, the components of the vector x(t+1) that have already been calculated get used to calculate the remaining components, i.e.

  y[i] =
  ( b[i]
      - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] +  # note the "y[]"!
          a[i,i+1] x[i+1] + ... + a[i,n] x[n] )  # note the "x[]"!
  ) / a[i,i]

💫 The Relaxation Method (RM)

In RM, the components of the vector x(t+1) get calculated by "mixing" old and new values (think cold water and hot water 🌊🔥), and the weight $weight decides the "aperture" of both the hot-water tap and the cold-water tap, according to:

  y[i] =
  ( b[i]
      - ( a[i,1] y[1] + ... + a[i,i-1] y[i-1] +  # note the "y[]"!
          a[i,i+1] x[i+1] + ... + a[i,n] x[n] )  # note the "x[]"!
  ) / a[i,i]
  y[i] = weight * y[i] + (1 - weight) * x[i]
and remember: the weight $weight should be greater than zero and less than two (!) 🫶

The three methods are supposed to have different levels of efficiency. Experiment! 🧪

And remember that in most cases it's probably a major slay to first call normalize() on your equation system before solving it! 👑

✨ OVERLOADED OPERATORS ✨

Okay so Perl lets you put actual math symbols on your matrix objects, and that's the whole point of this section. It's giving operator overload, and honestly? Iconic. 💅

SYNOPSIS

  • Unary operators: ""-"", ""~"", ""abs"", "test", ""!"", '""'
  • Binary operators: "".""
  • Binary (arithmetic) operators: ""+"", ""-"", ""*"", ""**"", ""+="", ""-="", ""*="", ""/="",""**=""
  • Binary (relational) operators: ""=="", ""!="", ""<"", ""<="", "">"", "">=""
  • And the wordy versions: ""eq"", ""ne"", ""lt"", ""le"", ""gt"", ""ge""
Just so you know: the wordy ones (""eq"", ""ne"", ... ) are literally just synonyms of the symbol ones (""=="", ""!="", ... ), and they exist for convenience only. No cap. 💖

DESCRIPTION

. — Concatenation

Sticks the two matrices together side by side and gives you the combo back.

Example: $c = $a . $b;

For example, if

$a=[ 1 2 ]   $b=[ 5 6 ]
   [ 3 4 ]      [ 7 8 ]

then

$c=[ 1 2 5 6 ]
   [ 3 4 7 8 ]
Heads up: only matrices with the same number of rows can be concatenated. Periodt.

- — Unary minus

Gives you the negative of the matrix, aka every element gets multiplied by the factor "-1".

Example:

$matrix = -$matrix;

~ — Transposition

Returns the transpose of the matrix. Flip it and reverse it. 🔄

Examples:

$temp = ~$vector * $vector;
$length = sqrt( $temp->element(1,1) );

if (~$matrix == $matrix) { # matrix is symmetric ... }

abs — Norm

Returns the "one"-Norm of the matrix.

Example:

$error = abs( $A * $x - $b );

test — Boolean test

Checks if there's at least one non-zero element hanging out in the matrix.

Example:

if ($xn_vector) { # result of iteration is not zero ... }

! — Negated boolean test

Checks if the matrix is all zeroes and nothing but zeroes.

Examples:

if (! $b_vector) { # heterogenous equation system ... }
else             { # homogenous equation system ... }

unless ($x_vector) { # not the null-vector! }

'""""' — "Stringify" operator

Turns the matrix into a string.

It uses scientific representation so you lose as little precision as possible when you want to read that string back in later with new_from_string().

By default it uses a 13-digit mantissa and a 20-character field for each element, which makes the lines wrap all nice on an 80-column screen.

Examples:

$matrix = Math::MatrixReal->new_from_string(<<"MATRIX");
[ 1  0 ]
[ 0 -1 ]
MATRIX
print "$matrix";

[  1.000000000000E+00  0.000000000000E+00 ]
[  0.000000000000E+00 -1.000000000000E+00 ]

$string = "$matrix";
$test = Math::MatrixReal->new_from_string($string);
if ($test == $matrix) { print ":-)\n"; } else { print ":-(\n"; }

+ — Addition

Adds the two matrices together and returns the sum.

Examples:

$matrix_S = $matrix_A + $matrix_B;

$matrix_A += $matrix_B;

- — Subtraction

Subtracts the second matrix from the first and returns the difference.

Examples:

$matrix_D = $matrix_A - $matrix_B;

$matrix_A -= $matrix_B;

This is the same as:

$matrix_S = $matrix_A + -$matrix_B;

$matrix_A += -$matrix_B;
(The latter way is less efficient, though. Just saying.)

* — Multiplication

Returns the matrix product of the two matrices, or the product of the matrix and a scalar factor.

Examples:

$matrix_P = $matrix_A * $matrix_B;

$matrix_A *= $matrix_B;

$vector_b = $matrix_A * $vector_x;

$matrix_B = -1 * $matrix_A;

$matrix_B = $matrix_A * -1;

$matrix_A *= -1;

/ — Division

Right now this is a shortcut for doing $a * $b ** -1 is $a / $b, which works for square matrices. You can also use 1/$a .

** — Exponentiation

Raises the matrix to an integer power. If you pass 0, you get the identity matrix back. If you pass a negative integer, it computes the inverse (if it exists) and then raises that inverse to the absolute value of the integer. The matrix must be quadratic.

Examples:

$matrix2 = $matrix ** 2;

$matrix **= 2;

$inv2 = $matrix ** -2;

$ident = $matrix ** 0;

== — Equality

Tests two matrices for equality.

Example:

if ( $A * $x == $b ) { print "EUREKA!\n"; }
Ngl, in most cases, because of numerical errors (finite precision of computer arithmetic and all that), comparing two matrices or vectors this way is a bad idea.

Better move: use the norm of the difference of the two matrices you want to compare, and compare that norm with a small number, like this:

if ( abs( $A * $x - $b ) < 1E-12 ) { print "BINGO!\n"; }

!= — Inequality

Tests two matrices for inequality.

Example:

while ($x0_vector != $xn_vector) { # proceed with iteration ... }

(Stops when the iteration becomes stationary)

Same warning as with ==: it's usually a bad idea to compare matrices or vectors this way. Compare the norm of the difference of the two matrices with a small number instead.

< — Less than

Examples:

if ( $matrix1 < $matrix2 ) { # ... }

if ( $vector < $epsilon ) { # ... }

if ( 1E-12 < $vector ) { # ... }

if ( $A * $x - $b < 1E-12 ) { # ... }

These are just shortcuts for saying:

if ( abs($matrix1) < abs($matrix2) ) { # ... }

if ( abs($vector) < abs($epsilon) ) { # ... }

if ( abs(1E-12) < abs($vector) ) { # ... }

if ( abs( $A * $x - $b ) < abs(1E-12) ) { # ... }

It uses the "one"-Norm for matrices and Perl's built-in abs() for scalars.

<= — Less than or equal

Just like the < operator, this is a shortcut for the same expression but with abs() wrapped around all the arguments.

Example:

if ( $A * $x - $b <= 1E-12 ) { # ... }

which is actually the same as:

if ( abs( $A * $x - $b ) <= abs(1E-12) ) { # ... }

It uses the "one"-Norm for matrices and Perl's built-in abs() for scalars.

> — Greater than

Just like < and <=, this

if ( $xn - $x0 > 1E-12 ) { # ... }

is a shortcut for:

if ( abs( $xn - $x0 ) > abs(1E-12) ) { # ... }

It uses the "one"-Norm for matrices and Perl's built-in abs() for scalars.

>= — Greater than or equal

Just like <, <=, and >, the following

if ( $LR >= $A ) { # ... }

is simply a shortcut for:

if ( abs($LR) >= abs($A) ) { # ... }

It uses the "one"-Norm for matrices and Perl's built-in abs() for scalars.

👀 SEE ALSO 👀

Math::VectorReal, Math::PARI, Math::MatrixBool, Math::Vec, DFA::Kleene, Math::Kleene, Set::IntegerRange, Set::IntegerFast .

📦 VERSION 📦

This man page documents Math::MatrixReal version 2.13

The latest code can be found at https://github.com/leto/math--matrixreal .

👑 AUTHORS 👑

Steffen Beyer <sb@engelschall.com>, Rodolphe Ortalo <ortalo@laas.fr>, Jonathan "Duke" Leto <jonathan@leto.net>.

Currently maintained by Jonathan "Duke" Leto, send all bugs/patches to Github Issues: https://github.com/leto/math--matrixreal/issues

💖 CREDITS 💖

Many thanks to Prof. Pahlings for stoking the fire of my enthusiasm for Algebra and Linear Algebra at the university (RWTH Aachen, Germany), and to Prof. Esser and his assistant, Mr. Jarausch, for their fascinating lectures in Numerical Analysis!

Copyright (c) 1996-2016 by various authors including the original developer Steffen Beyer, Rodolphe Ortalo, the current maintainer Jonathan "Duke" Leto and all the wonderful people in the AUTHORS file. All rights reserved.

📜 LICENSE AGREEMENT 📜

This package is free software; you can redistribute it and/or modify it under the same terms as Perl itself. Fuck yeah.

And that's the whole tea, bestie! This page comes from perl v5.30.3, dated 2020-06-13, and it's the Math::MatrixReal(3pm) manpage. Now go overload those operators like the main character you are. 💅✨🦋

📦 where this page comes from