Outer product of two rectangular matrices

341 Ansichten (letzte 30 Tage)
Dario Cortese
Dario Cortese am 19 Feb. 2019
Beantwortet: Branco Juan am 7 Dez. 2021
If I have a vector r , I can easily calculate its inner product
r=[1 2 3];
inner = r*r'
inner =
14
Same goes for the outer product (please see here for complete definition)
outer=r'*r
outer =
1 2 3
2 4 6
3 6 9
outer, as it should be, has components (where N is the total number of components, here 3). inner, on the other hand has components (where m is the number of rows, here 1).
I want to be able to do this standard thing to rectangular matrices too. The inner product of rectangular matrices is easy enough:
r =
1 2 3
1 1 1
inner=r*r'
inner =
14 6
6 3
inner has components (2x2=4) and this is what I expect from the matrix multiplication of r with its transponse.
Clearly, though, it is not clear how I should calculate the outer product of r with itself, because now the definition of "inner product with transponse" and "outer product with itself" have the same syntax in matlab. Indeed, if I try to repeat what I have done for the vector r, I obtain:
outer=r'*r
outer =
2 3 4
3 5 7
4 7 10
which is not the outer product of r with itself, as it's evident from the fact that it does not have =36, but only components (where n is the number of column). What matlab has interpreted my calculation to be is the inner product of r transponse and r.
How do I obtain the proper outer product, whose components are all the combinations of products between components of r?

Antworten (3)

Jan
Jan am 19 Feb. 2019
Bearbeitet: Jan am 19 Feb. 2019
If you have a [N x M] array and a [R x S x T] array, the output product becomes [N x M x R x S x T]. This can be done with nested for loops:
function C = OuterProduct(A, B) % version 1
sizeA = size(A);
sizeB = size(B);
C = zeros([sizeA, sizeB]);
if ndims(A) == 2 && ndims(B) == 2
for b2 = 1:sizeB(2)
for b1 = 1:sizeB(1)
for a2 = 1:sizeA(2)
for a1 = 1:sizeA(1)
C(a1, a2, b1, b2) = A(a1, a2) * B(b1, b2);
end
end
end
end
else
error('Not implemented yet.');
end
end
Is this correct so far? This is a primitive approach and it can be improved with respect to run-time. But this is a proof of concept yet. Now an improvement:
function C = OuterProduct(A, B) % version 2
sizeA = size(A);
sizeB = size(B);
C = zeros([sizeA, sizeB]);
if ndims(A) == 2 && ndims(B) == 2
for b2 = 1:sizeB(2)
for b1 = 1:sizeB(1)
%for a2 = 1:sizeA(2)
% for a1 = 1:sizeA(1)
C(:, :, b1, b2) = A * B(b1, b2);
% end
%end
end
end
else
error('Not implemented yet.');
end
end
A general solution is still ugly:
function C = OuterProduct(A, B) % version 3
C = zeros([size(A), size(B)]);
q = ndims(A) + 1;
index = cell(1, q);
index(1:q-1) = {':'}; % A cell as index vector of dynamic length
for k = 1:numel(B)
index{q} = k; % Linear index for elements of B
C(index{:}) = A * B(k);
end
end
Nicer and faster with bsxfun or the modern automatic expanding and the elementwise product:
function C = OuterProduct(A, B) % version 4
C = A .* reshape(B, [ones(1, ndims(A)), size(B)]);
% Matlab < R2016b:
% C = bsxfun(@times, A, reshape(B, [ones(1, ndims(A)), size(B)]))
end
Note: There is no reason to reshape A also:
AA = reshape(A, [size(A), ones(1, ndims(B))])
size(A)
size(AA) % Has still the same size!
Trailing dimensions of length 1 are ignored in Matlab.
Alternatively (see Matt J's solution):
function C = OuterProduct(A, B) % version 5
C = reshape(A(:) * B(:).', [size(A), size(B)]);
end
I love Matlab for the last two versions.

Matt J
Matt J am 19 Feb. 2019
Bearbeitet: Matt J am 19 Feb. 2019
outer = r(:)*r(:).'

Branco Juan
Branco Juan am 7 Dez. 2021
Notation: majuscule for Tensors, Matrices and Vectors; minuscule for their elements
In MatLab, the operator * is always the Matrix Product of matrices (tensors), which means it is the Dot Product for Vectors in Euclidean space (the Inner Product <V1;V2>=V1'.M.V2 with M being the identity).
It requires matrix dimensions to agree and suffices all due properties of such (preserving the order of the two tensors being operated, combining the dimensions involved and changing the number of elements)
Being for example V1 a raw vector of order 1 and dimension dim_v, that can be seen as a tensor of order 2 and dimensions (1;dim_v); V2 a column vector of order 1 and dimension dim_v, that can be seen as a tensor of order 2 and dimensions (dim_v;1)
DP = V1*V2 = sum( v1(i)*v2(i) ) for all i=1:dim_v, of order 2 (although it seems like order 1)
So that DP is a tensor (1x1) and the element dp(1,1) = v1(1,1)*v2(1,1) + v1(1,2)*v2(2,1) + ...v1(1,dim_v)*v2(dim_v,1)
So far, standard stuff.
Alternatively, being A a tensor of order 2 and dimensions (dim_1,dim_2); B a tensor of order 2 and dimensions (dim_2,dim_3), more general than your r and r' in the question above, where dim_1 = dim_3.
MP = A*B, of order 2
So that MP is a tensor (dim_1 x dim_3) and the elements mp(i,j) = sum( a(i,k)*b(k,j) ) for all k=1:dim_2
Do not confuse the operator * with the concept of Outer Product (well defined in the reference you've provided)
OP = A X B, of order 4
So that OP is a tensor (dim_1 x dim_2 x dim_2 x dim_3) and the elements op(i,j,k,g) = a(i,j)*b(k,g)
Matlab practitioners use the operator * to compute the Outer Product of Vectors because by chance it gives the same result, as you can see here:
MPV = V2*V1, of order 2
So that MPV is a tensor (dim_v x dim_v) and the elements mpv(i,j) = sum( a(i,k)*b(k,j) ) for all k=1:1;
thus mpv(i,j) = v1(i,1)*v2(1,j)
OPV = V1 X V2, of order 4
So that OPV is a tensor (1 x dim_v x dim_v x 1) and the elements opv(1,i,j,1) = v1(1,i)*v2(j,1) et voilà!!
You'll see the result as a (dim_v x dim_v) Matrix because we tend to ignore the dimensions of 1. But pay attention to the possible meaning of such possitions with respect to the final application.
So, once all the maths behind your problem are clear (I hope so), let's compute Outer Product of Tensors:
OPM = B X A
So that OPM should be a tensor of order 4, (dim_2 x dim_3 x dim_1 x dim_2) and the elements opm(k,g,i,j) = b(k,g)*a(i,j)
Solved via the new Matrix G of order 4 with the elements of B in its first two dimensions and repeated in the next two, and a new Matrix H of order 4 with the elements of A but in its first two dimensions, repeated in the next two, and finally permuted properly.
G should be (dim_2 x dim_3 x dim_1 x dim_2)
H should be (dim_2 x dim_3 x dim_1 x dim_2)
OPM = G.*H the element-wise product, or element-by-element multiplication.
So that OPM should be a tensor (dim_2 x dim_3 x dim_1 x dim_2) and the elements omp(m,n,i,j) = g(m,n)*h(i,j)
To solve your precise example:
A=[1,2,3;1,1,1]; % This is your r, a matrix (2x3)
B=A'; % This is your r', a matrix (3x2)
% you want to solve (r') Outer_Product (r), so Result = B Outer_Product A a
% Matrix (3x2x2x3) with 36 elements in total
% let's create the new matrices
G = repmat(B,[1 1 size(A)]); % this is a matrix (3x2x2x3)
Aux = repmat(A,[1 1 size(B)]); % this is a matrix (2x3x3x2), so we must reposition its elements
H = permute(Aux,[3 4 1 2]); % this is a matrix (3x2x2x3)
% the order [3 4 1 2] is applicable to any size matrices if they are of
% order 2.
Result = G.*H; % That's it
% Faster, without explanation
Result = B.*(permute(repmat(A,[1 1 size(B)]),[3 4 1 2]));
% No need to reshape B because MatLab is smart.
% Elapsed time is min 0.000481 seconds and max 0.001386 seconds
% reshape(B(:) * A(:).', [size(B), size(A)]); takes min 0.000244 seconds
% and max 0.000721 seconds
% B .* reshape(A, [ones(1, ndims(B)), size(A)]); takes min 0.000635 seconds
% and max 0.001464 seconds
% reshape(B, [size(B), ones(1, ndims(A))]); gives a wrong result...

Kategorien

Mehr zu Matrices and Arrays finden Sie in Help Center und File Exchange

Produkte


Version

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by