How MATLAB Stores Sparse Matrices

    To understand why the above examples are so slow, you need to understand how MATLAB stores its sparse matrices. An n-by-n MATLAB sparse matrix is stored as three arrays; I'll call them p, i, and x. These three arrays are not directly accessible from M, but they can be accessed by a mexFunction. The nonzero entries in column j are stored in i(p(j):p(j+1)-1) and x(p(j):p(j+1)-1), where x holds the numerical values and iholds the corresponding row indices. Below is a very small example. First, I create a full matrix and convert it into a sparse one. This is only so that you can easily see the matrix C and how it's stored in sparse form. You should never create a sparse matrix this way, except for tiny examples.

C = [4.5 0.0 3.2 0.03.1 2.9 0.0 0.90.0 1.7 3.0 0.03.5 0.4 0.0 1.0 ] ;C = sparse (C) 

C =   (1,1)       4.5000   (2,1)       3.1000   (4,1)       3.5000   (2,2)       2.9000   (3,2)       1.7000   (4,2)       0.4000   (1,3)       3.2000   (3,3)       3.0000   (2,4)       0.9000   (4,4)       1.0000

    Notice that the nonzero entries in C are stored in column order, with sorted row indices. The internal p, i, and xarrays can be reconstructed as follows. The find(C) statement returns a list of "triplets," where the kth triplet isi(k), j(k), and x(k). This specifies that C(i(k),j(k)) is equal to x(k). Next, find(diff(...)) constructs the column pointer array p (this only works if there are no all-zero columns in the matrix).

[i j x] = find (C) ;n = size(C,2) ;p = find (diff ([0 ; j ; n+1])) ;for col = 1:n    fprintf ('column %d:/n    k      row index     value/n', col) ;    disp ([(p(col):p(col+1)-1)' i(p(col):p(col+1)-1) x(p(col):p(col+1)-1)])end 

column 1:    k      row index     value    1.0000    1.0000    4.5000    2.0000    2.0000    3.1000    3.0000    4.0000    3.5000column 2:    k      row index     value    4.0000    2.0000    2.9000    5.0000    3.0000    1.7000    6.0000    4.0000    0.4000column 3:    k      row index     value    7.0000    1.0000    3.2000    8.0000    3.0000    3.0000column 4:    k      row index     value    9.0000    2.0000    0.9000   10.0000    4.0000    1.0000

Now consider what happens when one entry is added to C:

C(3,1) = 42 ;[i j x] = find (C) ;n = size(C,2) ;p = find (diff ([0 ; j ; n+1])) ;for col = 1:n    fprintf ('column %d:/n    k      row index     value/n', col) ;    disp ([(p(col):p(col+1)-1)' i(p(col):p(col+1)-1) x(p(col):p(col+1)-1)])end 

column 1:    k      row index     value    1.0000    1.0000    4.5000    2.0000    2.0000    3.1000    3.0000    3.0000   42.0000    4.0000    4.0000    3.5000column 2:    k      row index     value    5.0000    2.0000    2.9000    6.0000    3.0000    1.7000    7.0000    4.0000    0.4000column 3:    k      row index     value    8.0000    1.0000    3.2000    9.0000    3.0000    3.0000column 4:    k      row index     value   10.0000    2.0000    0.9000   11.0000    4.0000    1.0000

    and you can see that nearly every entry in C has been moved down by one in the i and x arrays. In general, the single statement C(3,1)=42 takes time proportional to the number of entries in matrix. Thus, looping nnz(A)times over the statement A(i,j)=A(i,j)+... takes time proportional to nnz(A)^2.

 

 

http://blogs.mathworks.com/loren/2007/03/01/creating-sparse-finite-element-matrices-in-matlab/