2

让x = [1,...,t]是一个带有t组件A和P数组的向量。我问自己是否有机会缩短这个，因为它看起来很麻烦：

for n = 1:t
    for m = 1:n
        H(n,m) = A(n,m) + x(n) * P(n,m)
    end
end
4

5 回答 5

4

我的建议：bsxfun(@times,x,P) + A;

例如

A = rand(3);
P = rand(3);
x = rand(3,1);

for n = 1:3
    for m = 1:3
        H(n,m) = A(n,m) + x(n) * P(n,m);
    end
end 


H2 = bsxfun(@times,x,P) + A;

%//Check that they're the same
all(H(:) == H2(:))

返回

ans =  1

编辑： 阿姆罗是对的！使第二个循环取决于第一次使用tril：

H2 = tril(bsxfun(@times,x,P) + A);

顺便说一句，矩阵是平方的吗，因为这也会产生其他问题

于 2013-09-20T12:46:04.683 回答
1

就像我在评论中指出的那样，除非是拼写错误，否则第二个 for 循环计数器取决于第一个 for 循环的计数器......

如果这是故意的，我想出了以下解决方案：

% some random data
t = 10;
x = (1:t)';
A = rand(t,t);
P = rand(t,t);

% double for-loop
H = zeros(t,t);
for n = 1:t
    for m = 1:n
        H(n,m) = A(n,m) + x(n) * P(n,m);
    end
end

% vectorized using linear-indexing
[a,b] = ndgrid(1:t,1:t);
idx = sub2ind([t t], nonzeros(tril(a)), nonzeros(tril(b)));
xidx = nonzeros(tril(a));
HH = zeros(t);
HH(tril(true(t))) = A(idx) + x(xidx).*P(idx);

% check the results are the same
assert(isequal(H,HH))

我更喜欢@Dan的解决方案。这里唯一的优点是我不计算不必要的值（因为矩阵的上半部分为零），而另一个解决方案计算完整的矩阵，然后减少额外的东西。

于 2013-09-20T13:08:24.363 回答
1

tril(A + P.*repmat(x',1,t))

编辑。这适用于当 x 是行向量时。如果 x 是列向量，则使用tril(A + P.*repmat(x,t,1))

如果您的示例代码是正确的，那么对于任何 j > i，H(i,j) = 0，例如 X(1,2)。例如t = 3，你会有。

H = 

'A(1,1) + x(1) * P(1,1)'                          []                          []
'A(2,1) + x(2) * P(2,1)'    'A(2,2) + x(2) * P(2,2)'                          []
'A(3,1) + x(3) * P(3,1)'    'A(3,2) + x(3) * P(3,2)'    'A(3,3) + x(3) * P(3,3)'
于 2013-09-20T13:09:31.767 回答
0

由于第二个循环的范围是1:n，您可以取矩阵的下三角部分A并P进行计算

H = bsxfun(@times,x(:),tril(P)) + tril(A);
于 2013-09-23T08:24:29.500 回答
0

一个好的开始将是

H = A + x*P

这可能不是一个可行的解决方案，您必须检查数组和向量的一致性，并确保您使用的是正确的乘法，但这应该足以为您指明正确的方向。如果您是 Matlab 新手，请注意向量可以是1xn或nx1，即行向量和列向量是不同的种类，这与许多编程语言不同。如果x不是你想要的 rhs，你可能想要它的转置，x'在 Matlab 中。

从某个角度来看，Matlab 是一种数组语言，显式循环通常是不必要的，而且通常甚至不是一个好方法。

于 2013-09-20T12:36:30.470 回答