Skip to content

tensor multiplication, powers, identity, inverses #1215

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pstaabp merged 6 commits into from
Apr 1, 2025
Merged

tensor multiplication, powers, identity, inverses #1215

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
a7e2649
tensor multiplication, powers, idenity, inverses
Alex-Jordan Mar 23, 2025
2612e10
Fix typo in Matrix doc and add a test.
pstaabp Mar 27, 2025
098a131
Merge pull request #17 from pstaabp/matrix-add-unit-test
Alex-Jordan Mar 28, 2025
f547673
cleanup from PR#1215 review
Alex-Jordan Mar 28, 2025
b97d8dd
correct POD for ->value
Alex-Jordan Mar 29, 2025
d7f6260
some more tests for matrices
Alex-Jordan Mar 29, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
242 changes: 176 additions & 66 deletions lib/Value/Matrix.pm
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,17 +64,17 @@ Allows complex entries

=back

=head2 Creation of Matrices
=head2 Creation of Matrix MathObjects

Examples:

$M1 = Matrix(1, 2, 3); # 1D (row vector)
$M1 = Matrix(1, 2, 3); # degree 1 (row vector)
$M1 = Matrix([1, 2, 3]);

$M2 = Matrix([1, 2], [3, 4]); # 2D (matrix)
$M2 = Matrix([1, 2], [3, 4]); # degree 2 (matrix)
$M2 = Matrix([[1, 2], [3, 4]]);

$M3 = Matrix([[1, 2], [3, 4]], [[5, 6], [7, 8]]); # 3D (tensor)
$M3 = Matrix([[1, 2], [3, 4]], [[5, 6], [7, 8]]); # degree 3 (tensor)
$M3 = Matrix([[[1, 2], [3, 4]], [[5, 6], [7, 8]]]);

$M4 = ...
Expand All @@ -83,21 +83,24 @@ Examples:

=head3 Conversion

$matrix->value produces an array of numbers (for a 1D tensor) or array refs representing the rows.
These are perl numbers and array refs, not MathObjects.
$matrix->value produces an array of references to nested arrays, except at
the deepest level, where there will be the more basic MathObjects that make
up the Matrix (e.g. Real, Complex, Fraction, a mix of these, etc)

$M1->value is (1, 2, 3)
$M2->value is ([1, 2], [3, 4])
$M3->value is ([[1, 2], [3, 4]], [[5, 6], [7, 8]])

$matrix->wwMatrix produces CPAN MatrixReal1 matrix, used for computation subroutines and can only
be used on 1D tensors (row vector) or 2D tensors (matrix).
be used on a degree 1 or degree 2 Matrix.

=head3 Information

$matrix->dimensions produces an array. For an n-dimensional tensor, the array has n entries for
$matrix->dimensions produces an array. For an n-degree Matrix, the array has n entries for
the n dimensions.

$matrix->degree returns the degree of the Matrix (the depth of the nested array).

=head3 Access values

row(i) : MathObjectMatrix
Expand Down Expand Up @@ -305,7 +308,7 @@ returns C<(2,2,2)>

#
# Recursively get the dimensions of the matrix.
# Returns (d_1,d_2,...,d_n) for an nD matrix.
# Returns (d_1,d_2,...,d_n) for a degree n Matrix.
#
sub dimensions {
my $self = shift;
Expand All @@ -315,6 +318,13 @@ sub dimensions {
return ($r, $v->[0]->dimensions);
}

sub degree {
my $self = shift;
my $v = $self->data;
return 1 unless Value::classMatch($v->[0], 'Matrix');
return 1 + $v->[0]->degree;
}

#
# Return the proper type for the matrix
#
Expand All @@ -327,24 +337,31 @@ sub typeRef {

=head3 C<isSquare>

Return true is the matrix is 1D of length 1 or 2D and square, false otherwise
Return true if the Matrix is degree 1 of length 1 or its last two dimensions are equal, false otherwise

Usage:

$A = Matrix([ [ 1, 2, 3, 4 ], [ 5, 6, 7, 8 ], [ 9, 10, 11, 12 ] ]);
$B = Matrix([ [ 1, 0, 0 ], [ 0, 1, 0 ], [ 0, 0, 1 ] ]);
$A = Matrix([ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]);
$B = Matrix([ [ 1, 0, 0 ], [ 0, 1, 0 ] ]);
$C = Matrix(1);
$D = Matrix(1, 2);
$E = Matrix([ [ [ 1, 2 ], [ 3, 4 ] ] ]);
$F = Matrix([ [ [ 1, 2 ] ], [ [ 3, 4 ] ] ]);

$A->isSquare; # is '' (false)
$B->isSquare; # is 1 (true);
$A->isSquare; # is 1 (true) because it is a 3x3 matrix
$B->isSquare; # is '' (false) because it is a 2x3 matrix
$C->isSquare; # is 1 (true) because it is a degree 1 matrix of length 1
$D->isSquare; # is '' (false) because it is a degree 1 matrix of length 2
$E->isSquare; # is 1 (true) because it is a 1x2x2 tensor
$F->isSquare; # is '' (false) because it is a 2x1x2 tensor

=cut

sub isSquare {
my $self = shift;
my @d = $self->dimensions;
return 1 if scalar(@d) == 1 && $d[0] == 1;
return 0 if scalar(@d) != 2;
return $d[0] == $d[1];
return $d[0] == 1 if scalar(@d) == 1;
return $d[-1] == $d[-2];
}

=head3 C<isRow>
Expand All @@ -369,7 +386,7 @@ sub isRow {

=head3 C<isOne>

Check for identity matrix.
Check for identity Matrix (for degree > 2, this means the last two dimensions hold identity matrices)

Usage:

Expand All @@ -379,19 +396,33 @@ Usage:
$B = Matrix([ [ 1, 0, 0 ], [ 0, 1, 0 ], [ 0, 0, 1 ] ]);
$B->isOne; # is true;

$C = Matrix(1);
$C->isOne; # is true

$D = Matrix([ [ [ 1, 0 ], [ 0, 1 ] ], [ [ 1, 0 ], [ 0, 1 ] ] ]);
$D->isOne; # is true

=cut

sub isOne {
my $self = shift;
return 0 unless $self->isSquare;
my $i = 0;
for my $row (@{ $self->{data} }) {
my $j = 0;
for my $k (@{ $row->{data} }) {
return 0 unless $k eq (($i == $j) ? "1" : "0");
$j++;
if ($self->degree <= 2) {
my $i = 0;
for my $row (@{ $self->{data} }) {
my $j = 0;
for my $k (@{ $row->{data} }) {
return 0 unless $k eq (($i == $j) ? "1" : "0");
$j++;
}
$i++;
}
} else {
for my $row (@{ $self->{data} }) {
if (!$row->isOne) {
return 0;
}
}
$i++;
}
return 1;
}
Expand Down Expand Up @@ -597,36 +628,73 @@ sub mult {
return $self->make(@coords);
}

# Make points and vectors into columns if they are on the right.
$r = !$flag && Value::classMatch($r, 'Point', 'Vector') ? ($self->promote($r))->transpose : $self->promote($r);
# Special case if $r is a Point or Vector and if $l is a degree 2 Matrix,
# we promote $r to a degree 1 Matrix and later restore the product to be Point or Vector
$r =
!$flag
&& Value::classMatch($r, 'Point', 'Vector')
&& $l->degree == 2 ? ($self->promote($r))->transpose : $self->promote($r);

if (!$flag && Value::classMatch($r, 'Point', 'Vector')) { $r = ($self->promote($r))->transpose }
else { $r = $self->promote($r) }
#
if ($flag) { my $tmp = $l; $l = $r; $r = $tmp }
if ($flag) { ($l, $r) = ($r, $l) }
my @dl = $l->dimensions;
my @dr = $r->dimensions;
if (scalar(@dl) == 1) { @dl = (1, @dl); $l = $self->make($l) }
if (scalar(@dr) == 1) { @dr = (@dr, 1); $r = $self->make($r)->transpose }
Value::Error("Can only multiply 2-dimensional matrices") if scalar(@dl) > 2 || scalar(@dr) > 2;
Value::Error("Matrices of dimensions %dx%d and %dx%d can't be multiplied", @dl, @dr) unless ($dl[1] == $dr[0]);
my @l = $l->value;
my @r = $r->value;

# Special case if $r and $l are both degree 1, perform dot product
if (scalar(@dl) == 1 && scalar(@dr) == 1) {
Value::Error("Can't multiply degree 1 matrices of different lengths") unless ($dl[0] == $dr[0]);
my $result = 0;
for my $i (0 .. $dl[0] - 1) {
$result += $l[$i] * $r[$i];
}
return $result;
}

# Perform matrix multiplication.
# Promote degree 1 Matrix to degree 2, as row or column depending on size
# Later restore result to degree 1 if appropriate
my $l1 = scalar(@dl) == 1;
my $r1 = scalar(@dr) == 1;
if ($l1) { @dl = (1, @dl); $l = $self->make($l); @l = $l->value }
if ($r1) { @dr = (@dr, 1); $r = $self->make($r)->transpose; @r = $r->value }

# Multiplication is only possible when dimensions are Zxmxn and Zxnxk
my $bad_dims;
if (scalar(@dl) != scalar(@dr)) {
$bad_dims = 1;
} elsif ($dl[-1] != $dr[-2]) {
$bad_dims = 1;
} else {
for my $i (0 .. scalar(@dl) - 3) {
if ($dl[$i] != $dr[$i]) {
$bad_dims = 1;
last;
}
}
}
Value::Error("Matrices of dimensions %s and %s can't be multiplied", join('x', @dl), join('x', @dr)) if $bad_dims;

my @l = $l->value;
my @r = $r->value;
# Perform matrix/tensor multiplication.
my @M = ();
for my $i (0 .. $dl[0] - 1) {
my @row = ();
for my $j (0 .. $dr[1] - 1) {
my $s = 0;
for my $k (0 .. $dl[1] - 1) { $s += $l[$i]->[$k] * $r[$k]->[$j] }
push(@row, $s);
if (scalar(@dl) == 2) {
my @row = ();
for my $j (0 .. $dr[1] - 1) {
my $s = 0;
for my $k (0 .. $dl[1] - 1) { $s += $l[$i]->[$k] * $r[$k]->[$j] }
push(@row, $s);
}
push(@M, $self->make(@row));
} else {
push(@M, $l->data->[$i] * $r->data->[$i]);
}
push(@M, $self->make(@row));
}
$self = $self->inherit($other) if Value::isValue($other);
return $self->make(@M);

if ($l1) { return $self->make(@M)->row(1) }
if ($r1) { return $self->make(@M)->transpose->row(1) }

return $other->new($self->make(@M));
}

sub div {
Expand All @@ -653,10 +721,23 @@ sub power {
$self->Error("Matrix is not invertible") unless defined($l);
}
Value::Error("Matrix powers must be non-negative integers") unless _isNumber($r) && $r =~ m/^\d+$/;
return $context->Package("Matrix")->I($l->length, $context) if $r == 0;
my $M = $l;
for my $i (2 .. $r) { $M = $M * $l }
return $M;
return $l if $r == 1;
my @powers = ($l);
my @d = $l->dimensions;
my $d = pop @d;
pop @d;
my $return = $context->Package("Matrix")->I(\@d, $d);
my $p = $r;

while ($p) {
if ($p % 2) {
$p -= 1;
$return *= $powers[-1];
}
push(@powers, $powers[-1] * $powers[-1]);
$p /= 2;
}
return $return;
}

# Do lexicographic comparison (row by row)
Expand Down Expand Up @@ -718,31 +799,51 @@ sub transpose {

=head3 C<I>

Get an identity matrix of the requested size
Get an identity Matrix of the requested size

Value::Matrix->I(n)

Usage:

Value::Matrix->I(3); # returns a 3 by 3 identity matrix.
$A->I; # return an n by n identity matrix, where n is the number of rows of A
Value::Matrix->I(3); # returns a 3x3 identity matrix.
Value::Matrix->I([2],3); # returns a 2x3x3 identity Matrix (tensor)
Value::Matrix->I([2,4],2); # returns a 2x4x2x2 identity Matrix (tensor)
$A->I; # return an identity matrix of the appropriate degree and dimensions such that I*A = A

=cut

sub I {
my $self = shift;
my $d = shift;
my $context = shift || $self->context;
$d = ($self->dimensions)[0] if !defined $d && ref($self);

Value::Error("You must provide a dimension for the Identity matrix") unless defined $d;
Value::Error("Dimension must be a positive integer") unless $d =~ m/^[1-9]\d*$/;
my $self = shift;
my @d;
my $d;
my $context;
if (ref($self) eq 'Value::Matrix') {
@d = $self->dimensions;
pop @d unless scalar(@d) == 1;
$d = pop @d;
$context = $self->context;
} else {
$d = shift; # array ref, or number
if (ref($d) eq 'ARRAY') {
@d = @{$d};
$d = shift;
}
Value::Error("You must provide a dimension for the Identity matrix") unless ($d || $d eq '0');
Value::Error("Dimension must be a positive integer") unless $d =~ m/^[1-9]\d*$/;
$context = shift || $self->context;
}

my @M = ();
my @M;
my $REAL = $context->Package('Real');

for my $i (0 .. $d - 1) {
push(@M, $self->make($context, map { $REAL->new(($_ == $i) ? 1 : 0) } 0 .. $d - 1));
if (!@d) {
for my $i (0 .. $d - 1) {
push(@M, $self->make($context, map { $REAL->new(($_ == $i) ? 1 : 0) } 0 .. $d - 1));
}
} else {
for my $i (1 .. $d[0]) {
push(@M, Value::Matrix->I([ @d[ 1 .. $#d ] ], $d));
}
}
return $self->make($context, @M);
}
Expand Down Expand Up @@ -1108,10 +1209,19 @@ sub det {

sub inverse {
my $self = shift;
$self->wwMatrixLR;
Value->Error("Can't take inverse of non-square matrix") unless $self->isSquare;
my $I = $self->{lrM}->invert_LR;
return (defined($I) ? $self->new($I) : $I);
my @d = $self->dimensions;
my @M;
for my $i (0 .. $d[0] - 1) {
if (scalar(@d) == 2) {
$self->wwMatrixLR;
Value->Error("Can't take inverse of non-square matrix") unless $self->isSquare;
my $I = $self->{lrM}->invert_LR;
return (defined($I) ? $self->new($I) : $I);
} else {
push(@M, $self->data->[$i]->inverse);
}
}
return $self->new($self->make(@M));
}

sub decompose_LR {
Expand Down
Loading