Replace for loop with matrix function, for function acting on matrix columns

2018 10 24
3 个回答
回答已采纳
更新时间：2018 10 24
5 次查看（30 天）

0 个投票

Let L be an n x n matrix, and let M be an n x p matrix. (In practice p << n/2, but n is of order 10^{5} or greater). I'd like to construct the n x p matrix N, with entries $N_{ij} = sum_k L_{ki}/M_{kj}$.

At present I'm using the for loop below. How can it be accelerated/avoided?

% Generate random matrices %
n = 100; p = 35;
L = rand(n,n);
M = rand(n,p);
%
% Produce N from L and M, where N_{ij} = sum_k Lki/Mkj %
N = inf(n,p);
for i = 1:n
  Li = L(:,i);
  for j = 1:p
      Mj = M(:,j);
      N(i,j) = sum(Li./Mj);
  end
end

请先登录，再回答此问题。

 采纳的回答

Stephen23
Stephen23 2018-10-24
编辑：Stephen23 2018-10-24

2 个投票

Using implicit scalar dimension expansion (version R2016b or later):
>> M = [1,2,2;4,5,4;8,10,1];
>> L = [1,2,3;1,2,3;1,2,3];
>> sum(permute(L,[2,3,1])./permute(M,[3,2,1]),3)
ans =
1.37500 0.80000 1.75000
2.75000 1.60000 3.50000
4.12500 2.40000 5.25000
For earlier versions replace the rdivide with bsxfun and rdivide.

更多回答（2 个）

Bruno Luong
Bruno Luong 2018-10-24

2 个投票

N = L.' * (1./M)

5 个评论
显示 3更早的评论 隐藏 3更早的评论

greengunpowder
greengunpowder 2018-10-24
Slick! However, I just tested the accuracy of the methods, and this results in an error of order 1e-11.
Bruno Luong
Bruno Luong 2018-10-24
编辑：Bruno Luong 2018-10-24
But do you know which method is more accurate?
greengunpowder
greengunpowder 2018-10-24
Not in general, because I've not explicitly calculated any solutions for big matrices. However, for the specific case with M and L as above an error arises in this method. (In the entry N(3,2).)
I did then assume that the discrepancy between solutions in other test cases would be down to the same (unknown) problem, which is not necessarily true...
(The for loop and the Stephen Cobeldick method have always produced the same N matrix in my tests.)
Bruno Luong
Bruno Luong 2018-10-24
The order of the summation is different when one uses SUM, or MTIMES and it also depends on the number of threads (so the CPU that MATLAB is running), which in return gives different summation results. So far the discrepancy is observed, but AFAIK no-one is able to question matrix MTIMES in a causual use cases.

请先登录，再进行评论。

madhan ravi
madhan ravi 2018-10-24
编辑：madhan ravi 2018-10-24

0 个投票

n = 100;
p = 35;
M = rand(n,p)
L = rand(n,n)
N = sum(L./M)

2 个评论
显示 无 隐藏 无

greengunpowder
greengunpowder 2018-10-24
sum(L./M) can't be calculated, because the matrix dimensions do not agree. If they did agree it would not the same function.
I would like entries of N, N_{ij} to be the sum of the element-wise division acting on the columns Li, Mj. See the difference below (apologies, can't work out how to copy content from my command line):
n = 3; p = 3;
M = [1 2 2; 4 5 4; 8 10 1];
L = [1 2 3; 1 2 3; 1 2 3];
%
sum(M./L) = [13.0000 8.5000 2.3333]
%
% Produce N from L and M, where N_{ij} = sum_k Lki/Mkj %
N = inf(n,p);
for i = 1:n
Li = L(:,i);
for j = 1:p
Mj = M(:,j);
N(i,j) = sum(Li./Mj);
end
end
%
N =
1.3750 0.8000 1.7500
2.7500 1.6000 3.5000
4.1250 2.4000 5.2500

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

产品

版本

R2018a

标签

Community Treasure Hunt

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

Start Hunting!

Translated by