Math-Matrix
view release on metacpan or search on metacpan
lib/Math/Matrix.pm view on Meta::CPAN
next;
}
croak "unknown parameter '$param'";
}
if ($debug) {
printf "\n";
printf "max_iter = %24d\n", $max_iter;
printf "rel_tol = %24.15e\n", $rel_tol;
printf "abs_tol = %24.15e\n", $abs_tol;
}
my $y_norm = _hypot(map { @$_ } @$y);
my $x = $y -> mldiv($A);
my $x_best;
my $iter_best;
my $abs_err_best;
my $rel_err_best;
for (my $iter = 1 ; ; $iter++) {
# Compute the residuals.
my $r = $A -> mmuladd($x, -$y);
# Compute the errors.
my $r_norm = _hypot(map @$_, @$r);
my $abs_err = $r_norm;
my $rel_err = $y_norm == 0 ? $r_norm : $r_norm / $y_norm;
if ($debug) {
printf "\n";
printf "iter = %24d\n", $iter;
printf "r_norm = %24.15e\n", $r_norm;
printf "y_norm = %24.15e\n", $y_norm;
printf "abs_err = %24.15e\n", $abs_err;
printf "rel_err = %24.15e\n", $rel_err;
}
# See if this is the first round or we have an new all-time best.
if ($iter == 1 ||
$abs_err < $abs_err_best ||
$rel_err < $rel_err_best)
{
$x_best = $x;
$iter_best = $iter;
$abs_err_best = $abs_err;
$rel_err_best = $rel_err;
}
if ($abs_err_best <= $abs_tol || $rel_err_best <= $rel_tol) {
last;
} else {
# If we still haven't got the desired result, but have reached
# the maximum number of iterations, display a warning.
if ($iter == $max_iter) {
carp "mldiv() stopped because the maximum number of",
" iterations (max. iter = $max_iter) was reached without",
" converging to any of the desired tolerances (",
"rel_tol = ", $rel_tol, ", ",
"abs_tol = ", $abs_tol, ").",
" The best iterate (iter. = ", $iter_best, ") has",
" a relative residual of ", $rel_err_best, " and",
" an absolute residual of ", $abs_err_best, ".";
last;
}
}
# Compute delta $x.
my $d = $r -> mldiv($A);
# Compute the improved solution $x.
$x -= $d;
}
return $x_best, $rel_err_best, $abs_err_best, $iter_best if wantarray;
return $x_best;
}
# If A is an M-by-M, compute A\y directly.
croak "mldiv(): sizes don't match" unless $y -> nrow() == $n;
# Create the augmented matrix.
my $x = $A -> catcol($y);
# Perform forward elimination.
my ($rowperm, $colperm);
eval { ($x, $rowperm, $colperm) = $x -> felim_fp() };
croak "mldiv(): matrix is singular or close to singular" if $@;
# Perform backward substitution.
eval { $x = $x -> bsubs() };
croak "mldiv(): matrix is singular or close to singular" if $@;
# Remove left half to keep only the augmented matrix.
$x = $x -> splicecol(0, $n);
# Reordering the rows is only necessary when full (complete) pivoting is
# used above. If partial pivoting is used, this reordeing could be skipped,
# but it executes so fast that it causes no harm to do it anyway.
@$x[ @$colperm ] = @$x;
return $x;
}
=pod
=item sldiv()
Scalar (left) division.
lib/Math/Matrix.pm view on Meta::CPAN
=item C<~>
Transpose.
$B = ~$A; # $B is the transpose of $A
=item C<abs>
Absolute value.
$B = abs $A; # $B contains absolute values of $A
=item C<int>
Truncate to integer.
$B = int $A; # $B contains only integers
=back
=head1 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.
=head2 Background
The linear system of equations
$A * $x = $y
can be solved for C<$x> with
$x = $y -> mldiv($A);
Ideally C<$A * $x> should equal C<$y>, but due to numerical errors, this is not
always the case. The following illustrates how to improve the solution C<$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.
=head2 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.
=head2 Parameters
=over 4
=item MaxIter
The maximum number of iterations to perform. The value must be a positive
integer. The default is 20.
=item RelTol
The limit for the relative error. The value must be a non-negative. The default
value is 1e-19 when perl is compiled with long doubles or quadruple precision,
and 1e-9 otherwise.
=item AbsTol
The limit for the absolute error. The value must be a non-negative. The default
value is 0.
=item Debug
If this parameter does not affect when the algorithm terminates, but when set to
non-zero, some information is displayed at each step.
=back
=head2 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 C<< $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 C<< $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 C<< $dx = $res -> mldiv($A); >> is
$dx = [[ -8.592098303124e-16, -2.86724066474914e-15 ],
( run in 0.538 second using v1.01-cache-2.11-cpan-71847e10f99 )