Basic Support for Sparse Matrices

by Pavel Holoborodko on July 31, 2013

Starting from version 3.5.0 toolbox has support for sparse matrices with multiple precision elements.

Please note, this is the first beta version with limited functionality – at the moment toolbox can only create and display multiprecision sparse matrices:

>> S = zeros(7,3);
>> S(2,1) = eps;
>> S(5,1) = sqrt(2);
>> S(3,2) = pi;
>> S(2,3) = exp(1);
>> S(5,3) = log(2);
>> S(6,3) = sqrt(-1);
>> S = sparse(S)
 
S =
 
   (2,1)          2.22044604925031e-16                         
   (5,1)               1.4142135623731                         
   (3,2)              3.14159265358979                         
   (2,3)              2.71828182845905                         
   (5,3)             0.693147180559945                         
   (6,3)                             0 +                     1i
 
>> A = mp(S)    % Convert S to multiprecision sparse matrix of 34 digits (default)
                % Note accuracy of all constants are extended automatically.
 
A =
 
   (2,1)    1.925929944387235853055977942584927e-34
   (5,1)    1.414213562373095048801688724209698
   (3,2)    3.141592653589793238462643383279503
   (2,3)    2.718281828459045235360287471352662
   (5,3)    0.693147180559945309417232121458177
   (6,3)    0 +     1i

From now on we will extend sparse matrices support to the point of full transparency & flexibility as we have done for dense matrices.

Next steps are: element-wise arithmetic and mathematical functions, basic matrix-matrix and matrix-vector operations. Major goal after that will be to add iterative Krylov subspace methods.

Although such (currently) basic support may seem to be trivial, its implementation is very far from elementary. In fact, sparse matrices support requires by order of magnitude more efforts than dense matrices. Below we outline the major reasons why it is so.

Implementation details

The main obstacle is that MATLAB has reduced object oriented programming capabilities. Besides strange restrictions on overloading destructor and pure assignment operator (see Undocumented MEX API in MATLAB), MATLAB also does not allow elements of sparse matrices to be of user-defined type (e.g. objects of mp type).

This imposes major restrictions for customization of sparse matrices in MATLAB.

In matrix algebra system there are two kinds of operations related to matrices / n-dimensional arrays: (a) mathematics (+,-,*, solvers, etc.) and (b) manipulation (indexing, referencing, concatenation, reflection, permutation, etc.).

Both kinds require significant efforts and expertise for efficient implementation, especially taking into account all the special cases and flexibility of indexing operations in MATLAB.

Operations from category (b) do not depend on elements itself (only their size), whereas (a) is very much dependent on nature of matrix elements (precision, etc).

In case of dense matrices (where elements are allowed to be of custom type in MATLAB), we were able to avoid implementation of (b)-type routines by borrowing them from MATLAB directly. Thus we could focus on polishing numerical methods.

However when it comes to sparse matrices (where MATLAB doesn’t allow elements to be of custom type) we have to introduce absolutely new data structure (unknown and unsupported by MATLAB) and re-write complete set of functions from (a) and (b) category as a consequence.

Although it might not be noticeable for a user, sparse matrices support pushed us to re-do good part of toolbox and add (b) category routines not only for sparse matrices but even for dense matrices / n-dim arrays too.

The bottom line

After countless hours of fighting with MATLAB’s shortcomings (undocumented API, limited OOP, etc) we found it is much easier to just re-create some inaccessible / restricted functionality from the scratch.

As a result, now toolbox re-implements major portion of core of MATLAB (efficient manipulation of arrays, matrices and mathematics with them) in multiple precision. Toolbox has became completely self-contained and capable of operating independently from MATLAB, all is needed – decent high-level language and interpreter.

Another important question is the values & cost of such multiple precision MATLAB-like engine.

{ 0 comments… add one now }

Leave a Comment

Previous post:

Next post: