Tuesday, January 28, 2014

fun with PERL Math::Matrix

Wanting to do some matrix math, I installed the Math::Matrix module from CPAN. There's some nice instructions on installing CPAN modules on about.com.

A common tool people use for matrix manipulation is MATLAB. I'm deeply embarrassed I have virtually no skills in MATLAB. However, this module is pretty cool.

What makes it so nice it is uses overloading of the basic mathematical operators. So that makes it easy to add, subtract, and multiply matrices, as well as a string conversion when it's accessed in a scalar context. Here's an example:

use Math::Matrix;
use strict;
my $a = new Math::Matrix ([1, 0, -1], [-1, 1, 0], [0, -1, 1]);
warn("a =\n$a\n");
my $at = $a->transpose;
warn("at =\n$at\n");
warn("at * a =\n", $at * $a, "\n");

The result:

a =
   1.00000    0.00000   -1.00000 
  -1.00000    1.00000    0.00000 
   0.00000   -1.00000    1.00000 

at =
   1.00000   -1.00000    0.00000 
   0.00000    1.00000   -1.00000 
  -1.00000    0.00000    1.00000 

at * a =
   2.00000   -1.00000   -1.00000 
  -1.00000    2.00000   -1.00000 
  -1.00000   -1.00000    2.00000 

I can also solve equations. Here's a trivial set of 3 linear equations where the solution vector is listed as a fourth column to a 3×3 matrix of coefficients:

print new Math::Matrix ([1, 0, 0, 2], [0, 1, 0, 3], [0, 0, 1, 4]) -> solve;

The result:

   2.00000 
   3.00000 
   4.00000 
Or I can solve for multiple solution vectors:
print new Math::Matrix ([1, 0, 0, 2, rand], [0, 1, 0, 3, rand], [0, 0, 1, 4, rand]) -> solve;

This yields:


   2.00000    0.76763 
   3.00000    0.38500 
   4.00000    0.92288 

You can also concatenate the solution vector to the coefficient matrix:

my $matrix = Math::Matrix->diagonal(1, 1, 1);
warn "matrix =\n$matrix\n";
my $vector = new Math::Matrix([1],[2],[3]);
warn "vector =\n$vector\n";
warn "solution =\n", $matrix->concat($vector)->solve(), "\n";

This generates:

matrix =
   1.00000    0.00000    0.00000 
   0.00000    1.00000    0.00000 
   0.00000    0.00000    1.00000 

vector =
   1.00000 
   2.00000 
   3.00000 

solution =
   1.00000 
   2.00000 
   3.00000 

Here I'm using "warn" instead of "print". It's my general practice in scripts to use "warn" for user messages, and "print" for data intended to piping into another command.

If I want to set particular components of the matrix, the following approach seems to work:

my @m;
for my $r ( 0 .. 4 ) {
  for my $c ( 0 .. 3 ) {
    $m[$r]->[$c] = 0;
  }
}
my $m = new Math::Matrix @m;

$m->[2]->[1] = 42;
print "m =\n$m\n";

I constructed an array of references to arrays, passed that to the class new method, and it returned a matrix object. I set elements of the matrix using the same notation I would us to set elements of a reference to that original array of array references.

The output:

m =
   0.00000    0.00000    0.00000    0.00000 
   0.00000    0.00000    0.00000    0.00000 
   0.00000   42.00000    0.00000    0.00000 
   0.00000    0.00000    0.00000    0.00000 
   0.00000    0.00000    0.00000    0.00000 

The first index was the row index (2) and the second index was the column references (1).

The big one for solving matrix problems is probably inversion. That's done with invert. For example:

my $m = Math::Matrix->diagonal(1, 2, 3, 4);
my $i = $m->invert;
print "m =\n$m\nm inverse =\n$i";

This inverts, as an example, a diagonal matrix.

So this is all good stuff for my present project.

No comments: