phdstudent
phdstudent

Reputation: 1144

Triple weighted sum

I was trying to vectorize a certain weighted sum but couldn't figure out how to do it. I have created a simple minimal working example below. I guess the solution involves either bsxfun or reshape and kronecker products but I still have not managed to get it working.

rng(1);
N = 200;
T1 = 5;
T2 = 7;
T3 = 10;


A = rand(N,T1,T2,T3);
w1 = rand(T1,1);
w2 = rand(T2,1);
w3 = rand(T3,1);

B = zeros(N,1);

for i = 1:N
 for j1=1:T1
  for j2=1:T2
   for j3=1:T3
    B(i) = B(i) + w1(j1) * w2(j2) * w3(j3) * A(i,j1,j2,j3);
   end
  end
 end
end

A = B;

For the two dimensional case there is a smart answer here.

Upvotes: 3

Views: 109

Answers (4)

CKT
CKT

Reputation: 781

If we're going the route of having functions anyway, and are favoring performance over elegance/brevity, then consider this:

function B = weightReduce(A, varargin)

    B = A;
    for i = length(varargin):-1:1
        N = length(varargin{i});
        B = reshape(B, [], N) * varargin{i};
    end

end

This is the performance comparison I see:

tic;
for i = 1:10000
    W = createWeights(w1,w2,w3);
    B = reshape(A, size(A,1), [])*W(:);
end
toc
Elapsed time is 0.920821 seconds.
tic; 
for i = 1:10000
    B2 = weightReduce(A, w1, w2, w3); 
end
toc
Elapsed time is 0.484470 seconds.

Upvotes: 1

Tasos Papastylianou
Tasos Papastylianou

Reputation: 22215

This is the logic behind it:

ww1 = repmat (permute (w1, [4, 1, 2, 3]), [N, 1,  T2, T3]);
ww2 = repmat (permute (w2, [3, 4, 1, 2]), [N, T1, 1,  T3]);
ww3 = repmat (permute (w3, [2, 3, 4, 1]), [N, T1, T2, 1 ]);

B = ww1 .* ww2 .* ww3 .* A;
B = sum (B(:,:), 2)

You can avoid permute by creating w1, w2, and w3 in the appropriate dimension in the first place. Also you can use bsxfun instead of repmat as appropriate for extra performance, I'm just showing the logic here and repmat is easier to follow.

EDIT: Generalised version for arbitrary input dimensions:

Dims   = {N, T1, T2, T3};  % add T4, T5, T6, etc as appropriate
Params = cell (1, length (Dims));

Params{1} = rand (Dims{:});
for n = 2 : length (Dims)
  DimSubscripts = ones (1, length (Dims));  DimSubscripts(n) = Dims{n};
  RepSubscripts = [Dims{:}];  RepSubscripts(n) = 1;
  Params{n} = repmat (rand (DimSubscripts), RepSubscripts);
end

B = times (Params{:});
B = sum (B(:,:), 2)

Upvotes: 1

Suever
Suever

Reputation: 65430

You can use an additional multiplication to modify the w1 * w2' grid from the previous answer to then multiply by w3 as well. You can then use matrix multiplication again to multiply with a "flattened" version of A.

W = reshape(w1 * w2.', [], 1) * w3.';
B = reshape(A, size(A, 1), []) * W(:);

You could wrap the creation of weights into it's own function and make this generalizable to N weights. Since this uses recursion, N is limited to your current recursion limit (500 by default).

function W = createWeights(W, varargin)
    if numel(varargin) > 0
        W = createWeights(W(:) * varargin{1}(:).', varargin{2:end});
    end
end

And use it with:

W = createWeights(w1, w2, w3);
B = reshape(A, size(A, 1), []) * W(:);

Update

Using part of @CKT's very good suggestion to use kron, we could modify createWeights just a little bit.

function W = createWeights(W, varargin)
    if numel(varargin) > 0
        W = createWeights(kron(varargin{1}, W), varargin{2:end});
    end
end

Upvotes: 5

CKT
CKT

Reputation: 781

Again, you couldn't generalize this that well for N-D unless you made some function to construct the Kronecker product vector, but how about

  A = reshape(A, N, []) * kron(w3, kron(w2, w1));

Upvotes: 1

Related Questions