Skip to content

Commit 3ffd905

Browse files
Alex-JordanAlex-Jordan
committed
normalized gram schmit method
1 parent b83c694 commit 3ffd905

2 files changed

Lines changed: 105 additions & 25 deletions

File tree

‎lib/Value/Matrix.pm‎

Lines changed: 67 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -669,36 +669,14 @@ sub isREQ {
669669
}
670670
return $lk * $r == $rk * $l;
671671
} elsif (scalar @ldim == 2) {
672-
# get Gram-Schmidt basis for each row space
673-
my @rrows = map { $r->new($_) } $r->value;
674-
@rrows = grep { !$_->isZero } @rrows;
675-
@rrows = map { $_ / ($_ * $_)**0.5 } @rrows;
676-
my @rgs = ($rrows[0]);
677-
for my $i (1 .. $#rrows) {
678-
my $proj = $r->new((0) x $rdim[1]);
679-
for my $v (@rgs) {
680-
$proj += ($rrows[$i] * $v) * $v;
681-
}
682-
my $gs = $rrows[$i] - $proj;
683-
push @rgs, $gs / ($gs * $gs)**0.5 unless $rrows[$i] == $proj;
684-
}
685-
my @lrows = map { $l->new($_) } $l->value;
686-
@lrows = grep { !$_->isZero } @lrows;
687-
@lrows = map { $_ / ($_ * $_)**0.5 } @lrows;
688-
my @lgs = ($lrows[0]);
689-
for my $i (1 .. $#lrows) {
690-
my $proj = $r->new((0) x $ldim[1]);
691-
for my $v (@lgs) {
692-
$proj += ($lrows[$i] * $v) * $v;
693-
}
694-
my $gs = $lrows[$i] - $proj;
695-
push @lgs, $gs / ($gs * $gs)**0.5 unless $lrows[$i] == $proj;
696-
}
672+
my @rgs = $r->NGS;
673+
my @lgs = $l->NGS;
697674

698675
return 0 if scalar @lgs != scalar @rgs;
699676

700677
# project each row from $rrows onto @lgs;
701678
# if the complement is nonzero, the row spaces disagree
679+
my @rrows = map { $r->new($_) } $r->value;
702680
for my $v (@rrows) {
703681
my $proj = $r->new((0) x $ldim[1]);
704682
for my $w (@lgs) {
@@ -707,6 +685,7 @@ sub isREQ {
707685
return 0 unless $v == $proj;
708686
}
709687
# and vice versa
688+
my @lrows = map { $l->new($_) } $l->value;
710689
for my $v (@lrows) {
711690
my $proj = $r->new((0) x $rdim[1]);
712691
for my $w (@rgs) {
@@ -726,6 +705,69 @@ sub isREQ {
726705
}
727706
}
728707

708+
=head3 C<NGS>
709+
710+
Normalized Gram Schmidt basis for a given list of vectors. This uses "fuzzy" equivalency to avoid machine rounding
711+
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
712+
method with the option C<< cols => 1 >>, then the columns are used as the list of vectors.
713+
714+
In array context, returns an array of 1D Matrix objects. In scalar context, returns a Matrix with those rows.
715+
716+
If C<< cols => 1 >> was used, return the Matrix objects in column form.
717+
718+
719+
Usage:
720+
721+
$A = Matrix([1, 2, 2], [2, 2, 1]);
722+
723+
@X = $A->NGS # @X is (Matrix(1, 2, 2)/3, Matrix(10, 2, -7)/sqrt(153))
724+
@X = Value::Matrix->NGS([1, 2, 2], [2, 2, 1]) # @X is the same as above
725+
$X = $A->NGS # $X is Matrix([1/3, 2/3, 2/3], [10/sqrt(153), 2/sqrt(153), -7/sqrt(153)])
726+
$X = Value::Matrix->NGS([1, 2, 2], [2, 2, 1]) # $X is the same as above
727+
728+
@X = $A->NGS(cols => 1) # @X is (Matrix([1/sqrt(5)], [2/sqrt(5)]), Matrix([2/sqrt(5)], [-1/sqrt(5)]))
729+
$X = $A->NGS(cols => 1) # $X is Matrix([1/sqrt(5), 2/sqrt(5)], [2/sqrt(5), -1/sqrt(5)])
730+
731+
=cut
732+
733+
sub NGS {
734+
my ($self, @args) = @_;
735+
my %options;
736+
my @rows;
737+
if (ref($self) eq 'Value::Matrix') {
738+
%options = @args;
739+
@rows = $options{cols} ? $self->transpose->value : $self->value;
740+
} else {
741+
@rows = @args;
742+
}
743+
Value::Error('You must provide vectors to apply Gram Schmidt to') if !@rows;
744+
@rows = map { Value::Matrix->new($_) } @rows;
745+
@rows = grep { !$_->isZero } @rows;
746+
Value::Error('You must provide at least one nonzero row for Gram Schmidt') unless @rows;
747+
my $n = ($rows[0]->dimensions)[0];
748+
for my $r (@rows) {
749+
Value::Error('Rows provided for Gram Schmidt should not be nested arrays') unless $r->degree == 1;
750+
Value::Error('All rows provided for Gram Schmidt should have the same length') unless ($r->dimensions)[0] == $n;
751+
}
752+
@rows = map { $_ / ($_ * $_)**0.5 } @rows;
753+
my @gs = ($rows[0]);
754+
for my $i (1 .. $#rows) {
755+
my $proj = Value::Matrix->new((0) x $n);
756+
for my $v (@gs) {
757+
$proj += ($rows[$i] * $v) * $v;
758+
}
759+
my $gs = $rows[$i] - $proj;
760+
push @gs, $gs / ($gs * $gs)**0.5 unless $rows[$i] == $proj;
761+
}
762+
return
763+
$options{cols}
764+
? wantarray
765+
? map { $_->transpose } @gs
766+
: Value::Matrix->new(@gs)->transpose
767+
: wantarray ? @gs
768+
: Value::Matrix->new(@gs);
769+
}
770+
729771
sub _isNumber {
730772
my $n = shift;
731773
return Value::isNumber($n) || Value::classMatch($n, 'Fraction');

‎t/math_objects/matrix.t‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,44 @@ subtest 'Check if two matrices are (fuzzy) row equivalent' => sub {
224224

225225
};
226226

227+
subtest 'Normalized Gram Schmidt' => sub {
228+
my @A = ([ 1, 2, 2 ], [ 2, 2, 1 ]);
229+
my $A = Matrix(@A);
230+
my @NGSA = $A->NGS;
231+
my @NGSAc = $A->NGS(cols => 1);
232+
my @NGSVMA = Value::Matrix->NGS(@A);
233+
my $NGSA = $A->NGS;
234+
my $NGSAc = $A->NGS(cols => 1);
235+
my $NGSVMA = Value::Matrix->NGS(@A);
236+
237+
is "@NGSA", "[0.333333,0.666667,0.666667] [0.808452,0.16169,-0.565916]",
238+
'Test that NGS finds a correct normalized GS basis';
239+
is "@NGSAc", "[[0.447214],[0.894427]] [[0.894427],[-0.447214]]",
240+
'Test that NGS finds a correct normalized GS basis';
241+
is "@NGSVMA", "[0.333333,0.666667,0.666667] [0.808452,0.16169,-0.565916]",
242+
'Test that NGS finds a correct normalized GS basis';
243+
is "$NGSA", "[[0.333333,0.666667,0.666667],[0.808452,0.16169,-0.565916]]",
244+
'Test that NGS finds a correct normalized GS basis';
245+
is "$NGSAc", "[[0.447214,0.894427],[0.894427,-0.447214]]", 'Test that NGS finds a correct normalized GS basis';
246+
is "$NGSVMA", "[[0.333333,0.666667,0.666667],[0.808452,0.16169,-0.565916]]",
247+
'Test that NGS finds a correct normalized GS basis';
248+
249+
like dies {
250+
Value::Matrix->NGS();
251+
}, qr/You must provide vectors to apply Gram Schmidt to/, 'Test that you must pass something as an argument';
252+
like dies {
253+
Value::Matrix->NGS([ 0, 0 ], [ 0, 10**-16 ]);
254+
}, qr/You must provide at least one nonzero row for Gram Schmidt/,
255+
'Test that you must pass at least one nonzero row';
256+
like dies {
257+
Value::Matrix->NGS([ [ 1, 0 ], [ 1, 0 ] ]);
258+
}, qr/Rows provided for Gram Schmidt should not be nested arrays/,
259+
'Test that rows are rows, not nested matrices';
260+
like dies {
261+
Value::Matrix->NGS([ 1, 0 ], [ 1, 0, 0 ]);
262+
}, qr/All rows provided for Gram Schmidt should have the same length/, 'Test that rows have same length';
263+
};
264+
227265
subtest 'Transpose a Matrix' => sub {
228266
my $A = Matrix([ [ 1, 2, 3, 4 ], [ 5, 6, 7, 8 ], [ 9, 10, 11, 12 ] ]);
229267
my $B = Matrix([ [ 1, 5, 9 ], [ 2, 6, 10 ], [ 3, 7, 11 ], [ 4, 8, 12 ] ]);

0 commit comments

Comments
 (0)