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

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 

Answer 1

[Stephen23]

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.

[madhan ravi]

+1 like a rockstar @stephen

[greengunpowder]

Excellent. Thank you :-)

Answer 2

[Bruno Luong]

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

[Stephen23]

+1 best use of MATLAB

[greengunpowder]

Slick! However, I just tested the accuracy of the methods, and this results in an error of order 1e-11.

[Bruno Luong]

But do you know which method is more accurate?

[greengunpowder]

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]

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.

Answer 3

[madhan ravi]

 n = 100; 

 p = 35;

M = rand(n,p)

L = rand(n,n)

N = sum(L./M)

[madhan ravi]

You don't need loop to do this operation

[greengunpowder]

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