Triple weighted sum

Question

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.

do you need to generalise? because if so I'd rephrase your N, T1, T2, T3 as an array. — Tasos Papastylianou, Aug 23 '16 at 19:47
I actually just wanted the three dimensional case. But the generalization might be useful for someone else :) — phdstudent, Aug 23 '16 at 20:27
Thought you might be interested in benchmarks. For 1000 iterations performed in Octave, the looped version took 20 minutes, my version took 1.66s with repmat and 0.92s without, and Suever's took 0.77s. :) — Tasos Papastylianou, Aug 23 '16 at 22:47
PS. and 0.40s in Julia (with my method adapted to Julia :p ) — Tasos Papastylianou, Aug 23 '16 at 23:19
Yes, at some point I should learn julia. I program in Matlab and Fortran,learning Julia should not have a huge fixed cost. — phdstudent, Aug 23 '16 at 23:20

Suever · Accepted Answer · 2016-08-23T19:57:56.710

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

Tasos Papastylianou · Answer 2 · 2016-08-23T22:24:49.990

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)

score 1 · Answer 3 · answered Aug 23 '16 at 19:52

1

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));

answered Aug 23 '16 at 19:52

CKT

781
5
6

score 1 · Answer 4 · answered Aug 24 '16 at 15:45

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.

Triple weighted sum

4 Answers4