From b83c694b64f5110f32ccd734d73e2e3b3a3bf76d Mon Sep 17 00:00:00 2001 From: Alex Jordan Date: Mon, 1 Jun 2026 20:39:28 -0700 Subject: [PATCH 1/2] 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 277d92f83..f7257ebe3 100644 --- a/lib/Value/Matrix.pm +++ b/lib/Value/Matrix.pm @@ -631,6 +631,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 c793daa31..dadaf3472 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 ] ]); From 3ffd9052423503d624fe3f97062c77075da1668e Mon Sep 17 00:00:00 2001 From: Alex Jordan Date: Sat, 20 Jun 2026 14:21:45 -0700 Subject: [PATCH 2/2] normalized gram schmit method --- lib/Value/Matrix.pm | 92 ++++++++++++++++++++++++++++++----------- t/math_objects/matrix.t | 38 +++++++++++++++++ 2 files changed, 105 insertions(+), 25 deletions(-) diff --git a/lib/Value/Matrix.pm b/lib/Value/Matrix.pm index f7257ebe3..e14d2ce4f 100644 --- a/lib/Value/Matrix.pm +++ b/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,69 @@ sub isREQ { } } +=head3 C + +Normalized Gram Schmidt 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. If called as a +method with the option C<< cols => 1 >>, then the columns are used as the list of vectors. + +In array context, returns an array of 1D Matrix objects. In scalar context, returns a Matrix with those rows. + +If C<< cols => 1 >> was used, return the Matrix objects in column form. + + +Usage: + + $A = Matrix([1, 2, 2], [2, 2, 1]); + + @X = $A->NGS # @X is (Matrix(1, 2, 2)/3, Matrix(10, 2, -7)/sqrt(153)) + @X = Value::Matrix->NGS([1, 2, 2], [2, 2, 1]) # @X is the same as above + $X = $A->NGS # $X is Matrix([1/3, 2/3, 2/3], [10/sqrt(153), 2/sqrt(153), -7/sqrt(153)]) + $X = Value::Matrix->NGS([1, 2, 2], [2, 2, 1]) # $X is the same as above + + @X = $A->NGS(cols => 1) # @X is (Matrix([1/sqrt(5)], [2/sqrt(5)]), Matrix([2/sqrt(5)], [-1/sqrt(5)])) + $X = $A->NGS(cols => 1) # $X is Matrix([1/sqrt(5), 2/sqrt(5)], [2/sqrt(5), -1/sqrt(5)]) + +=cut + +sub NGS { + my ($self, @args) = @_; + my %options; + my @rows; + if (ref($self) eq 'Value::Matrix') { + %options = @args; + @rows = $options{cols} ? $self->transpose->value : $self->value; + } else { + @rows = @args; + } + Value::Error('You must provide vectors to apply Gram Schmidt 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 Schmidt') unless @rows; + my $n = ($rows[0]->dimensions)[0]; + for my $r (@rows) { + Value::Error('Rows provided for Gram Schmidt should not be nested arrays') unless $r->degree == 1; + Value::Error('All rows provided for Gram Schmidt 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 + $options{cols} + ? wantarray + ? map { $_->transpose } @gs + : Value::Matrix->new(@gs)->transpose + : wantarray ? @gs + : Value::Matrix->new(@gs); +} + 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 dadaf3472..27dbf4811 100644 --- a/t/math_objects/matrix.t +++ b/t/math_objects/matrix.t @@ -224,6 +224,44 @@ subtest 'Check if two matrices are (fuzzy) row equivalent' => sub { }; +subtest 'Normalized Gram Schmidt' => sub { + my @A = ([ 1, 2, 2 ], [ 2, 2, 1 ]); + my $A = Matrix(@A); + my @NGSA = $A->NGS; + my @NGSAc = $A->NGS(cols => 1); + my @NGSVMA = Value::Matrix->NGS(@A); + my $NGSA = $A->NGS; + my $NGSAc = $A->NGS(cols => 1); + 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 "@NGSAc", "[[0.447214],[0.894427]] [[0.894427],[-0.447214]]", + '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'; + is "$NGSA", "[[0.333333,0.666667,0.666667],[0.808452,0.16169,-0.565916]]", + 'Test that NGS finds a correct normalized GS basis'; + is "$NGSAc", "[[0.447214,0.894427],[0.894427,-0.447214]]", '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 Schmidt 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 Schmidt/, + 'Test that you must pass at least one nonzero row'; + like dies { + Value::Matrix->NGS([ [ 1, 0 ], [ 1, 0 ] ]); + }, qr/Rows provided for Gram Schmidt 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 Schmidt 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 ] ]);