Is it possible to vectorize this 'for loop' involving multiple matrix inversions

조회 수: 6 (최근 30일)
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.

채택된 답변

James Tursa
James Tursa 2021년 8월 5일

추가 답변 (1개)

Walter Roberson
Walter Roberson 2021년 8월 6일
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))
ans = '@(cov_Ip_b,cov_Ip_g,cov_Ip_r,eps,var_I_bb,var_I_gb,var_I_gg,var_I_rb,var_I_rg,var_I_rr)[(cov_Ip_r.*eps.^2-cov_Ip_r.*var_I_gb.^2-cov_Ip_b.*eps.*var_I_rb+cov_Ip_r.*eps.*var_I_bb-cov_Ip_g.*eps.*var_I_rg+cov_Ip_r.*eps.*var_I_gg+cov_Ip_b.*var_I_gb.*var_I_rg-cov_Ip_b.*var_I_gg.*var_I_rb-cov_Ip_g.*var_I_bb.*var_I_rg+cov_Ip_g.*var_I_gb.*var_I_rb+cov_Ip_r.*var_I_bb.*var_I_gg)./(eps.^2.*var_I_bb-eps.*var_I_gb.^2+eps.^2.*var_I_gg-eps.*var_I_rb.^2-eps.*var_I_rg.^2+eps.^2.*var_I_rr-var_I_bb.*var_I_rg.^2-var_I_gg.*var_I_rb.^2-var_I_gb.^2.*var_I_rr+eps.^3+eps.*var_I_bb.*var_I_gg+eps.*var_I_bb.*var_I_rr+eps.*var_I_gg.*var_I_rr+var_I_bb.*var_I_gg.*var_I_rr+var_I_gb.*var_I_rb.*var_I_rg.*2.0),(cov_Ip_g.*eps.^2-cov_Ip_g.*var_I_rb.^2-cov_Ip_b.*eps.*var_I_gb+cov_Ip_g.*eps.*var_I_bb+cov_Ip_g.*eps.*var_I_rr-cov_Ip_r.*eps.*var_I_rg-cov_Ip_b.*var_I_gb.*var_I_rr+cov_Ip_b.*var_I_rb.*var_I_rg+cov_Ip_g.*var_I_bb.*var_I_rr-cov_Ip_r.*var_I_bb.*var_I_rg+cov_Ip_r.*var_I_gb.*var_I_rb)./(eps.^2.*var_I_bb-eps.*var_I_gb.^2+eps.^2.*var_I_gg-eps.*var_I_rb.^2-eps.*var_I_rg.^2+eps.^2.*var_I_rr-var_I_bb.*var_I_rg.^2-var_I_gg.*var_I_rb.^2-var_I_gb.^2.*var_I_rr+eps.^3+eps.*var_I_bb.*var_I_gg+eps.*var_I_bb.*var_I_rr+eps.*var_I_gg.*var_I_rr+var_I_bb.*var_I_gg.*var_I_rr+var_I_gb.*var_I_rb.*var_I_rg.*2.0),(cov_Ip_b.*eps.^2-cov_Ip_b.*var_I_rg.^2+cov_Ip_b.*eps.*var_I_gg-cov_Ip_g.*eps.*var_I_gb+cov_Ip_b.*eps.*var_I_rr-cov_Ip_r.*eps.*var_I_rb+cov_Ip_b.*var_I_gg.*var_I_rr-cov_Ip_g.*var_I_gb.*var_I_rr+cov_Ip_g.*var_I_rb.*var_I_rg+cov_Ip_r.*var_I_gb.*var_I_rg-cov_Ip_r.*var_I_gg.*var_I_rb)./(eps.^2.*var_I_bb-eps.*var_I_gb.^2+eps.^2.*var_I_gg-eps.*var_I_rb.^2-eps.*var_I_rg.^2+eps.^2.*var_I_rr-var_I_bb.*var_I_rg.^2-var_I_gg.*var_I_rb.^2-var_I_gb.^2.*var_I_rr+eps.^3+eps.*var_I_bb.*var_I_gg+eps.*var_I_bb.*var_I_rr+eps.*var_I_gg.*var_I_rr+var_I_bb.*var_I_gg.*var_I_rr+var_I_gb.*var_I_rb.*var_I_rg.*2.0)]'

카테고리

Help CenterFile Exchange에서 Logical에 대해 자세히 알아보기

제품


릴리스

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by