提高求解方程的代码性能

Question

我对 MATLAB 很陌生，我正在努力寻找这个方程的解。

其中矩阵的维度是 {λ} = N×K、{Y} = N×D、{π} = 1×K 和 {μ} = D×K。

我创建的东西看起来更像是一个怪物，而不是高效的 MATLAB 代码，这是由于严重缺乏 MATLAB 技能。

实际上，我决定一步一步解决这个问题，以便让自己熟悉 MATLAB 的工作原理。我最终得到了一个效率极低的代码，我什至无法将它组合起来。甚至最终结果也是错误的，因为我想创建一个 NxK 矩阵。任何指导将不胜感激。

%Dimensions:
nn = 10;
dd = 7;
kk = 5;

%Initial variables:
lambda0 = rand(nn,kk);
sigma = rand(1);
Y = rand(nn,dd);
mu = rand(dd,kk);

%Calculate pies:
First_part = log(pie(:)./(1.-pie(:))); % Kx1 vector

%Calculate Second part:
for n = 1:nn
for d = 1:dd
lambda_mu(d,:) = lambda0(n,:).*mu(d,:); %lambdamu is a DxK matrix
end
lambda_mu2(n,:) = sum(lambda_mu,2); %This is a NXD matrix
end

%Y-lambdamu2:
for n = 1:nn
YY(n,:) = Y(n,:)-lambda_mu2(n,:); % This is a NxD vector
end

%YY*mu:
Second = (YY*mu)./sigma^2; % NxK

%Third:
Third = (mu'*mu)./2*sigma^2; % KxK

%Final:
for n=1:nn
Final_part = transpose(First_part(:))+Second(n,:)-Third;
end

Update1：好吧，我欠自己一杯啤酒。 所以，我通过如下修改代码来创建 {λ} = N×K 矩阵取得了一些进展：

更新[2]：花太多时间编程，你会错过细节。我修复了错误，因为我设置了 {Y} = D×K，而不是 {Y} = D×N。

%Dimensions:
nn = 10;
dd = 7;
kk = 5;

%Initial variables:
lambda0 = rand(nn,kk);
sigma = rand(1);
Y = rand(dd,nn);
mu = rand(dd,kk);
pie = rand(1,kk);
sigma = rand(1);

%Calculate the equation.
for n = 1:nn
for k = 1:kk
    lambda(n,k) = log(pie(k)/(1-pie(k))) + (transpose(Y(:,n)-... 
        (sum(lambda0(n,1:end ~= k).*mu(:,1:end ~= k),2)))*mu(:,k))/...
        sigma^2 - transpose(mu(:,k))*mu(:,k)/2*sigma^2;
end
end

但是，现在的问题是它只能达到 n=5，例如，当我尝试循环 n=6 时，我收到以下消息：“索引超出矩阵维度”，所以实际上我只得到一个 5×5 矩阵。请给点建议？

P.S：我什至尝试将 lambda 从 lambda(n,k) 更改为 lambda(k,n) 但结果是一样的。

Answer 1

为了使尺寸工作，乘以ones(nn,kk) 和ones(1,kk)

nn = 10;
dd = 7;
kk = 5;

%Initial variables:
lambda0 = rand(nn,kk);
sigma = rand(1);
Y = rand(nn,dd);
mu = rand(dd,kk);
pie = rand(kk,1);

%Calculate pies:
First_part = ones(nn,kk)*log(pie./(1.-pie))*ones(1,kk); % Kx1 vector

%Calculate Second part:
for n = 1:nn
    for d = 1:dd
        lambda_mu(d,:) = lambda0(n,:).*mu(d,:); %lambdamu is a DxK matrix
    end
    lambda_mu2(n,:) = sum(lambda_mu,2); %This is a NXD matrix
end

%Y-lambdamu2:
for n = 1:nn
    YY(n,:) = Y(n,:)-lambda_mu2(n,:); % This is a NxD vector
end

%YY*mu:
Second_part = (YY*mu)./sigma^2; % NxK

%Third:
Third_part = ones(nn,kk)*(mu'*mu)./2*sigma^2; % KxK


%Final:
for n=1:nn
    Final_part(n,:) = First_part(n,:)+Second_part(n,:)-Third_part(n,:);
end

编辑：我的想法更像是：

%Initial variables:
lambda = rand(nn,kk);
sigma = rand(1);
y = rand(nn,1);
mu = rand(nn,kk);
pie = rand(1,kk);

%Calculate pies:
First_part = ones(nn,1)*log(pie./(1-pie)); % NxK vector

避免循环。

也许使用类似的东西：

ones(kk)-diag([ones(kk,1)])

当对i 求和时不等于j，或者更确切地说使用第三维来表示j。