######################################################################## ## ## Copyright (C) 2022 Kasper H. Filtenborg ## ## This program is free software: you can redistribute it and/or modify ## it under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## This program is distributed in the hope that it will be useful, ## but WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with this program. If not, see . ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{C} =} tensorprod (@var{A}, @var{B}, @var{dimA}, @var{dimB}) ## @deftypefnx {} {@var{C} =} tensorprod (@var{A}, @var{B}, @var{dim}) ## @deftypefnx {} {@var{C} =} tensorprod (@var{A}, @var{B}) ## @deftypefnx {} {@var{C} =} tensorprod (@var{A}, @var{B}, "all") ## @deftypefnx {} {@var{C} =} tensorprod (@var{A}, @var{B}, @dots{}, "NumDimensionsA", @var{value}) ## Compute the tensor product between tensors @var{A} and @var{B}. ## ## The dimensions of @var{A} and @var{B} that are contracted are defined by ## @var{dimA} and @var{dimB}, respectively. @var{dimA} and @var{dimB} are ## vectors with equal length and define the dimensions to match up. The ## matched dimensions of @var{A} and @var{B} must have equal length. ## ## When @var{dim} is used, it is equivalent to @var{dimA} = @var{dimB} = ## @var{dim}. ## ## Running without an additional argument is equivalent to @var{dimA} = ## @var{dimB} = []. This computes the outer product between @var{A} and ## @var{B}. ## ## Using the "all" option results in the inner product between @var{A} and ## @var{B}. For this, it is required that size(@var{A}) == size(@var{B}). ## ## Use the property-value pair with the property name "NumDimensionsA" ## when @var{A} has trailing singleton dimensions that should be transfered to ## @var{C}. The specified value should be the number of dimensions of @var{A}. ## ## @seealso{kron, dot, mtimes} ## @end deftypefn function C = tensorprod (A, B, varargin) if (nargin == 0) print_usage (); elseif (nargin < 2) error ("tensorprod: too few inputs given"); elseif (nargin > 6) error ("tensorprod: too many inputs given"); endif ## Check that A and B are single or double if (! isfloat (A)) error ("tensorprod: A must be a single or double precision array"); endif if (! isfloat (B)) error ("tensorprod: B must be a single or double precision array"); endif ## Check for misplaced NumDimensionsA property if (nargin > 2) if (strcmpi (varargin{end}, "NumDimensionsA")) error (["tensorprod: a value for the NumDimensionsA property must ", ... "be provided"]); elseif (strcmpi ( strtok (inputname (nargin, false)), "NumDimensionsA")) # FIXME: Add support for keyword=value syntax error (["tensorprod: NumDimensionsA=ndimsA syntax is not yet ", ... "supported in Octave - provide the value as a property-value pair"]); endif endif ## Check for NumDimensionsA property if (nargin > 3) if (strcmpi (varargin{end - 1}, "NumDimensionsA")) if (! (isnumeric (varargin{end}) && isscalar (varargin{end}))) error (["tensorprod: value for NumDimensionsA must be a ", ... "numeric scalar"]); elseif (varargin{end} < 1 || mod (varargin{end}, 1) != 0) error (["tensorprod: value for NumDimensionsA must be a ", ... "positive integer"]); endif NumDimensionsA = varargin{end}; endif endif existNumDimensionsA = exist ("NumDimensionsA"); ndimargs = nargin - 2 - 2 * existNumDimensionsA; ## Set dimA and dimB if (ndimargs == 0) ## Calling without dimension arguments dimA = []; dimB = []; elseif (ndimargs == 1) ## Calling with dim or "all" option if (isnumeric (varargin{1})) if (! (isvector (varargin{1}) || isnull (varargin{1}))) error ("tensorprod: dim must be a numeric vector of integers or []"); endif ## Calling with dim dimA = transpose ([varargin{1}(:)]); elseif (ischar (varargin{1})) if (strcmpi (varargin{1}, "all")) if (! size_equal (A, B)) error (["tensorprod: size of A and B must be identical when ", ... "using the \"all\" option"]); endif else error ("tensorprod: unknown option \"%s\"", varargin{1}); endif ## Calling with "all" option dimA = 1:ndims(A); else error (["tensorprod: third argument must be a numeric vector of ", ... "integers, [], or \"all\""]); endif dimB = dimA; elseif (ndimargs == 2) ## Calling with dimA and dimB if (! (isnumeric (varargin{1}) && (isvector (varargin{1}) || ... isnull (varargin{1})))) error("tensorprod: dimA must be a numeric vector of integers or []"); endif if (! (isnumeric (varargin{2}) && (isvector (varargin{2}) || ... isnull (varargin{2})))) error ("tensorprod: dimB must be a numeric vector of integers or []"); endif if (length (varargin{1}) != length (varargin{2})) error (["tensorprod: an equal number of dimensions must be ", ... "matched for A and B"]); endif dimA = transpose ([varargin{1}(:)]); dimB = transpose ([varargin{2}(:)]); else ## Something is wrong - try to find the error for i = 1:ndimargs if (ischar (varargin{i})) if (strcmpi (varargin{i}, "NumDimensionsA")) error ("tensorprod: misplaced \"NumDimensionsA\" option"); elseif (strcmpi (varargin{i}, "all")) error ("tensorprod: misplaced \"all\" option"); else error ("tensorprod: unknown option \"%s\"", varargin{i}); endif elseif (! isnumeric (varargin{i})) error (["tensorprod: optional arguments must be numeric vectors ", ... "of integers, [], \"all\", or \"NumDimensionsA\""]); endif endfor error ("tensorprod: too many dimension inputs given"); endif ## Check that dimensions are positive integers ([] will also pass) if (any ([dimA < 1, dimB < 1, (mod (dimA, 1) != 0), (mod (dimB, 1) != 0)])) error ("tensorprod: dimension(s) must be positive integer(s)"); endif ## Check that the length of matched dimensions are equal if (any (size (A, dimA) != size (B, dimB))) error (["tensorprod: matched dimension(s) of A and B must have the ", ... "same length(s)"]); endif ## Find size and ndims of A and B ndimsA = max ([ndims(A), max(dimA)]); sizeA = size (A, 1:ndimsA); ndimsB = max ([ndims(B), max(dimB)]); sizeB = size (B, 1:ndimsB); ## Take NumDimensionsA property into account if (existNumDimensionsA) if (NumDimensionsA < ndimsA) if (ndimargs == 1) error (["tensorprod: highest dimension of dim must be less than ", ... "or equal to NumDimensionsA"]); elseif (ndimargs == 2) error (["tensorprod: highest dimension of dimA must be less ", ... "than or equal to NumDimensionsA"]); else error (["tensorprod: NumDimensionsA cannot be smaller than the ", ... "number of dimensions of A"]); endif elseif (NumDimensionsA > ndimsA) sizeA = [sizeA, ones(1, NumDimensionsA - ndimsA)]; ndimsA = NumDimensionsA; endif endif ## Interchange the dimension to sum over the end of A and the front of B ## Prepare for A remainDimA = setdiff (1:ndimsA, dimA); # Dimensions of A to keep newDimOrderA = [remainDimA, dimA]; # New order of dimensions (dimensions to keep first, dimensions to contract last) newSizeA = [prod(sizeA(remainDimA)), prod(sizeA(dimA))]; # Size of temporary 2D representation of A remainSizeA = sizeA(remainDimA); # Contribution to size of C from remaining dimensions of A ## Prepare for B remainDimB = setdiff (1:ndimsB, dimB); # See comments for A newDimOrderB = [remainDimB, dimB]; newSizeB = [prod(sizeB(remainDimB)), prod(sizeB(dimB))]; # In principle, prod(sizeB(dimB)) should always be equal to prod(sizeA(dimA)) remainSizeB = sizeB(remainDimB); ## Do reshaping into 2D array newA = reshape (permute (A, newDimOrderA), newSizeA); newB = reshape (permute (B, newDimOrderB), newSizeB); ## Compute C = newA * transpose (newB); ## If not an inner product, reshape back to tensor if (! isscalar (C)) C = reshape (C, [remainSizeA, remainSizeB]); endif endfunction %!test %! rand ("seed", 0); %! A = rand (3, 2, 5); %! B = rand (2, 4, 5); %! v = rand (4, 1); %! assert ( tensorprod (A, B)); %! assert ( tensorprod (A, A, 1)); %! assert ( tensorprod (A, A, 4)); %! assert ( tensorprod (A, A, [])); %! assert ( tensorprod (A, A, [2, 3])); %! assert ( tensorprod (A, B, 2, 1)); %! assert ( tensorprod (A, B, 4, 4)); %! assert ( tensorprod (A, B, [2, 3], [1, 3])); %! assert ( tensorprod (A, A, [], [])); %! assert ( A(1), 0.9999996423721313330669073875470); %! assert ( tensorprod (A, A, "all"), 11.258348100536915329070518200751); %! assert ( tensorprod (A, A, "all")); %! assert ( tensorprod (A, A, 1, "NumDimensionsA", 4)); %! assert ( tensorprod (A, A, 4, "NumDimensionsA", 4)); %! assert ( tensorprod (A, A, 4, "numdimensionsa", 4)); %! assert ( tensorprod (A, A, [2, 3], "NumDimensionsA", 4)); %! assert ( tensorprod (A, A, [], "NumDimensionsA", 4)); %! assert ( tensorprod (A, B, 2, 1, "NumDimensionsA", 4)); %! assert ( tensorprod (A, B, [2, 3], [1, 3], "NumDimensionsA", 4)); %! assert ( tensorprod (A, B, [2, 3], [1; 3], "NumDimensionsA", 4)); %! assert ( tensorprod (A, B, [], [], "NumDimensionsA", 4)); %! assert ( tensorprod (1, 2), 2); %! assert ( tensorprod (v, v, "all"), dot (v, v)); %! assert ( tensorprod (v, v), reshape (v * transpose (v), [4, 1, 4])); ## Test empty inputs %!assert ( tensorprod ([], []), zeros (0, 0, 0, 0)) %!assert ( tensorprod ([], 1), []) %!assert ( tensorprod (1, []), zeros (1, 1, 0, 0)) %!assert ( tensorprod (zeros (0, 0, 0), zeros (0, 0, 0)), zeros (0, 0, 0, 0, 0, 0)) %!assert ( tensorprod ([], [], []), zeros (0, 0, 0, 0)) %!assert ( tensorprod ([], [], 1), []) %!assert ( tensorprod ([], [], 2), []) %!assert ( tensorprod ([], [], 3), zeros (0, 0, 0, 0)) %!assert ( tensorprod ([], [], 4), zeros (0, 0, 1, 0, 0)) %!assert ( tensorprod ([], [], 5), zeros (0, 0, 1, 1, 0, 0)) %!assert ( tensorprod ([], [], 3, "NumDimensionsA", 4), zeros (0, 0, 1, 0, 0)) %!assert ( tensorprod ([], [], 3, 4, "NumDimensionsA", 5), zeros (0, 0, 1, 1, 0, 0)) ## Test input validation %!error tensorprod () %!error tensorprod (1) %!error tensorprod ("foo", 1) %!error tensorprod (1, "bar") %!error tensorprod (1, 1, "foo") %!error tensorprod (1, 1, "foo", 1) %!error tensorprod (1, 1, 1, "bar") %!error tensorprod (1, 1, 1, "foo", 1) %!error tensorprod (1, 1, 1, "all", 1) %!error tensorprod (1, 1, "NumDimensionsA", 1, 1) %!error tensorprod (1, 1, 1, {}, 1) %!error tensorprod (ones (3, 4), ones (4, 3), 1) %!error tensorprod (1, 1, 0) %!error tensorprod (1, 1, -1) %!error tensorprod (1, 1, 1.5) %!error tensorprod (1, 1, NaN) %!error tensorprod (1, 1, Inf) %!error tensorprod (1, 1, {}) %!error tensorprod (1, 1, zeros(0,0,0)) %!error tensorprod (1, 1, zeros(0,0,0), []) %!error tensorprod (1, 1, [], zeros(0,0,0)) %!error tensorprod (ones (3, 4), ones (4, 3), 1, 1) %!error tensorprod (ones (3, 4), ones (4, 3), 1, [1, 2]) %!error tensorprod (ones (3, 4), ones (4, 3), "all") %!error tensorprod (1, 1, "NumDimensionsA") %!error tensorprod (ones (2, 2, 2), 1, "NumDimensionsA", 2) %!error tensorprod (1, 1, 5, "NumDimensionsA", 4) %!error tensorprod (1, 1, 5, 5, "NumDimensionsA", 4) %!error tensorprod (1, 1, NumDimensionsA=4) %!error tensorprod (1, 1, numdimensionsa=4) %!error tensorprod (1, 1, 2, 1, 1) %!error tensorprod (1, 1, 2, 1, 1, 1) %!error tensorprod (1, 1, 2, 1, 1, 1, 1) %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", "foo") %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", {}) %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", -1) %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", 0) %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", 1.5) %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", NaN) %!error tensorprod (1, 1, 2, 1, "NumDimensionsA", Inf)