NAME

Math::Matrix - multiply and invert matrices

SYNOPSIS

use Math::Matrix;

# Generate a random 3-by-3 matrix.
srand(time);
my $A = Math::Matrix -> new([rand, rand, rand],
                            [rand, rand, rand],
                            [rand, rand, rand]);
$A -> print("A\n");

# Append a fourth column to $A.
my $x = Math::Matrix -> new([rand, rand, rand]);
my $E = $A -> concat($x -> transpose);
$E -> print("Equation system\n");

# Compute the solution.
my $s = $E -> solve;
$s -> print("Solutions s\n");

# Verify that the solution equals $x.
$A -> multiply($s) -> print("A*s\n");

DESCRIPTION

This module implements various constructors and methods for creating and manipulating matrices.

All methods return new objects, so, for example, $X->add($Y) does not modify $X.

$X -> add($Y);         # $X not modified; output is lost
$X = $X -> add($Y);    # this works

Some operators are overloaded (see "OVERLOADING") and allow the operand to be modified directly.

$X = $X + $Y;          # this works
$X += $Y;              # so does this

METHODS

Constructors

Identify matrices

Identify elements

This section contains methods for identifying and locating elements within an array. See also ["Scalar comparison"](#scalar-comparison).

Basic properties

Manipulate matrices

These methods are for combining matrices, splitting matrices, extracing parts of a matrix, inserting new parts into a matrix, deleting parts of a matrix etc. There are also methods for shuffling elements around (relocating elements) inside a matrix.

These methods are not concerned with the values of the elements.

Mathematical functions

Addition

Subtraction

Negation

Multiplication

Powers

Inversion

Factorisation

Miscellaneous matrix functions

Elementwise mathematical functions

These method work on each element of a matrix.

Columnwise or rowwise mathematical functions

These method work along each column or row of a matrix. Some of these methods return a matrix with the same size as the invocand matrix. Other methods collapse the dimension, so that, e.g., if the method is applied to the first dimension a p-by-q matrix becomes a 1-by-q matrix, and if applied to the second dimension, it becomes a p-by-1 matrix. Others, like ["diff()"](#diff), reduces the length along the dimension by one, so a p-by-q matrix becomes a (p-1)-by-q or a p-by-(q-1) matrix.

Comparison

Matrix comparison

Methods matrix comparison. These methods return a scalar value.

Scalar comparison

Each of these methods performs scalar (element by element) comparison and returns a matrix of ones and zeros. Scalar expansion is performed if necessary.

Vector functions

Conversion

Matrix utilities

Apply a subroutine to each element

Forward elimination

These methods take a matrix as input, performs forward elimination, and returns a matrix where all elements below the main diagonal are zero. In list context, four additional arguments are returned: an array with the row permutations, an array with the column permutations, an integer with the number of row swaps and an integer with the number of column swaps performed during elimination.

The permutation vectors can be converted to permutation matrices with ["to_permmat()"](#to_permmat).

Back-substitution

Miscellaneous methods

OVERLOADING

The following operators are overloaded.

IMPROVING THE SOLUTION OF LINEAR SYSTEMS

The methods that do an explicit or implicit matrix left division accept some additional parameters. If these parameters are specified, the matrix left division is done repeatedly in an iterative way, which often gives a better solution.

Background

The linear system of equations

$A * $x = $y

can be solved for $x with

$x = $y -> mldiv($A);

Ideally $A * $x should equal $y, but due to numerical errors, this is not always the case. The following illustrates how to improve the solution $x computed above:

$r = $A -> mmuladd($x, -$y);    # compute the residual $A*$x-$y
$d = $r -> mldiv($A);           # compute the delta for $x
$x -= $d;                       # improve the solution $x

This procedure is repeated, and at each step, the absolute error

||$A*$x - $y|| = ||$r||

and the relative error

||$A*$x - $y|| / ||$y|| = ||$r|| / ||$y||

are computed and compared to the tolerances. Once one of the stopping criteria is satisfied, the algorithm terminates.

Stopping criteria

The algorithm stops when at least one of the errors are within the specified tolerances or the maximum number of iterations is reached. If the maximum number of iterations is reached, but noen of the errors are within the tolerances, a warning is displayed and the best solution so far is returned.

Parameters

Example

If

$A = [[  8, -8, -5,  6, -1,  3 ],
      [ -7, -1,  5, -9,  5,  6 ],
      [ -7,  8,  9, -2, -4,  3 ],
      [  3, -4,  5,  5,  3,  3 ],
      [  9,  8, -3, -4,  1,  6 ],
      [ -8,  9, -1,  3,  5,  2 ]];

$y = [[  80, -13 ],
      [  -2, 104 ],
      [ -57, -27 ],
      [  47, -28 ],
      [   5,  77 ],
      [  91, 133 ]];

the result of $x = $y -> mldiv($A);, using double precision arithmetic, is the approximate solution

$x = [[ -2.999999999999998, -5.000000000000000 ],
      [ -1.000000000000000,  3.000000000000001 ],
      [ -5.999999999999997, -8.999999999999996 ],
      [  8.000000000000000, -2.000000000000003 ],
      [  6.000000000000003,  9.000000000000002 ],
      [  7.999999999999997,  8.999999999999995 ]];

The residual $res = $A -> mmuladd($x, -$y); is

$res = [[  1.24344978758018e-14,  1.77635683940025e-15 ],
        [  8.88178419700125e-15, -5.32907051820075e-15 ],
        [ -1.24344978758018e-14,  1.77635683940025e-15 ],
        [ -7.10542735760100e-15, -4.08562073062058e-14 ],
        [ -1.77635683940025e-14, -3.81916720471054e-14 ],
        [  1.24344978758018e-14,  8.43769498715119e-15 ]];

and the delta $dx = $res -> mldiv($A); is

$dx = [[   -8.592098303124e-16, -2.86724066474914e-15 ],
       [ -7.92220125658508e-16, -2.99693950082398e-15 ],
       [ -2.22533360993874e-16,  3.03465504177947e-16 ],
       [  6.47376093198353e-17, -1.12378127899388e-15 ],
       [  6.35204502123966e-16,  2.40938179521241e-15 ],
       [  1.55166908001001e-15,  2.08339859425849e-15 ]];

giving the improved, and in this case exact, solution $x -= $dx;,

$x = [[ -3, -5 ],
      [ -1,  3 ],
      [ -6, -9 ],
      [  8, -2 ],
      [  6,  9 ],
      [  8,  9 ]];

SUBCLASSING

The methods should work fine with any kind of numerical objects, provided that the assignment operator = is overloaded, so that Perl knows how to create a copy.

You can check the behaviour of the assignment operator by assigning a value to a new variable, modify the new variable, and check whether this also modifies the original value. Here is an example:

$x = Some::Class -> new(0);           # create object $x
$y = $x;                              # create new variable $y
$y++;                                 # modify $y
print "it's a clone\n" if $x != $y;   # is $x modified?

The subclass might need to implement some methods of its own. For instance, if each element is a complex number, a transpose() method needs to be implemented to take the complex conjugate of each value. An as_string() method might also be useful for displaying the matrix in a format more suitable for the subclass.

Here is an example showing Math::Matrix::Complex, a fully-working subclass of Math::Matrix, where each element is a Math::Complex object.

use strict;
use warnings;

package Math::Matrix::Complex;

use Math::Matrix;
use Scalar::Util 'blessed';
use Math::Complex 1.57;     # "=" didn't clone before 1.57

our @ISA = ('Math::Matrix');

# We need a new() method to make sure every element is an object.

sub new {
    my $self = shift;
    my $x = $self -> SUPER::new(@_);

    my $sub = sub {
        defined(blessed($_[0])) && $_[0] -> isa('Math::Complex')
          ? $_[0]
          : Math::Complex -> new($_[0]);
    };

    return $x -> sapply($sub);
}

# We need a transpose() method, since the transpose of a matrix
# with complex numbers also takes the conjugate of all elements.

sub transpose {
    my $x = shift;
    my $y = $x -> SUPER::transpose(@_);

    return $y -> sapply(sub { ~$_[0] });
}

# We need an as_string() method, since our parent's methods
# doesn't format complex numbers correctly.

sub as_string {
    my $self = shift;
    my $out = "";
    for my $row (@$self) {
        for my $elm (@$row) {
            $out = $out . sprintf "%10s ", $elm;
        }
        $out = $out . sprintf "\n";
    }
    $out;
}

1;

BUGS

Please report any bugs or feature requests via https://github.com/pjacklam/p5-Math-Matrix/issues.

Old bug reports and feature requests can be found at http://rt.cpan.org/NoAuth/Bugs.html?Dist=Math-Matrix.

SUPPORT

You can find documentation for this module with the perldoc command.

perldoc Math::Matrix

You can also look for information at:

LICENSE AND COPYRIGHT

Copyright (c) 2020-2021, Peter John Acklam

Copyright (C) 2013, John M. Gamble jgamble@ripco.com, all rights reserved.

Copyright (C) 2009, oshalla https://rt.cpan.org/Public/Bug/Display.html?id=42919

Copyright (C) 2002, Bill Denney gte273i@prism.gatech.edu, all rights reserved.

Copyright (C) 2001, Brian J. Watson bjbrew@power.net, all rights reserved.

Copyright (C) 2001, Ulrich Pfeifer pfeifer@wait.de, all rights reserved. Copyright (C) 1995, Universität Dortmund, all rights reserved.

Copyright (C) 2001, Matthew Brett matthew.brett@mrc-cbu.cam.ac.uk

This program is free software; you may redistribute it and/or modify it under the same terms as Perl itself.

AUTHORS

Peter John Acklam pjacklam@gmail.com (2020-2021)

Ulrich Pfeifer pfeifer@ls6.informatik.uni-dortmund.de (1995-2013)

Brian J. Watson bjbrew@power.net

Matthew Brett matthew.brett@mrc-cbu.cam.ac.uk