Why does MATLAB calculates when using speye?

Hello,

i just noticed a weird thing

>> tic;speye(10000)*speye(10000);toc
Elapsed time is 0.000753 seconds.

>> tic;speye(100000)*speye(100000);toc
Elapsed time is 0.007663 seconds.

>> tic;speye(1000000)*speye(1000000);toc
Elapsed time is 0.073734 seconds.

>> tic;speye(10000000)*speye(10000000);toc
Elapsed time is 0.737171 seconds.

Where is the sense in using an implemented sparse unit matrix for calculation? Same goes for a multiplication between a dense matrix and a sparse unit matrix.

>> tic;eye(1000)*speye(1000);toc
Elapsed time is 0.032003 seconds.
>> tic;eye(10000)*speye(10000);toc
Elapsed time is 1.331953 seconds.

and between zero matrices ...

>> tic;sparse(10000,10000)*sparse(10000,10000);toc
Elapsed time is 0.000416 seconds.
>> tic;sparse(100000,100000)*sparse(100000,100000);toc
Elapsed time is 0.001833 seconds.
>> tic;sparse(1000000,1000000)*sparse(1000000,1000000);toc
Elapsed time is 0.012763 seconds.
>> tic;sparse(10000000,10000000)*sparse(10000000,10000000);toc
Elapsed time is 0.111892 seconds.
>> tic;eye(1000,1000)*sparse(1000,1000);toc
Elapsed time is 0.013433 seconds.
>> tic;eye(10000,10000)*sparse(10000,10000);toc
Elapsed time is 0.753174 seconds.

I thought this datatype is also designed to avoid such unnecessary performance lacks.

This does really matter for me, because in my calculation some matrices will be unit matrices in special cases. So i created an independent calculation function which does not care about the matrices, it only calculates them. In the other part i check the conditions and set some matrices to unit sparse matrices, because i want to avoid the calculation time and send them to the calculation function. Unfortunately this does not work as expected and i get an unnecessary performance lack.

Answer 1

[Jan]

You do not measure the time for the multiplication, but the time for the creation of the matrix is included also. Please post the reuslts of:

tic;speye(10000)*speye(10000);toc
A = speye(10000);
tic; A * A; toc

tic;speye(10000000)*speye(10000000);toc
A = speye(10000000);
tic; A * A; toc

The sense of the sparse arrays is that e.g. the multiplication considers the 0 elements and does not perform a multiplication for them. In addition the storage of sparse matrices is more dense than a full matrix, such that speye(10000000) is possible even when eye(10000000) would cause an out-of-memory already. It is not the idea of sparse matrices to short-cut multiplications with the identity matrix and it is easy to avoid such actions programmatically.

[wma]

tic;speye(10000)*speye(10000);toc
Elapsed time is 0.000496 seconds.

A = speye(10000);
tic; A * A; toc
Elapsed time is 0.000406 seconds.

tic;speye(10000000)*speye(10000000);toc
Elapsed time is 0.808250 seconds.

A = speye(10000000);
tic; A * A; toc
Elapsed time is 0.380131 seconds.

My complain is that speye(x) and speye(y) have should never differ in creation/multiplication/inverse time because it is wasted time and memory. Same thing for completely empty matrices.

If A and speye is multiplied like

   B = A * speye(x);

then B should be set as a pointer to A without doing any multiplication. There is also no need for a copy of A as long as no manipulation was done in B. Many functions like inv, sqrt should be implemented to avoid wasting memory and time by using simple math rules instead of silly one by one multiplication.

If two speye(x) and speye(x) were multiplied there is also no need to calculate something, the result is speye(x) in any size!

If one of the matrices is sparse(x,y) the result for the whole term A*B*C*D*...*Z*sparse(x,y) = sparse(size(A,1),y)! There is also no need to waste time at this calculation.

[Jan]

But this is not the intention of sparse arrays, neither in Matlab nor when you use any other sparse library (at least I do not know any). While the sparse matrices are optimized for a giantic number of different possible values, you are talking about 2 exceptions: The empty and the identity matrix.

It is very easy to carry an extra variable as flag to catch these 2 rare exceptions:

if rand > 0
  A = speye(1000);
  AisEye = true;
else
  A = rand(1000);
  AisEye = false;
end
...
if AisEye
  B = A;
  BisEye = true;
else
  B = A * A;
  BisEye = false;
end

etc. If you want to avoid the computations performed by the strange and rare exception inv(speye()), better avoid calling inv() expesically on identity matrices. If you need this for your work, an oop class, which considers such exceptions, would be easy to implement. Therefore I do not see a reason to complain.

[wma]

It is a bad code guide to turnover common situations to the end user and waste intelligent information.

To begin a multiplication with an speye(x) or sparse(x,y) in the chain is bad design from the root and should be fixed there.

"... (at least I do not know any) ..."

Yes of course, because it is still not common to calculate with sparse types. Your example needs 14 lines for two matrices. What is about

km = -(B*Q*B')\(A*x + B*Q*R'*kr - B*v + psi)

this common least squares calculation, where Q can be speye(b) and B can be something like

B = [ speye(a),...,usual matrix,speye(a) ... ];

And now consider you have about 300 iterations for 10 - 15 equations like this. Your code, for just for avoiding a missing optimization, will rise exponentially.

I understand your point, that it can be done by manual hand, but you have to know that your approach is a "premature optimization" and its the root of all evil :)

Anyway, thanks for your answers.

[Jan]

Of course using sparse arrays is common. When I started programming in the early 80th, sparse data types have been even more important than today, because 64MB RAM have been very expensive.

I understand your point of premature optimization. Exactly the same argument matchs, when the handling of special values is included in the libraries of a basic data type. "Sparse" matrices are defined such, that their values are not stored sequentially, but only the non-zero values are stored together with a list of their indices. Incorporating specific operations for rarely used exceptional values would mean to embed high-level functionality in the low-level data type.

Exploiting the block-structure of an array is definitely a high-level function and it is usual to catch such exceptions manually. My code example would be very ugly, when implemented directly. But when encapsulated in a user-defined class as "SmartSparse", which has e.g. the fields "isEye", "isNull" and "data", and these flags are checked in dedicated inv(), times(), mtimes(), rdiv(), ldiv() etc., you can use your "km = ..." line without changes. In addtion you can compare the timings and results of your function just by defining B as standard sparse double array, or as your SmartSparse class.

I'd call the including this feature in the basic datatype "premature", because it matters your 300*15*8 mutliplications and divisions, while the billions of multiplications of users with standard cases would waste time for the additional tiny checks.

Adjusted linear algebra solvers, which consider the specific block-sparse structure of the system, are common in scientific software for e.g. the simulation of multi-body-systems. E.g. "natural coordinates" (Garcia de Jalon) lead to large sparse mass matrices, with a 1 in the half of the diagonal elements. Then the "isEye" and "isNull" flags would be too trivial, but only hand-coded LA solvers can be efficient. Therefore I repeat: User-specific exceptions should be managed in user-specific code.

Btw, how should speye()-eps be handled? Or 2*speye()-speye()? The number of exceptions will grow exponentially also, when you really catch all possible optimizations.