=head1 NAME
Math::Orthonormalize - Gram-Schmidt Orthonormalization of vectors
=head1 SYNOPSIS
use Math::Orthonormalize qw(:all);
my @base_of_r_2 = (
[2, 1],
[1, 3]
);
my $vector = [1, 2, 3];
my @orthonormalized = orthonormalize(@base_of_r_2);
my @orthogonalized = orthogonalize(@base_of_r_2);
my $normalized = normalize($vector);
my $scaled = scale(2, $vector);
my $scalar = scalar_product($vector1, $vector2);
=head1 DESCRIPTION
Math::Orthonormalize offers subroutines to compute normalized or
non-normalized orthogonal bases of Euclidean vector spaces. That means:
Given a vector base of R^n, it computes a new base of R^n whose individual
vectors are all orthogonal. If those new base vectors all have a length of
1, the base is orthonormalized.
The module uses the Gram-Schmidt Algorithm.
=head2 EXPORT
No subroutines are exported by default, but the standart Exporter semantics are
in place, including the ':all' tag that imports all of the exportable
subroutines which are listed below.
=cut
package Math::Orthonormalize;
use 5.006;
use strict;
use warnings;
use Math::Symbolic qw/parse_from_string/;
use Carp;
our $VERSION = '1.00';
require Exporter;
our @ISA = qw(Exporter);
our %EXPORT_TAGS = ( 'all' => [ qw(
orthonormalize
orthogonalize
scalar_product
normalize
scale
) ] );
our @EXPORT_OK = ( @{ $EXPORT_TAGS{'all'} } );
our @EXPORT = ();
=head1 SUBROUTINES
=head2 orthonormalize
Takes any number (>1) of vectors (array refs of vector components) as argument
which form a base (that is, they are linearly independent) and returns
an orthogonalized and normalized base of the same vector space
(that is, n new array references).
=cut
sub orthonormalize {
my @vectors = @_;
croak "No arguments to orthogonalize"
if not @vectors;
croak "Arguments to orthogonalize must be array refs (vectors)"
if grep {ref($_) ne 'ARRAY'} @vectors;
my $dim = @{$vectors[0]};
croak "Vectors must have same dimension for orthogonalization"
if grep {@$_ != $dim} @vectors;
my @base = orthogonalize(@vectors);
@base = map {normalize($_)} @base;
return @base;
}
=head2 orthogonalize
Takes any number (>1) of vectors (array refs of vector components) as argument
which form a base (that is, they are linearly independent) and returns
an orthogonalized base of the same vector space (that is, n new array
references).
=cut
sub orthogonalize {
my @vectors = @_;
croak "No arguments to orthogonalize"
if not @vectors;
croak "Arguments to orthogonalize must be array refs (vectors)"
if grep {ref($_) ne 'ARRAY'} @vectors;
my $dim = @{$vectors[0]};
croak "Vectors must have same dimension for orthogonalization"
if grep {@$_ != $dim} @vectors;
my @newbase;
push @newbase, $vectors[0];
my @squares;
push @squares, scalar_product($vectors[0], $vectors[0]) if @vectors > 1;
foreach my $i (1..$#vectors) {
my $new = $vectors[$i];
foreach my $j (0..$#newbase) {
my $vec = scale(
scalar_product($vectors[$i], $newbase[$j])
/ $squares[$j],
$newbase[$j]
);
@$new = map {$_ -= shift @$vec} @$new;
}
push @newbase, $new;
push @squares, scalar_product($new, $new) if $i < $#vectors;
}
return @newbase;
}
=head2 normalize
Normalizes a vector. That is, it changes the vector length to 1 without
changing the vector's direction.
Takes an array reference with the vector components as argument and returns
a new array reference containing the normalized vector components.
=cut
sub normalize {
my $v = shift;
croak "Argument to normalize() must be an array ref (vector)"
if not ref($v) eq 'ARRAY';
croak "Cannot normalize 0-dimensional vectors"
if @$v == 0;
my $sum = $v->[0]**2;
$sum += $v->[$_]**2 for 1..$#$v;
$sum = sqrt($sum);
croak "Cannot normalize 0-vector"
if $sum == 0;
return [map {$_ / $sum} @$v];
}
=head2 scale
Takes a scalar and a vector (array reference of vector components) as
arguments. Multiplies every component of the vector by the specified
scalar and returns a new array reference containing the scaled vector
components.
=cut
sub scale {
croak "scale() takes a scalar and a vector (ary ref) as arguments"
if not @_ == 2 or not ref($_[1]) eq 'ARRAY';
croak "Cannot handle 0-dimensional vectors"
if @{$_[1]} == 0;
my $new = [];
foreach (@{$_[1]}) {
push @$new, $_[0] * $_;
}
return $new;
}
=head2 scalar_product
Computes the scalar product of two vectors. Expects two array references
with vector components (same number of components) as argument and
returns their scalar product.
=cut
sub scalar_product {
croak "Invalid number of arguments to scalar_product()"
if not @_ == 2;
my ($v1, $v2) = @_;
croak "Cannot compute the scalar product of vectors of different ",
"dimensions"
if @$v1 != @$v2;
croak "Cannot deal with 0-dimensional vectors"
if @$v1 == 0;
my $sum = $v1->[0] * $v2->[0];
return $sum if @$v1 == 1;
$sum += $v1->[$_] * $v2->[$_] for 1..$#$v1;
return $sum;
}
1;
__END__
=head1 AUTHOR
Steffen Mueller, orthonormalize-module at steffen-mueller dot net
=head1 SEE ALSO
(German) Merziger, Wirth: "Repetitorium der Höheren Mathematik" (Binomi, 1999)
You may find the current versions of this module at http://steffen-mueller.net/
or on CPAN.
=head1 COPYRIGHT AND LICENSE
Copyright (C) 2004-2005 by Steffen Mueller
This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.6.1 or,
at your option, any later version of Perl 5 you may have available.
=cut