From e9319c22656963ef907c6f201ab2042bc604053b Mon Sep 17 00:00:00 2001 From: Alex Jordan Date: Mon, 1 Jun 2026 20:39:28 -0700 Subject: [PATCH] fuzzy test for row equivalency --- lib/Value/Matrix.pm | 95 +++++++++++++++++++++++++++++++++++++++++ t/math_objects/matrix.t | 44 +++++++++++++++++++ 2 files changed, 139 insertions(+) diff --git a/lib/Value/Matrix.pm b/lib/Value/Matrix.pm index 5b1892cdd..8f720a308 100644 --- a/lib/Value/Matrix.pm +++ b/lib/Value/Matrix.pm @@ -627,6 +627,101 @@ sub isRREF { return 1; } +=head3 C + +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'); diff --git a/t/math_objects/matrix.t b/t/math_objects/matrix.t index 32dfe3d45..5b2f2f70c 100644 --- a/t/math_objects/matrix.t +++ b/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 ] ]);