how to write the right code to measure the elapsed time of this method d1 = diag(1./diag(A))

how to write the right code to measure the elapsed time of this method

d1 = diag(1./diag(A))

Answer 1

[John D'Errico]

Why do you care? The computations will be almost inifinitessimal, even for a huge matrix. (Do you understand the notation O(n)?) There are only n operations to be performed here, so O(n). The only cost as the matrix gets large will be allocating and storing the array, almost all of which will be zeros.

A = rand(10000);
timeit(@() diag(1./diag(A)))

ans = 0.0141

More important may be what are you doing with that matrix! The result will be a diagonal matrix. And almost the only reason you are doing that is to multiply it with another matrix. And there are better ways to do that operation! The point is, a matrix multiply scales the rows (or columns as you do it) of the final array. But a matrix multiply is the expensive part. For example...

A = rand(5000);
B = rand(5000);
Adiag = 1./diag(A);
d = diag(1./diag(A));
timeit(@() B*d)

ans = 0.4235

timeit(@() B.*Adiag.')

ans = 0.0220

Why are the two results that much different? THINK ABOUT IT! What is the time cost for a matrix multiply of two nxn matrices? If you don't know, then you really need to do some reading. Hint: O(n^3)

But all you want to do is perform n^2 multiplies. What is the time cost of the second form I show? HInt: O(n^2)

My point is, if you are asking this question, it is most likely you are asking the wrong question. Just because you CAN scale the rows or coluns of matrix by a constant does not mean that a good way to accomplish the task by using a matrix multiply with a diagonal matrix. And just because you see it written in a formula does not mean you should always write code that looks just like your equation. At least, that is if you care about speed.

[Matt J]

If you want the convenience of scaling row and columns with matrix-multiplication notation, using spdiag() might be a tolerable compromise,

A = rand(8000);
B = rand(8000);
Adiag = 1./diag(A);
d = diag(1./diag(A));
ds= diagsp(1./(diag(A)));

timeit(@() diag(1./diag(A)))

ans = 0.0131

timeit(@() diagsp(1./diag(A)))

ans = 1.7778e-04

timeit(@() B*d)

ans = 2.8423

timeit(@() B*ds)

ans = 0.2889

timeit(@() B.*Adiag.')

ans = 0.0712

function D=diagsp(v)

  N=numel(v);
  D=spdiags(v,0,N,N);

end

Answer 2

[David Goodmanson]

Hi Amna,

A common way to do this is with tic and toc. So:

a = rand(1e4,1e4);
tic
b = diag(1./diag(a));
toc
  Elapsed time is 0.018291 seconds.

Not so much with this particular calculation, but there tends to be variation in the time reported by tic and toc. To average this out, you can put the calculation into a do-nothing-else for loop (being careful not to put the calculation of 'a' inside the loop):

a = rand(1e4,1e4);
tic
for k = 1:200
    b = diag(1./diag(a));
end
toc
Elapsed time is 3.527247 seconds.

i.e. .0176 sec per operation