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

3 ビュー (過去 30 日間)
Benjamin Tilbury 2021 年 8 月 5 日

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 2021 年 8 月 5 日
##### 1 件のコメント表示非表示 なし
Benjamin Tilbury 2021 年 8 月 6 日
Brilliant. This is perfect.

サインインしてコメントする。

### その他の回答 (1 件)

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)]'

サインインしてコメントする。

R2021a

### Community Treasure Hunt

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

Start Hunting!