Is it possible to vectorize this 'for loop' involving multiple matrix inversions
8 views (last 30 days)
Show older comments
Benjamin Tilbury
on 5 Aug 2021
Answered: Walter Roberson
on 6 Aug 2021
I have a vendetta against for loops.
I'm trying to remove them all from a particular project but I've hit a roadblock.
Below is part of the code for an implementation of the color fast guided filter:
% eps is scalar
% All other variables are 1xN vectors
a = zeros(N, 3);
for q = 1:N
sigma = [var_I_rr(q), var_I_rg(q), var_I_rb(q);
var_I_rg(q), var_I_gg(q), var_I_gb(q);
var_I_rb(q), var_I_gb(q), var_I_bb(q)];
cov_Ip = [cov_Ip_r(q), cov_Ip_g(q), cov_Ip_b(q)];
a(q,:) = cov_Ip / (sigma + eps * eye(3));
% Alternately
% a(q,:) = cov_Ip * inv(sigma + eps * eye(3));
end
In essence we've got a 1x3 vector being multiplied by the inverse of a 3x3 matrix, n number of times.
The following attempt obiously doesn't work because 'A/B' or 'A*inv(B)' is only defined for 2-D matrices, but hopefully it illustrates what I'm trying to do.
top = cat(2,reshape(var_I_rr,1,1,[]),reshape(var_I_rg,1,1,[]),reshape(var_I_rb,1,1,[]));
mid = cat(2,reshape(var_I_rg,1,1,[]),reshape(var_I_gg,1,1,[]),reshape(var_I_gb,1,1,[]));
bot = cat(2,reshape(var_I_rb,1,1,[]),reshape(var_I_gb,1,1,[]),reshape(var_I_bb,1,1,[]));
sigma = cat(1,top,mid,bot); % I'm sure there's a more efficient way to stack like this but it's not really important
cov_Ip = cat(2,reshape(cov_Ip_r,1,1,[]),reshape(cov_Ip_g,1,1,[]),reshape(cov_Ip_b,1,1,[]));
eps = repmat(eps*reshape(eye(3),3,3,1),1,1,N);
% sigma is now a 3x3xN matrix
% eps is now a 3x3xN matrix
% cov_Ip is now a 1x3xN matrix
b = cov_Ip / (sigma+eps); % Does not work.
If sigma+eps was already inverted (along the first two dimensions) I could do the following:
% invSigma is a 3x3xN matrix
b = squeeze(pagemtimes(cov_Ip, invSigma))';
To illustrate the problem more clearly, here's simpler example of what I'm effectively trying to vectorize.
% A is a 1 x 3 x N matrix
% B is a 3 x 3 x N matrix
% C is a N x 3 matrix
for i = 1:N
a = A(:,:,i); % 1 x 3
b = B(:,:,i); % 3 x 3
c = a / b; % 1 x 3
% c = a * inv(b);
C(i,:) = c;
end
Is there any way to vectorize this. Is it possible to invert the first two dimensions of a 3-D matrix independently over the last dimension without using a for loop? Like 'pagemtimes' but for division/inverting.
0 Comments
Accepted Answer
More Answers (1)
Walter Roberson
on 6 Aug 2021
q = 1;
syms eps
syms cov_Ip_b cov_Ip_g cov_Ip_r
syms var_I_bb
syms var_I_gb var_I_gg
syms var_I_rb var_I_rg var_I_rr
sigma = [var_I_rr(q), var_I_rg(q), var_I_rb(q);
var_I_rg(q), var_I_gg(q), var_I_gb(q);
var_I_rb(q), var_I_gb(q), var_I_bb(q)];
cov_Ip = [cov_Ip_r(q), cov_Ip_g(q), cov_Ip_b(q)];
a(q,:) = cov_Ip / (sigma + eps * eye(3));
func2str(matlabFunction(a))
0 Comments
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!