Implements the data type "matrix of reals" (and consequently also
"vector of reals")
Implements the data type "matrix of reals", which can be used almost
like any other basic Perl type thanks to OPERATOR OVERLOADING, i.e.,
does what you would like it to do (a matrix multiplication).
Also features many important operations and methods: matrix norm,
matrix transposition, matrix inverse, determinant of a matrix, order
and numerical condition of a matrix, scalar product of vectors, vector
product of vectors, vector length, projection of row and column vectors,
a comfortable way for reading in a matrix from a file, the keyboard or
your code, and many more.
Allows to solve linear equation systems using an efficient algorithm
known as "L-R-decomposition" and several approximative (iterative) methods.
Features an implementation of Kleene's algorithm to compute the minimal
costs for all paths in a graph with weighted edges (the "weights" being
the costs associated with each edge).
-
use MatrixReal1;
Makes the methods and overloaded operators of this module available
to your program.
-
use MatrixReal1 qw(min max);
-
use MatrixReal1 qw(:all);
Use one of these two variants to import (all) the functions which the module
offers for export; currently these are "min()" and "max()".
-
$new_matrix = new MatrixReal1($rows,$columns);
The matrix object constructor method.
Note that this method is implicitly called by many of the other methods
in this module!
-
$new_matrix = MatrixReal1->new($rows,$columns);
An alternate way of calling the matrix object constructor method.
-
$new_matrix = $some_matrix->new($rows,$columns);
Still another way of calling the matrix object constructor method.
Matrix "$some_matrix" is not changed by this in any way.
-
$new_matrix = MatrixReal1->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).
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 = MatrixReal1->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 = MatrixReal1->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 = MatrixReal1->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 = MatrixReal1->new_from_string($string); };
if ($@)
{
print "$@";
# ...
# (error handling)
}
else
{
# continue...
}
or as follows:
eval { $matrix = MatrixReal1->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 =
MatrixReal1->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:
MatrixReal1::new_from_string(): syntax error in input string
MatrixReal1::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:
MatrixReal1::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.
-
$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".
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"!).
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".
Matrix "$some_matrix" is not changed by this in any way.
-
$row_vector = $matrix->row($row);
This is a projection method which returns an object reference to
a NEW matrix (which in fact is a (row) vector since it has only
one row) to which row number "$row" of matrix "$matrix" has
already been copied.
Matrix "$matrix" is not changed by this in any way.
-
$column_vector = $matrix->column($column);
This is a projection method which returns an object reference to
a NEW matrix (which in fact is a (column) vector since it has
only one column) to which column number "$column" of matrix
"$matrix" has already been copied.
Matrix "$matrix" is not changed by this in any way.
-
$matrix->zero();
Assigns a zero to every element of the matrix "$matrix", i.e.,
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 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!)
-
$matrix->assign($row,$column,$value);
Explicitly assigns a value "$value" to a single element of the
matrix "$matrix", located in row "$row" and column "$column",
thereby replacing the value previously stored there.
-
$value = $matrix->element($row,$column);
Returns the value of a specific element of the matrix "$matrix",
located in row "$row" and column "$column".
-
($rows,$columns) = $matrix->dim();
Returns a list of two items, representing the number of rows
and columns 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.
-
$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.
-
$matrix1->transpose($matrix2);
Calculates the transposed matrix of matrix "$matrix2" 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.
Transposition is a symmetry operation: imagine you rotate the matrix
along the axis of its main diagonal (going through elements (1,1),
(2,2), (3,3) and so on) by 180 degrees.
Another way of looking at it is to say that rows and columns are
swapped. In fact the contents of element (i,j) are swapped
with those of element (j,i).
Note that (especially for vectors) it makes a big difference if you
have 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!
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!
Hint: throughout this module, whenever a vector is explicitly required
for input, a COLUMN vector is expected!
-
$matrix1->add($matrix2,$matrix3);
Calculates the sum of matrix "$matrix2" and matrix "$matrix3"
and stores the result in matrix "$matrix1" (which 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.
-
$matrix1->subtract($matrix2,$matrix3);
Calculates the difference of matrix "$matrix2" minus matrix "$matrix3"
and stores the result in matrix "$matrix1" (which 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.
-
$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".
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.
-
$minimum = MatrixReal1::min($number1,$number2);
Returns the minimum of the two numbers "number1" and "number2".
-
$minimum = MatrixReal1::max($number1,$number2);
Returns the maximum of the two numbers "number1" and "number2".
-
$minimal_cost_matrix = $cost_matrix->kleene();
Copies the matrix "$cost_matrix" (which has to be quadratic!) to
a new matrix of the same size (i.e., "clones" the input matrix) and
applies Kleene's algorithm to it.
See the Math::Kleene(3) manpage for more details about this algorithm!
The method returns an object reference to the new matrix.
Matrix "$cost_matrix" is not changed by this method in any way.
-
($norm_matrix,$norm_vector) = $matrix->normalize($vector);
This method is used to improve the numerical stability when solving
linear equation systems.
Suppose you have a matrix "A" and a vector "b" and you want to find
out a vector "x" so that A * x = b, i.e., the vector "x" which
solves the equation system represented by the matrix "A" and the
vector "b".
Applying this method to the pair (A,b) yields a pair (A',b') where
each row has been divided by (the absolute value of) the greatest
coefficient appearing in that row. So this coefficient becomes equal
to "1" (or "-1") in the new pair (A',b') (all others become smaller
than one and greater than minus one).
Note that this operation does not change the equation system itself
because the same division is carried out on either side of the equation
sign!
The method requires a quadratic (!) matrix "$matrix" and a vector
"$vector" for input (the vector must be a column vector with the same
number of rows as the input matrix) and returns a list of two items
which are object references to a new matrix and a new vector, in this
order.
The output matrix and vector are clones of the input matrix and vector
to which the operation explained above has been 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 to solve
equation systems (explained immediately below following this method):
Consider the following little program:
#!perl -w
use MatrixReal1 qw(new_from_string);
$A = MatrixReal1->new_from_string(<<"MATRIX");
[ 1 2 3 ]
[ 5 7 11 ]
[ 23 19 13 ]
MATRIX
$b = MatrixReal1->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 ]
You can see that in the second example (where "normalize()" has been used),
the result is "better", i.e., more accurate!
-
$LR_matrix = $matrix->decompose_LR();
This method is needed to solve linear equation systems.
Suppose you have a matrix "A" and a vector "b" and you want to find
out a vector "x" so that A * x = b, i.e., the vector "x" which
solves the equation system represented by the matrix "A" and the
vector "b".
You might also have a matrix "A" and a whole bunch of different
vectors "b1".."bk" for which you need to find vectors "x1".."xk"
so that A * xi = bi, for i=1..k.
Using Gaussian transformations (multiplying a row or column with
a factor, swapping two rows or two columns and adding a multiple
of one row or column to another), it is 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 of 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, as follows:
[ 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 (these are 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 information
loss! 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, that "LR" and "L * R" are not the same!!!
Note also that for the 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 help 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, simply divide "b[i]" by "A[i,i]"
to get "x[i]"!).
To find the solution to our problem "A * x = b", we divide this
problem in 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 is 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 have many vectors "b1".."bk"
for which you are searching solutions to A * xi = bi - this scheme
is much more efficient than a straightforward, "brute force" approach.
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, i.e.,
a matrix that has 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 concerning the hidden
permutations of its rows and columns).
However, as soon as you are applying any method that alters the contents
of the matrix, its "magical" properties are stripped off, and the matrix
immediately reverts to an "ordinary" matrix (with the values it just happens
to contain at that moment, be they 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 returns a list of three items if a solution exists or an
empty list otherwise (!).
Therefore, you should always use this method like this:
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 items returned are: the dimension "$dimension" of the solution
space (which is zero if only one solution exists, one if the solution is
a straight line, two if the solution is a plane, and so on), the solution
vector "$x_vector" (which is the vector "x" of your equation system
A * x = b) and a matrix "$base_matrix" representing a base of the
solution space (a set of vectors which put 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.
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!
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.
Therefore, 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.
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).
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.
-
$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!) vector "$vector".
Note that the "length" calculated by this method is in fact the
"two"-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 "MatrixReal1::vector_norm(): vector is not a column vector"
unless ($cols == 1);
croak "MatrixReal1::vector_norm(): norm index must be > 0"
unless ($n > 0);
croak "MatrixReal1::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
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 "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]
There is one major restriction, though: a fix-point iteration is
guaranteed to converge only if the first derivative of the iteration
function has an absolute value less than one in an area around the
point "x(*)" for which "Phi( x(*) ) = x(*)" is to be true, and
if the start vector "x(0)" lies within that area!
This is best verified grafically, which unfortunately is impossible
to do in this textual documentation!
See literature on Numerical Analysis for details!
In our case, this restriction translates to the following three conditions:
There must exist a norm so that the norm of the matrix of the iteration
function, ( En - A ), has a value less than one, the matrix "A" may
not have any zero value on its main diagonal and the initial vector
"x(0)" must be "good enough", i.e., "close enough" to the solution
"x(*)".
(Remember school math: the first derivative of a straight line given by
"y = a * x + b" is "a"!)
The three methods expect a (quadratic!) matrix "$matrix" as their
first argument, a start vector "$x0_vector", a vector "$b_vector"
(which is the vector "b" in your equation system "A * x = b"), in the
case of the "Relaxation Method" ("RM"), a real number "$weight" best
between zero and two, and finally an error limit (real number) "$epsilon".
(Note that 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 conditions of the three
conditions listed above and return an empty list if these conditions
are not fulfilled.
Therefore, you should always test their return value using some
code like:
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 iterate until abs( Phi(x) - x ) < epsilon.
(Beware that theoretically, infinite loops might result if the starting
vector is too far "off" the solution! In practice, this shouldn't be
a problem. Anyway, you can always press <ctrl-C> if you think that the
iteration takes too long!)
The difference between the three methods is the following:
In the "Global Step Method" ("GSM"), the new vector "x(t+1)"
(called "y" here) is calculated from the vector "x(t)"
(called "x" here) according to the 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]
In the "Single Step Method" ("SSM"), the components of the vector
"x(t+1)" which have already been calculated are 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]
In the "Relaxation method" ("RM"), the components of the vector
"x(t+1)" are calculated by "mixing" old and new value (like
cold and hot water), and the weight "$weight" determines the
"aperture" of both the "hot water tap" as well as of the "cold
water tap", according to the formula:
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]
Note that the weight "$weight" should be greater than zero and
less than two (!).
The three methods are supposed to be of different efficiency.
Experiment!
Remember that in most cases, it is probably advantageous to first
"normalize()" your equation system prior to solving it!