Different run time of eig() for the same matrix (calculated in 2 different ways)

I hope, a good answer will help understanding of how to use MATLAB in a truly efficient way...

From the same input (stored in input.mat in the attachment), I calculate the matrix M in 2 different ways and then calculate its evals and evecs.

Case 1:

clear;
load('input.mat');
D = d.^(-1/2);
D = diag(D);
M = D*L*D;
[V, v] = eig(M);

Case 2:

clear;
load('input.mat');
d = sqrt(d);
D = d * d';
M = L./D;
[V, v] = eig(M);

Clearly, the 2nd way is more efficient, BUT the resulting matrix M is the same for the 2 cases. So the eig() function should probably have the same running time in the both cases. However, for the first case the run time of eig() is about 12.5 s, and for the second case the run time of eig() is about 1.8 s (the numbers were taken from the profile results).

The input matrix L is a matrix with many zeros (but it's a full matrix). The input vector d is such that isequal(diag(L),d) is true.

I've uploaded the inputs because the above mysterious behaviour is apparently not present for a randomly generated input square matrix.

I guess, it might have something to do with memory efficiency or cache usage, but it would be interesting to know a more precise reason.

[the cyclist]

I don't see the uploaded matrix.

[Ilya]

Yes, indeed, I've uploaded now. By the way, I think I now know why.. issymmetric(M) is true for the 2nd case, but false in the 1st case. Due to very small numerical errors (~1e-16), M fails to qualify as a symmetric matrix in the 1st case..

Answer 1

[John D'Errico]

I think what you are missing is that in the second case, it turned out that M was EXACTLY symmetric.

norm(M1-M2) ans = 2.4365e-14

norm(M1-M') ans = 7.1248e-15

norm(M2-M2') ans = 0

I named the two versions of M as M1 and M2.

See that the latter one was symmetric, whereas the former, not.

MANY times in MATLAB the code can take advantage of symmetry. And ALMOST symmetric is not symmetric. It looks like you saw the difference when I read your comments. But MATLAB does not know how close to a symmetric matrix it must be to consider it symmetric. So the code will use the symmetric branch only when that is exactly true.

[Ilya]

Yes, I've printed eigenvalues and saw that many are complex, which should not be the case for a symmetric matrix. To conclude this question, it would be nice to have a function that gives a metric of "proximity" to symmetry. Then if matrix is very near (e.g. only because of numeric roundoff), it's possible to copy the lower triangular part and obtain a significant boost in eig() speed...

[Ilya]

OK, to conclude my question, I've found that there is no such obvious metric for "near symmetry" implemented in MATLAB, so one way is to use something like this which is not very fast (i.e. it distorts the "true calculation time" of an algorithm). Alternatively it's possible to copy the lower triangular part as M = triu(M,1) + diag(diag(M)) + triu(M,1).' if ~issymmetric(M) but it's expected to be symmetric.

[Walter Roberson]

For real square matrices, (M+M.')/2 is symmetric

[Ilya]

@Walter: Thanks, this is a good (fast) one