Multiplication of very large matrix

Hello everyone,

Could anyone help me in re-writing the following equation to make the processing faster. N here is the number of points and it could be 50,000. D is a diagonal matrix where only the diagonal contains values and the other elements are zeros.

 M=eye(N)-((1./max((ones(N)'*D*ones(N)),eps))*(ones(N)*ones(N)'*D));

[Rik]

Why do you have all those multiplications with ones(N) in there? Most of them don't actually do anything.

Can you confirm D is N by N?

[MA]

Yes. D is N by N.

[David Hill]

Looks like just:

M=eye(N)-1./D^2;

[Jan]

I'm curious also about the multiplications with ones() . The purpose of transposing a ones() is magic also.

Please post some realistic test data, such that the reader have a chance to understand, what you are calculating.

[MA]

I am trying to implement equation 10 in the following paper:

https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4711055

I recieved different results when ignoring the onces( ).

[David Goodmanson]

Hi MA,

for those of us not in IEEE, the reference is not of much help. The

(ones(N)*ones(N)'*D))

factor is in the numerator of the rhs, is that correct? Also, what is meant by ones(N)? In Matlab it's an NxN matrix, but in context it makes at least as much sense if it were an Nx1 column vector.

[MA]

Thanks very much for your help guys. I now know how to simplify it. So, if anyone is interesting. Here is how to simplify it:

M=eye(N)-repmat(((N.*N.*diag(D))./sum(diag(D)))',[N,1]); 

[MA]

You are right David Goodmanson. It is a vector of size N by 1. So, multiplying by its transpose make sense. Thanks alot for the help.

Answer 1

[Jan]

How strange: The comments above have not been displayed on my other computer. So this answer was written 1 hour after MA's comment.

In ones(N)*ones(N)' the transposition is completely meaningless. The result can be obtained much cheaper by: repmat(N, [N, N]). The matrox multiplication with this matrix can be formulated much cheaper:

The multiplication by () is an extremly expensive way to calculate:

A = ones(N) * ones(N)' * D
B = sum(D, 1) * N

Now B is a row vector only, but A is only a blownup version with N repetitions.

Because D is a diagonal matrix, you can omit the sum also.

Equivalently for ones(N)'*D*ones(N) : Of course you cannot simply omit the multiplication, but you can express it much cheaper by: sum(D(:)) or cheaper: sum(diag(D)) .

DD = diag(D);
M  = eye(N) - N^2 * DD.' ./ max(sum(DD), eps);

Note, that this matrix is extremely redundant: All columns contain the same value except for the diagonal, which is 1 smaller. An efficient implementation would exploit this instead of creating a large matrix.