3 に答える
In a multi-dimensional (tensor) case of y = u⊗v, I believe that you need to shift the dimensions of the second operand like so:
v_t = permute(v, circshift(1:(ndims(u) + ndims(v)), [0, ndims(u)]));
and then multiply them with bsxfun:
y = bsxfun(@times, u, v_t);
The regular matrix multiplication is defined only for vector and 2-D matrices, so we couldn't use it in the general case.
Also note that this computation still fails if the second operand is a 1-D vector, because ndims returns 2 instead of 1 for vectors. For this purpose, lets define our own function that counts dimensions:
my_ndims = @(x)(isvector(x) + ~isvector(x) * ndims(x));
To complete the answer, you can define a new function (e.g. an anonymous function), like so:
outprod = @(u, v)bsxfun(@times, u, permute(v, circshift(1:(my_ndims(u) + my_ndims(v)), [0, my_ndims(u)])));
and then use it as many times as you want. For example, y = a×a×a would be computed like so:
y = outprod(outprod(a, a), a);
Of course, you can write a better function that takes a variable number of arguments to save you some typing. Something along these lines:
function y = outprod(u, varargin)
my_ndims = @(x)(isvector(x) + ~isvector(x) * ndims(x));
y = u;
for k = 1:numel(varargin)
v = varargin{k};
v_t = permute(v, circshift(1:(my_ndims(y) + my_ndims(v)),[0, my_ndims(y)]));
y = bsxfun(@times, y, v_t);
end
I hope I got the math right!
You can use as well the kron function:
kron(a * a', a)
or when four outer (kronecker tensor) products needed:
kron(kron(a * a', a), a)
and so on. The last one gives you a m x n matrix, where m = n * n * n.
If adding dimensions is desired as going on with the products, you may use the reshape function:
reshape(kron(a * a', a), [n, n, n])
or
reshape(kron(kron(a * a', a), a), [n, n, n, n])
and so on. The last one gives you a n x n x n x n tensor.
The problem with using kron as in a previous solution is that it throws off canonical indexing of the outerproduct.
Instead, ndgrid is ideal for this scenario:
a = [1; 2; 3];
b = [4; 5];
c = [6; 7; 8; 9];
[xx, yy, zz] = ndgrid(1:length(a), 1:length(b), 1:length(c));
% desired outerproduct
M = a(xx) .* b(yy) .* c(zz);
On paper, we can check that the desired solution M is the datacube:
M(:,:,1) = | M(:,:,2) = | M(:,:,3) = | M(:,:,4) =
| | |
24 30 | 28 35 | 32 40 | 36 45
48 60 | 56 70 | 64 80 | 72 90
72 90 | 84 105 | 96 120 | 108 135
Using the Kronecker product approach
M2 = reshape(kron(a * b', c), [length(a), length(b), length(c)]);
we would get:
M2(:,:,1) = | M2(:,:,2) = | M2(:,:,3) = | M2(:,:,4) =
| | |
24 36 | 64 84 | 30 45 | 80 105
28 48 | 72 96 | 35 60 | 90 120
32 56 | 72 108 | 40 70 | 90 135
Datacube M2 has the same elements as M, but these elements are rearranged. This is because kron(a * b', c) does not contain the slices of M in contiguous blocks to facilitate direct application of the reshape function. To compute outerproduct this way, we would need to apply a rearrangement operation/function (determinable, but labrious and time consuming) to the elements of kron(a * b', c).
A further advantage of using ndgrid is that it generalizes easily to higher orders.