fuzzy test for row equivalency

Author: Alex-Jordan

lib/Value/Matrix.pm

@@ -631,6 +631,101 @@ sub isRREF {
 	return 1;
 }
 
+=head3 C<isREQ>
+
+Fuzzy, probabilistic check for if two matrices are row equivalent.
+
+Usage:
+
+   $A = Matrix([3, 1], [1, 0.3333]);
+   $B = Matrix([1, 1/3], [0, 0]);
+   $A->isREQ($B);   # true
+
+=cut
+
+sub isREQ {
+	my ($l, $r) = @_;
+	my @ldim = $l->dimensions;
+	my @rdim = $r->dimensions;
+	Value::Error('Cannot compare row equivalency of matrices of different degree') unless scalar @ldim == scalar @rdim;
+	for my $i (0 .. $#ldim) {
+		Value::Error('Cannot compare row equivalency because matrices differ in the '
+				. $l->NameForNumber($i + 1)
+				. ' dimension')
+			unless $ldim[$i] == $rdim[$i];
+	}
+	return 1 if $l->isZero && $r->isZero;
+	return 0 if $l->isZero && !$r->isZero || !$l->isZero && $r->isZero;
+
+	if (scalar @ldim == 1) {
+		# simply need to determine if rows are parallel
+		my ($lk, $rk);
+		for my $i (1 .. $ldim[0]) {
+			if ($l->element($i) != 0) {
+				$lk = $l->element($i);
+				$rk = $r->element($i);
+				last;
+			}
+		}
+		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;
+		}
+
+		return 0 if scalar @lgs != scalar @rgs;
+
+		# project each row from $rrows onto @lgs;
+		# if the complement is nonzero, the row spaces disagree
+		for my $v (@rrows) {
+			my $proj = $r->new((0) x $ldim[1]);
+			for my $w (@lgs) {
+				$proj += ($v * $w) * $w;
+			}
+			return 0 unless $v == $proj;
+		}
+		# and vice versa
+		for my $v (@lrows) {
+			my $proj = $r->new((0) x $rdim[1]);
+			for my $w (@rgs) {
+				$proj += ($v * $w) * $w;
+			}
+			return 0 unless $v == $proj;
+		}
+
+		# if we got this far, the row spaces are (fuzzy) equal
+		# so the matrices are row equivalent
+		return 1;
+	} else {
+		for my $i (1 .. $ldim[0]) {
+			return 0 unless $l->new($l->{data}[ $i - 1 ])->isREQ($r->new($r->{data}[ $i - 1 ]));
+		}
+		return 1;
+	}
+}
+
 sub _isNumber {
 	my $n = shift;
 	return Value::isNumber($n) || Value::classMatch($n, 'Fraction');

t/math_objects/matrix.t

@@ -180,6 +180,50 @@ subtest 'Test if Matrix is in (R)REF' => sub {
 	ok !$B4->isRREF, "$B4 is not in RREF";
 };
 
+subtest 'Check if two matrices are (fuzzy) row equivalent' => sub {
+	my $A1 = Matrix(1, 2, 3);
+	my $A2 = Matrix([ 1, 2, 3 ], [ 4, 5, 6 ]);
+	my $A3 = Matrix([ [ 1, 2 ], [ 3, 4 ] ], [ [ 5, 6 ], [ 7, 8 ] ]);
+	my $B1 = Matrix(2,        4,        6);
+	my $B2 = Matrix([ 1, 2 ], [ 3, 4 ], [ 5, 6 ]);
+	my $C1 = Matrix(0,        4,        6);
+	my $Z1 = Matrix([ 0, 0, 0 ], [ 0, 0, 0 ]);
+	my $Z2 = Matrix([ 0, 0, 0 ], [ 0, 0, 10**-14 ]);
+
+	my $M1 = Matrix([ 3, 1 ],      [ 1, 1 / 3 ]);
+	my $M2 = Matrix([ 3, 1 ],      [ 0, 0 ]);
+	my $M3 = Matrix([ 1, 0.3333 ], [ 0, 0 ]);
+	my $M4 = Matrix([ 1, 0.33 ],   [ 0, 0 ]);
+	my $M5 = Matrix([ 6, 2 ],      [ 9, 3 ]);
+	my $M6 = Matrix([ 3, 1 ],      [ 1, 3 ]);
+
+	my $N1 = Matrix([ [ 1, 2 ], [ 2, 4 ] ], [ [ 1, 2 ], [ 3, 4 ] ]);
+	my $N2 = Matrix([ [ 3, 6 ], [ 0, 0 ] ], [ [ 1, 0 ], [ 0, 1 ] ]);
+	my $N3 = Matrix([ [ 3, 6 ], [ 1, 0 ] ], [ [ 1, 0 ], [ 0, 1 ] ]);
+
+	ok $Z1->isREQ($Z2),  'Zero matrices are row equivalent';
+	ok !$Z1->isREQ($A2), 'A zero matrix is not row equivalent to a nonzero matrix';
+	ok $A1->isREQ($B1),  'Parallel degree 1 matrices are row equivalent';
+	ok !$C1->isREQ($B1), 'Non-parallel degree 1 matrices are not row equivalent';
+	ok $M2->isREQ($M5),  'Row equivalent matrices are row equivalent';
+	ok $M2->isREQ($M1),  'Row equivalent matrices are row equivalent, even with some machine rounding';
+	ok $M2->isREQ($M3),  'Row equivalent matrices are row equivalent, even with student rounding';
+	ok !$M2->isREQ($M4), 'Matrices are not row equivalent if rounding is too much';
+	ok !$M2->isREQ($M6), 'Non-row equivalent matrices are not row equivalent';
+	ok $N1->isREQ($N2),  'Row equivalent matrices are row equivalent, even at degree above 2';
+	ok !$N1->isREQ($N3), 'Non-row equivalent matrices are not row equivalent, even at degree above 2';
+
+	like dies {
+		$A2->isREQ($A3);
+	}, qr/Cannot compare row equivalency of matrices of different degree/,
+		'Test that error is thrown for comparing matrices of differing degrees';
+	like dies {
+		$A2->isREQ($B2);
+	}, qr/Cannot compare row equivalency because matrices differ in the first dimension/,
+		'Test that error is thrown for comparing matrices of differing dimensioon';
+
+};
+
 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 ] ]);