Math-MatrixReal
view release on metacpan or search on metacpan
lib/Math/MatrixReal.pm view on Meta::CPAN
for (my $k = 0; $k < $i; $k++)
{
$g += $Q->[0][$i][$k] * $Q->[0][$k][$j];
}
for (my $k = 0; $k < $i; $k++)
{
$Q->[0][$k][$j] -= $g * $Q->[0][$k][$i];
}
}
$d->[$i] = $Q->[0][$i][$i];
# Reset row and column of Q for next iteration
$Q->[0][$i][$i] = 1.0;
for (my $j = 0; $j < $i; $j++)
{
$Q->[0][$i][$j] = $Q->[0][$j][$i] = 0.0;
}
}
return ($d, $e);
}
# Computes sqrt(a*a + b*b), but more carefully...
sub _pythag ($$)
{
my ($a, $b) = @_;
my $aa = abs($a);
my $ab = abs($b);
if ($aa > $ab)
{
# NB: Not needed!: return 0.0 if ($aa == 0.0);
my $t = $ab / $aa;
return ($aa * sqrt(1.0 + $t*$t));
} else {
return 0.0 if ($ab == 0.0);
my $t = $aa / $ab;
return ($ab * sqrt(1.0 + $t*$t));
}
}
# QL algorithm with implicit shifts to determine the eigenvalues
# of a tridiagonal matrix. Internal routine.
sub _tridiagonal_QLimplicit
{
my ($EV, $d, $e) = @_;
my ($rows, $cols) = ($EV->[1], $EV->[2]);
$e->[$rows-1] = 0.0;
# Start real computation
for (my $l = 0; $l < $rows; $l++)
{
my $iter = 0;
my $m;
OUTER:
do {
for ($m = $l; $m < ($rows - 1); $m++)
{
my $dd = abs($d->[$m]) + abs($d->[$m+1]);
last if ((abs($e->[$m]) + $dd) == $dd);
}
if ($m != $l)
{
## why only allow 30 iterations?
croak("Too many iterations!") if ($iter++ >= 30);
my $g = ($d->[$l+1] - $d->[$l])
/ (2.0 * $e->[$l]);
my $r = _pythag($g, 1.0);
$g = $d->[$m] - $d->[$l]
+ $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
my ($p,$s,$c) = (0.0, 1.0,1.0);
for (my $i = ($m-1); $i >= $l; $i--)
{
my $ii = $i + 1;
my $f = $s * $e->[$i];
my $t = _pythag($f, $g);
$e->[$ii] = $t;
if ($t == 0.0)
{
$d->[$ii] -= $p;
$e->[$m] = 0.0;
next OUTER;
}
my $b = $c * $e->[$i];
$s = $f / $t;
$c = $g / $t;
$g = $d->[$ii] - $p;
my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
$p = $s * $t2;
$d->[$ii] = $g + $p;
$g = $c * $t2 - $b;
for (my $k = 0; $k < $rows; $k++)
{
my $t1 = $EV->[0][$k][$ii];
my $t2 = $EV->[0][$k][$i];
$EV->[0][$k][$ii] = $s * $t2 + $c * $t1;
$EV->[0][$k][$i] = $c * $t2 - $s * $t1;
}
}
$d->[$l] -= $p;
$e->[$l] = $g;
$e->[$m] = 0.0;
}
} while ($m != $l);
}
return;
}
# Core householder reduction routine (when eigenvector
# are NOT wanted).
sub _householder_values ($)
{
my ($Q) = @_; # NB: Q is destroyed on output...
my ($rows, $cols) = ($Q->[1], $Q->[2]);
# Creates tridiagonal
# Set up tridiagonal needed elements
my $d = []; # N Diagonal elements 0...n-1
my $e = []; # N-1 Off-Diagonal elements 0...n-2
my @p = ();
for (my $i = ($rows - 1); $i > 1; $i--)
{
my $scale = 0.0;
for (my $k = 0; $k < $i; $k++)
lib/Math/MatrixReal.pm view on Meta::CPAN
{
my $g = 0.0;
for (my $k = 0; $k <= $j; $k++)
{
$g += $Q->[0][$j][$k] * $Q->[0][$i][$k];
}
for (my $k = $j+1; $k < $i; $k++)
{
$g += $Q->[0][$k][$j] * $Q->[0][$i][$k];
}
# Form elements of P
$p[$j] = $g / $h;
$f += $p[$j] * $Q->[0][$i][$j];
}
my $hh = $f / ($h + $h);
for (my $j = 0; $j < $i; $j++)
{
my $t = $Q->[0][$i][$j];
my $g = $p[$j] - $hh * $t;
$p[$j] = $g;
for (my $k = 0; $k <= $j; $k++)
{
$Q->[0][$j][$k] -= $t * $p[$k]
+ $g * $Q->[0][$i][$k];
}
}
}
}
# Updates for i==1
$e->[0] = $Q->[0][1][0];
# Updates diagonal elements
for (my $i = 0; $i < $rows; $i++)
{
$d->[$i] = $Q->[0][$i][$i];
}
return ($d, $e);
}
# QL algorithm with implicit shifts to determine the
# eigenvalues ONLY. This is O(N^2) only...
sub _tridiagonal_QLimplicit_values
{
my ($M, $d, $e) = @_; # NB: M is not touched...
my ($rows, $cols) = ($M->[1], $M->[2]);
$e->[$rows-1] = 0.0;
# Start real computation
for (my $l = 0; $l < $rows; $l++)
{
my $iter = 0;
my $m;
OUTER:
do {
for ($m = $l; $m < ($rows - 1); $m++)
{
my $dd = abs($d->[$m]) + abs($d->[$m+1]);
last if ((abs($e->[$m]) + $dd) == $dd);
}
if ($m != $l)
{
croak("Too many iterations!") if ($iter++ >= 30);
my $g = ($d->[$l+1] - $d->[$l])
/ (2.0 * $e->[$l]);
my $r = _pythag($g, 1.0);
$g = $d->[$m] - $d->[$l]
+ $e->[$l] / ($g + (($g >= 0.0) ? abs($r) : -abs($r)));
my ($p,$s,$c) = (0.0, 1.0,1.0);
for (my $i = ($m-1); $i >= $l; $i--)
{
my $ii = $i + 1;
my $f = $s * $e->[$i];
my $t = _pythag($f, $g);
$e->[$ii] = $t;
if ($t == 0.0)
{
$d->[$ii] -= $p;
$e->[$m] = 0.0;
next OUTER;
}
my $b = $c * $e->[$i];
$s = $f / $t;
$c = $g / $t;
$g = $d->[$ii] - $p;
my $t2 = ($d->[$i] - $g) * $s + 2.0 * $c * $b;
$p = $s * $t2;
$d->[$ii] = $g + $p;
$g = $c * $t2 - $b;
}
$d->[$l] -= $p;
$e->[$l] = $g;
$e->[$m] = 0.0;
}
} while ($m != $l);
}
return;
}
# Householder reduction of a real, symmetric matrix A.
# Returns a tridiagonal matrix T and the orthogonal matrix
# Q effecting the transformation between A and T.
sub householder ($)
{
my ($A) = @_;
my ($rows, $cols) = ($A->[1], $A->[2]);
croak "Matrix is not quadratic"
unless ($rows = $cols);
croak "Matrix is not symmetric"
unless ($A->is_symmetric());
# Copy given matrix TODO: study if we should do in-place modification
my $Q = $A->clone();
# Do the computation of tridiagonal elements and of
# transformation matrix
my ($diag, $offdiag) = $Q->_householder_vectors();
# Creates the tridiagonal matrix
my $T = $A->shadow();
for (my $i = 0; $i < $rows; $i++)
{ # Set diagonal
lib/Math/MatrixReal.pm view on Meta::CPAN
[ 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!
Hint: throughout this module, whenever a vector is explicitly required
for input, a B<COLUMN> vector is expected!
=item *
C<$trace = $matrix-E<gt>trace();>
This returns the trace of the matrix, which is defined as
the sum of the diagonal elements. The matrix must be
quadratic.
=item *
C<$minor = $matrix-E<gt>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 specified and returning
the remaining rows and columns as a matrix. This method is used by C<cofactor()>.
=item *
C<$cofactor = $matrix-E<gt>cofactor();>
The cofactor matrix is constructed as follows:
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: Although the cofactor matrix is 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 C<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)
=item *
C<$adjoint = $matrix-E<gt>adjoint();>
The adjoint is just the transpose of the cofactor matrix. This method is
just an alias for C< ~($matrix-E<gt>cofactor)>.
=back
=item *
C<$part_of_matrix = $matrix-E<gt>submatrix(x1,y1,x2,Y2);>
Submatrix permit to select only part of existing matrix in order to produce a new one.
This method take four arguments to define a selection area:
=over 6
=item - firstly: Coordinate of top left corner to select (x1,y1)
=item - secondly: Coordinate of bottom right corner to select (x2,y2)
=back
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 ]
=back
=head2 Arithmetic Operations
=over 4
=item *
C<$matrix1-E<gt>add($matrix2,$matrix3);>
Calculates the sum of matrix "C<$matrix2>" and matrix "C<$matrix3>"
lib/Math/MatrixReal.pm view on Meta::CPAN
Note that the "length" calculated by this method is in fact the
"two"-norm (also know as the Euclidean norm) of a vector "C<$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!
=item *
C<$xn_vector = $matrix-E<gt>>C<solve_GSM($x0_vector,$b_vector,$epsilon);>
=item *
C<$xn_vector = $matrix-E<gt>>C<solve_SSM($x0_vector,$b_vector,$epsilon);>
=item *
C<$xn_vector = $matrix-E<gt>>C<solve_RM($x0_vector,$b_vector,$weight,$epsilon);>
In some cases it might not be practical or desirable to solve an
equation system "C<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 "C<x(t+1) = Phi( x(t) )>" which has
the property:
Phi(x) = x <==> A * x = b
We can define "C<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
This last step is true because
x[i] - ( a[i,1] x[1] + ... + a[i,i] x[i] + ... + a[i,n] x[n] ) + b[i]
is the same as
( -a[i,1] x[1] + ... + (1 - a[i,i]) x[i] + ... + -a[i,n] x[n] ) + b[i]
qed
Note that actually solving the equation system "C<A * x = b>" means
to calculate
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]
( run in 0.823 second using v1.01-cache-2.11-cpan-71847e10f99 )