CODEWITHSUNDEEP
matlabperformancematrixvectorvectorization
hipHopMetropolisHastings
hipHopMetropolisHastings

Reputation: 277

Vectorization of Matrix Quadratics in MATLAB

I am trying to "vectorize" this loop in Matlab for computational efficiency

for t=1:T
   j=1;
   for m=1:M
      for n=1:N
         y(t,j) = v{m,n} + data(t,:)*b{m,n} +  data(t,:)*f{m,n}*data(t,:)';
         j=j+1;
      end
   end
end

Where v is a (M x N) cell of scalars. b is a (M x N) cell of (K x 1) vectors. f is a (M x N) cell of (K x K) matrices. data is a (T x K) array.

To give an example of what I mean the code I used to vectorize the same loop without the quadratic term is:

B = [reshape(cell2mat(v)',1,N*M);cell2mat(reshape(b'),1,M*N)];
X = [ones(T,1),data];
y = X*B;

Thanks!

Upvotes: 1

Views: 87

Answers (2)

Divakar
Divakar

Reputation: 221614

Here's the most vectorized form targeted for performance -

% Extract as multi-dim arrays
vA = reshape([v{:}],M,N);
bA = reshape([b{:}],K,M,N);
fA = reshape([f{:}],K,K,M,N);

% Perform : data(t,:)*f{m,n} for all iterations
data_f_mult = reshape(data*reshape(fA,K,[]),T,K,M,N);

% Now there are three parts :
% v{m,n}
% data(t,:)*b{m,n}
% data(t,:)*f{m,n}*data(t,:)';

% Compute those parts one by one
parte1 = vA(:).';
parte2 = data*reshape(bA,[],M*N);

parte3 = zeros(T,M*N);
for t = 1:T
    parte3(t,:) = data(t,:)*reshape(data_f_mult(t,:,:),K,[]);
end

% Finally sum those up and to present in desired format permute dims
sums = bsxfun(@plus, parte1, parte2 + parte3);
out = reshape(permute(reshape(sums,T,M,N),[1,3,2]),[],M*N);

Upvotes: 0

hipHopMetropolisHastings
hipHopMetropolisHastings

Reputation: 277

For those interested here was the solution I found

f = f';
tMat = blkdiag(f{:})+(blkdiag(f{:}))';
y2BB = [reshape(cell2mat(v)',1,N*M);...
        cell2mat(reshape(b',1,M*N));...
        reshape(diag(blkdiag(f{:})),K,N*M);...
        reshape(tMat((tril(tMat,-1)~=0)),sum(1:K-1),M*N)]; 
y2YBar = [ones(T,1),data,data.^2];

jj=1;
kk=1;
ll=1;
for k=1:sum(1:K-1)
    y2YBar = [y2YBar,data(:,jj).*data(:,kk+jj)];
    if kk<(K-ll)
        kk=kk+1;
    else
        kk=1;
        jj=jj+1;
        ll=ll+1;
    end
end
y = y2YBar*y2BB;

Upvotes: 2

Related Questions