normalized gram schmit method

Author: Alex-Jordan

lib/Value/Matrix.pm

@@ -669,36 +669,14 @@ sub isREQ {
 		}
 		return $lk * $r == $rk * $l;
 	} elsif (scalar @ldim == 2) {
-		# get Gram-Schmidt basis for each row space
-		my @rrows = map { $r->new($_) } $r->value;
-		@rrows = grep { !$_->isZero } @rrows;
-		@rrows = map  { $_ / ($_ * $_)**0.5 } @rrows;
-		my @rgs = ($rrows[0]);
-		for my $i (1 .. $#rrows) {
-			my $proj = $r->new((0) x $rdim[1]);
-			for my $v (@rgs) {
-				$proj += ($rrows[$i] * $v) * $v;
-			}
-			my $gs = $rrows[$i] - $proj;
-			push @rgs, $gs / ($gs * $gs)**0.5 unless $rrows[$i] == $proj;
-		}
-		my @lrows = map { $l->new($_) } $l->value;
-		@lrows = grep { !$_->isZero } @lrows;
-		@lrows = map  { $_ / ($_ * $_)**0.5 } @lrows;
-		my @lgs = ($lrows[0]);
-		for my $i (1 .. $#lrows) {
-			my $proj = $r->new((0) x $ldim[1]);
-			for my $v (@lgs) {
-				$proj += ($lrows[$i] * $v) * $v;
-			}
-			my $gs = $lrows[$i] - $proj;
-			push @lgs, $gs / ($gs * $gs)**0.5 unless $lrows[$i] == $proj;
-		}
+		my @rgs = $r->NGS;
+		my @lgs = $l->NGS;
 
 		return 0 if scalar @lgs != scalar @rgs;
 
 		# project each row from $rrows onto @lgs;
 		# if the complement is nonzero, the row spaces disagree
+		my @rrows = map { $r->new($_) } $r->value;
 		for my $v (@rrows) {
 			my $proj = $r->new((0) x $ldim[1]);
 			for my $w (@lgs) {
@@ -707,6 +685,7 @@ sub isREQ {
 			return 0 unless $v == $proj;
 		}
 		# and vice versa
+		my @lrows = map { $l->new($_) } $l->value;
 		for my $v (@lrows) {
 			my $proj = $r->new((0) x $rdim[1]);
 			for my $w (@rgs) {
@@ -726,6 +705,45 @@ sub isREQ {
 	}
 }
 
+=head3 C<NGS>
+
+Normalized Gram Schmit basis for a given list of vectors. This uses "fuzzy" equivalency to avoid machine rounding
+issues. If called as a method on a Matrix, the rows of that Matrix are used as the list of vectors.
+
+Usage:
+
+   $A = Matrix([1, 2, 2], [2, 2, 1]);
+   $A->NGS                                   # returns [Matrix(1, 2, 2)/3, Matrix(10, 2, -7)/sqrt(153)]
+
+   Value::Matrix->NGS([1, 2, 2], [2, 2, 1])  #returns the same as above
+
+=cut
+
+sub NGS {
+	my ($self, @rows) = @_;
+	@rows = $self->value if !@rows && ref($self) eq 'Value::Matrix';
+	Value::Error('You must provide vectors to apply Gram Schmit to') if !@rows;
+	@rows = map  { Value::Matrix->new($_) } @rows;
+	@rows = grep { !$_->isZero } @rows;
+	Value::Error('You must provide at least one nonzero row for Gram Schmit') unless @rows;
+	my $n = ($rows[0]->dimensions)[0];
+	for my $r (@rows) {
+		Value::Error('Rows provided for Gram Schmit should not be nested arrays')     unless $r->degree == 1;
+		Value::Error('All rows provided for Gram Schmit should have the same length') unless ($r->dimensions)[0] == $n;
+	}
+	@rows = map { $_ / ($_ * $_)**0.5 } @rows;
+	my @gs = ($rows[0]);
+	for my $i (1 .. $#rows) {
+		my $proj = Value::Matrix->new((0) x $n);
+		for my $v (@gs) {
+			$proj += ($rows[$i] * $v) * $v;
+		}
+		my $gs = $rows[$i] - $proj;
+		push @gs, $gs / ($gs * $gs)**0.5 unless $rows[$i] == $proj;
+	}
+	return @gs;
+}
+
 sub _isNumber {
 	my $n = shift;
 	return Value::isNumber($n) || Value::classMatch($n, 'Fraction');

t/math_objects/matrix.t

@@ -224,6 +224,33 @@ subtest 'Check if two matrices are (fuzzy) row equivalent' => sub {
 
 };
 
+subtest 'Normalized Gram Schmit' => sub {
+	my @A      = ([ 1, 2, 2 ], [ 2, 2, 1 ]);
+	my $A      = Matrix(@A);
+	my @NGSA   = $A->NGS;
+	my @NGSVMA = Value::Matrix->NGS(@A);
+
+	is "@NGSA", "[0.333333,0.666667,0.666667] [0.808452,0.16169,-0.565916]",
+		'Test that NGS finds a correct normalized GS basis';
+	is "@NGSVMA", "[0.333333,0.666667,0.666667] [0.808452,0.16169,-0.565916]",
+		'Test that NGS finds a correct normalized GS basis';
+
+	like dies {
+		Value::Matrix->NGS();
+	}, qr/You must provide vectors to apply Gram Schmit to/, 'Test that you must pass something as an argument';
+	like dies {
+		Value::Matrix->NGS([ 0, 0 ], [ 0, 10**-16 ]);
+	}, qr/You must provide at least one nonzero row for Gram Schmit/,
+		'Test that you must pass at least one nonzero row';
+	like dies {
+		Value::Matrix->NGS([ [ 1, 0 ], [ 1, 0 ] ]);
+	}, qr/Rows provided for Gram Schmit should not be nested arrays/,
+		'Test that rows are rows, not nested matrices';
+	like dies {
+		Value::Matrix->NGS([ 1, 0 ], [ 1, 0, 0 ]);
+	}, qr/All rows provided for Gram Schmit should have the same length/, 'Test that rows have same length';
+};
+
 subtest 'Transpose a Matrix' => sub {
 	my $A = Matrix([ [ 1, 2, 3, 4 ], [ 5, 6, 7, 8 ], [ 9, 10, 11, 12 ] ]);
 	my $B = Matrix([ [ 1, 5, 9 ], [ 2, 6, 10 ], [ 3, 7, 11 ], [ 4, 8, 12 ] ]);