Matlab 编码器导致 bsxfun 中的矩阵扩展问题

Posted 2023-02-16

【中文标题】Matlab 编码器导致 bsxfun 中的矩阵扩展问题

【英文标题】：Matlab coder causing problems with matrix expansion in bsxfun

【发布时间】：2016-11-28 22:43:08

【问题描述】：

我一直在尝试从一个本科生为我们编写的代码中矢量化一些操作。他使用了丑陋的嵌套循环，我试图将其更改为 bsxfun 调用。当我按原样使用此函数时，此方法有效，但当我尝试使用编码器对其进行 mex 编译时，出现错误。

编辑：我意识到我应该澄清代码的某些部分。 S、W 和 N 是表示成对比较 ij 的偏好的 nxn 矩阵。 mu 是包含绝对评分重建的（预测）分数的行向量。

要检查代码的作用，您可能需要检查页面的最底部。

原代码如下：

for i = 1:n
    for j = i:n
        if i~= j
            mu = x(i) - x(j);
            md1 = normcdf(-delta1, mu); %minus d0
            md0 = normcdf(-delta0, mu);
            d0 = normcdf( delta0, mu);
            d1 = normcdf( delta1, mu);

            Y = Y + S(i,j) .* log(1  - d1) + ...
                    W(i,j) .* log(d1 - d0) + ...
                    N(i,j) .* log(d0 - md0) + ...
                    W(j,i) .* log(md0 - md1) + ...
                    S(j,i) .* log(md1);

        end
    end
end
Y = -Y; %this will be used in fminsearch, therefore the negative sign

我为此写了一个快速而肮脏的矢量化解决方案：

diffs = bsxfun(@minus,x,x');

% create the deltas array augmented by -Inf and Inf for easy diff calculation
deltass = sort([-Inf, Inf, deltas, -1*deltas]);
deltass = reshape(deltass, [1, 1, length(deltass)]);
normcdfs = diff(1-bsxfun(@normcdfupper, diffs, deltass), 1, 3);
normcdfs(repmat(logical(tril(ones(size(normcdfs(:,:,1))))),[1,1,length(deltass)-1])) = 1;

Y = -sum(sum(sum(normcdfs_input.*log(normcdfs))));

编辑：norcdfupper 是我编写的一个包装类，它调用 normcdf(diff, delta, 'upper') 因为我遇到了问题，尾部可能会导致错误。

现在就像我说的，这在我的日常使用中效果很好。但是由于混合旧的嵌套 for 循环版本带来了巨大的加速，我尝试对这个加速版本做同样的事情，只是在 coder 中遇到以下错误：

Expansion is only supported along dimensions where one input argument or the other has a fixed length of 1.

Error in lik_ml (line 27)
normcdfs = diff(1-bsxfun(@normcdfupper,diffs,deltass),1,3);

diffs 是一个 :infx:inf 双精度数，而 deltass 是一个 1x1x:inf 双精度数。 Coder 告诉我：

diffs deltass

如果有人可以帮助我解决此错误消息，我将不胜感激。

编辑：一些代码来测试它：

clear all;

% rng(322);
delta0 = 1;
delta1 = 2;
deltas = [delta0, delta1];

sigma = 0.5;
options = 5;

maxVotes = 10000;
voteStep = 3;

initialVoteStep = 3;

raw_mu = rand(1,options);
mu = sort(zscore(raw_mu));

S = zeros(options);
W = zeros(options);
N = zeros(options);

for i = 1:options
    for j = i:options
        if i ~= j
            S(i,j) = 1 - normcdf(delta1, mu(i)-mu(j));
            W(i,j) = normcdf( delta1, mu(i)-mu(j)) - normcdf( delta0, mu(i)-mu(j));
            N(i,j) = normcdf( delta0, mu(i)-mu(j)) - normcdf(-delta0, mu(i)-mu(j));
            W(j,i) = normcdf(-delta0, mu(i)-mu(j)) - normcdf(-delta1, mu(i)-mu(j));
            S(j,i) = normcdf(-delta1, mu(i)-mu(j));
        end
    end
end

diffs = bsxfun(@minus,mu,mu');
deltass = sort([-Inf, Inf, deltas, -1*deltas]);
deltass = reshape(deltass,[1,1,length(deltass)]);
normcdfs = diff(bsxfun(@(x,y) normcdf(y,x),diffs,deltass),1,3);

【问题讨论】：

你在哪里定义x？

哦，对不起，我稍微修改了代码。 x 的意思是 mu。我将编辑上面的评论。

【参考方案1】：

这里是for循环的矢量化版本

[I,J] = meshgrid(1:n);
idx = I<J;
I = I(idx);
J = J(idx);
IJ = sub2ind([n n],I,J);
JI = sub2ind([n n],J,I);
mu = x(I)-x(J);
md1 = normcdf(-delta1, mu).'; 
md0 = normcdf(-delta0, mu).';
d0 = normcdf( delta0, mu).';
d1 = normcdf( delta1, mu).';
Y = sum(S(IJ) .* log(max(0,1  - d1)) + ...
W(IJ) .* log(max(0,d1 - d0)) + ...
N(IJ) .* log(max(0,d0 - md0)) + ...
W(JI) .* log(max(0,md0 - md1)) + ...
S(JI) .* log(max(0,md1)));

【讨论】：

您的代码存在一些问题。首先，有一个矩阵乘法错误：size(S(IJ)) ans = 20 1 size(log(1-d1)) ans = 1 20 所以我改变了要转置的 normcdfs 的维度。但随后结果不同。使用 rng(322) 作为我原始帖子底部的代码 sn-p 并将相应的变量输入到您的原始解决方案中：lik5_ml_delta_known(mu,S,W,N,delta0,delta1) ans = 10.3631 lik5_ml_delta_known3(mu,S,W,N,delta0,delta1) ans = 17.6042 :(

此外，我遇到的真正问题是理解为什么我的实现不能被混合。给出的错误消息不足以让我理解它。

@FranzHahn 向量化版本是正确的，除了应该应用您提到的一个 traspose。结果是不同的，因为循环版本与向量化版本随机数的生成方式不同。在矢量化版本中，前 20 个随机数存储在 md1 中，但在循环版本中，前 4 个随机数存储在 md1 md0 d0 d1 然后是第二个 4 个数字 ...

对不起，我可能是盲人，但我看不出随机数在这里是如何起作用的。我们只是在给定模型（具有均值 mu 的正态分布）的情况下计算数据的对数似然（S/W/N 矩阵）。在计算中没有调用随机数生成器，它只发生在数据生成中，我输入的变量完全相同。

我相信你的实现中的错误在于原来的双重嵌套 for 循环从 i=1:options 和 j=**i**:options 运行。您的实施中似乎没有考虑到这一点。