Creating & Using Sparse Matrix & Vectors in MATLAB - Concepts & Examples

So the first question, What is a sparse Matrix /Vector or Signal & why should we bother about this?

The answer is, a sparse matrix is any matrix that contains very few non-zero elements or conversely  a matrix whose most elements are zero. In real world situations, mostly the signal are sparse in nature, so there is a need of them in our test scenarios in order to efficiently develop algorithms & methods to deal with them in real world.

For example, the below 7x7 is sparse, with a sparsity nearly 30% (i.e., 3 out of 10 elements are non zero)

0     2     0     0     0     7     0
6     0     0     5     0     7     0
0     0     0     0     0     0     0
0     0     0     9     5     0     0
5     0     0     9     1     0     0
2     0     0     0     0     0     2
4     0     0     9     0     9     0

Another example would be a 10x1 column vector, with a sparsity of 30% (i.e., 3 out of 10 elements are non zero)
     3
     0
     0
     5
     0
     1
     0
     0
     0
     0

Sparse signal in MATLAB is designated as a special data type. 

Main reason for that:
  • The data type stores the sparse signal such that it knows what was the actual dimension of it, & wherever the matrix elements are non zero it will conserve its value & the corresponding coordinates. Thus saving a lot of memory space.
  • Vector and matrix operations can be sped up significantly by not processing the elements that are zero. Thus saving a lot of CPU time.

So, how MATLAB's sparse datatype works?

Suppose we have created a 10x1 sparse vector named 'A'  with a sparsity of 30% the MATLAB Sparse datatype will interpret it as,

A <10x1 sparse double>  %Matlab's sparse data type stores the original dimension of sparse vector
val =
   (1,1)       0.9427
            % Only the coordinates & the corresponding values in them is retained
   (8,1)       0.4177             % Where the value is non zero
  (10,1)      0.9831

Taking another example of a 7x7 sparse matrix, MATLAB's sparse datatype will store it as,
A <7x7 sparse double> %Matlab's sparse data type stores the original dimension of sparse matrix
val =
   (1,1)       0.2518
   (1,2)       0.2904
   (4,2)       0.7302
   (7,2)       0.2607
   (2,3)       0.6171
   (3,3)       0.9827
   (5,3)       0.5841
   (7,3)       0.5944
   (2,4)       0.2653
   (5,4)       0.1078
   (2,5)       0.8244
   (4,5)       0.3439
   (5,5)       0.9063
   (5,6)       0.8797
   (5,7)       0.8178
NOTE: Just like a normal matrix, on a Sparse matrix you can apply arithmetic operation, moreover all the matrix in that operation need not to be a sparse one.

Example:
Suppose Sparse MATRIX A is having the value, A is 7x7 Sparse MATRIX

A =
   (1,1)       0.2518
   (1,2)       0.2904
   (4,2)       0.7302
   (7,2)       0.2607
   (2,3)       0.6171
   (3,3)       0.9827
   (5,3)       0.5841
   (7,3)       0.5944
   (2,4)       0.2653
   (5,4)       0.1078
   (2,5)       0.8244
   (4,5)       0.3439
   (5,5)       0.9063
   (5,6)       0.8797
   (5,7)       0.8178

Suppose Sparse MATRIX B random 7x7 matrix,
B=
   -0.5583    0.6076   -0.4470    0.1269   -0.3086    0.1093    0.6001
   -0.3114   -0.1178    0.1097   -0.6568   -1.0966    1.8140    0.5939
   -0.5700    0.6992    1.1287   -1.4814   -0.4930    0.3120   -2.1860
   -1.0257    0.2696   -0.2900    0.1555   -0.1807    1.8045   -1.3270
   -0.9087    0.4943    1.2616    0.8186    0.0458   -0.7231   -1.4410
   -0.2099   -1.4831    0.4754   -0.2926   -0.0638    0.5265    0.4018
   -1.6989   -1.0203    1.1741   -0.5408    0.6113   -0.2603    1.4702

Then the arithmetic operation can be easily performed by typing C=A+B(in this case no conversion is ever needed, as matlab interpret the sparse data type at root level a matrix)

C=
   -0.3065    0.8980   -0.4470    0.1269   -0.3086    0.1093    0.6001
   -0.3114   -0.1178    0.7267   -0.3915   -0.2722    1.8140    0.5939
   -0.5700    0.6992    2.1114   -1.4814   -0.4930    0.3120   -2.1860
   -1.0257    0.9999   -0.2900    0.1555    0.1631    1.8045   -1.3270
   -0.9087    0.4943    1.8456    0.9263    0.9521    0.1565   -0.6233
   -0.2099   -1.4831    0.4754   -0.2926   -0.0638    0.5265    0.4018
   -1.6989   -0.7595    1.7685   -0.5408    0.6113   -0.2603    1.4702

Similarly it is applicable to Matrix operations like subtraction (if dimensions are agreed), multiplication(if dimensions are agreed), transpose, inverse(if non singular).


So, after understanding all these stuffs, the question will be How to generate A sparse Matrix or Vector in MATLAB?

Well its very simple, MATLAB have some inbuilt functions to achieve this task! (Yes there are many other & efficient ways to generate rather using this inbuilt commands that we will discuss in our upcoming articles.)

Function 1: speye(N) --> This command will create a sparse NxN identity matrix & store it as a sparse data type.
e.g.,
speye(6) will give the result as:
ans =

   (1,1)        1
   (2,2)        1
   (3,3)        1
   (4,4)        1
   (5,5)        1
   (6,6)        1


Function 2: sprand(M,N,% sparsity) -->This command will generate a sparse MxN random sparse matrix whose element's value varies from 0 to 1 & store it as a sparse data type. %sparsity has to be between 0 to 1.
For making a sparse column vector, set M=1 & for sparse row vector, set N=1.
e.g.,
sprand(6,6,0.4) will give the result as: (40% sparsity as 0.4*100=40
ans =
   (2,1)       0.4162
   (2,2)       0.8419
   (4,3)       0.6135
   (5,3)       0.5407
   (1,4)       0.6620
   (2,4)       0.8329
   (4,4)       0.5822
   (3,5)       0.2564
   (5,5)       0.8699
   (6,6)       0.2648


Function 3: sprandn(M,N,% sparsity) -->This command will generate a normally distributed sparse MxN random matrix whose element's value can be negative. %sparsity has to be between 0 to 1.
For making a sparse column vector, set M=1 & for sparse row vector, set N=1.
e.g.,
sprandn(6,6,0.4) will give the result as: (40% sparsity as 0.4*100=40
ans =
   (1,1)       0.0668
   (2,1)       0.0355
   (6,1)       0.3165
   (4,2)      -0.5073
   (3,3)      -0.0692
   (4,3)       0.2358
   (6,3)      -1.3429
   (4,4)       0.2458
   (2,5)       2.2272
   (4,5)       0.0700
   (4,6)      -0.6086
   (5,6)      -1.2226



Interconversion functions (From Sparse to Matrix representation & vice versa): 

If the user is having trouble understanding the sparse datatype or require a full matrix of the sparse dignal thus generated the function he/she has to use is, full( A )
* Here A  is a Sparse matrix.
Taking previous example, of sprandn(6,6,0.4),
if we want it to see it as full matrix, we would have used it as,
A=sprandn(6,6,0.4);
X=full(A);
or simply "X=full(sprandn(6,6,0.4))"  the result of which is,
X =

         0   -0.4320         0         0         0             0
         0         0             0         0    0.6489         0
         0         0      -0.3601      0         0       0.7059
         0         0             0         0         0             0
         0    1.4158         0         0   -1.6045    1.0289
         0         0             0         0         0             0


For getting the it back to sparse datatype, we need to use another function "sparse(X)", here X is a matrix which is sparse in nature (e.g, above matrix X).

So type the command "Y=sparse(X)", you will get the result as,
Y =

   (1,2)      -0.4320
   (5,2)       1.4158
   (3,3)      -0.3601
   (2,5)       0.6489
   (5,5)      -1.6045
   (3,6)       0.7059
   (5,6)       1.0289



NOTE:(Special case) If you are converting Zero Matrix or Zero Vector to a sparse datatype variable it will look like as, All zero sparse: M-by-N, example of which,

X=zeros(6,6)
X=
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0
     0     0     0     0     0     0
Y=sparse(X);
Y =
   All zero sparse: 6-by-6



0 comments:

Matlab Code & Implementation of Jacobi Method or Jacobi Iterative Method For Solving the System Of Linear Equations

In linear algebra, the Jacobi method (or Jacobi iterative method) is an algorithm for determining the solutions of a diagonally dominant system of linear equations.
The process is iterated until it converges to the solution. This algorithm is a stripped-down version of the Jacobi transformation method of matrix diagonalization. The method is named after Carl Gustav Jacob Jacobi. (Source: Wikipedia)

Suppose the system of Linear equation is denoted as 'Ax=b'
In order to solve it we need to find the solution of 'x'
To apply Jacobi Iteration we first need to check for the diagonal dominance of the measurement matrix 'A'. If the matrix 'A' is not diagonally dominant then the solution will not converge.
Example of Diagonal Dominant Matrix is:
A= 5 3 2 
      1 6 2
      2 1 7
Although if the given measurement matrix is not diagonally dominant we need to perform elementary row transform to make it diagonally domimant.
An example of matrix which is not diagonal dominant is:
A= 2 7 3
      4 1 -1
      1 -3 12
In order to make it diagonally dominant we need to exchange 1st & 2nd row, thus the matrix will be:
A= 4 1 -1
      2 7  3
      1 -3 12
Its simple to check for the diagonal dominance of a matrix, please refer this link. Check for the diagonaly Dominance Of a Matrix.
NOTE: We need to apply the same elementary transformation in the vector 'b'

NOTE: We have implemented the Jacobi Iteration method for  a 3x3 measurement matrix, the code will be applicable for higher order system also, but the check for the diagonal dormancy will become tedious.

In normal method we solve for 'x' as x= inv(A)*y (inv(A) is the inverse of matrix A)
But in Jacobi Iterative Method for 'Ax=y' what we do is as Follow:
We define matrix 'A' as the sum of two matrix 'D' & 'E'. Here 'D' is the diagonal matrix & 'E' is the remainder matrix of same order after subtracting 'D' from 'A'.
Then Ax=b will become (D+E)x=y 
=> Dx= -Ex+y
=> x= -inv(D)*E*x+inv(D)*y
*we start with a random value of x (say in MATLAB term it will be rand(3,1) for a 3x3 matrix A)
Then we apply it iteratively for K+1th step the formula, x(K+1)= -inv(D)*E*x(K)+inv(D)*y
The above method is known as Jacobi iteration for solving the system of Linear Equation.

MATLAB Code for Jacobi Method or Jacobi Iterative Method, ( Code in BOLD Text)

function jacobi(A, b, N) 
%function Jacobi with 3 Params A, b from Ax=b & N is the number of iterations
test=all((2*abs(diag(A)))- sum(abs(A),2)>=0);  %1st check of Diagonal Dominance
if test==0
    A([1 2],:) = A([2 1],:);
  
%Elementary Row Transformation of A
    b([1 2]) = b([2 1]);     %Elementary Row Transformation of b
end

test=all((2*abs(diag(A)))- sum(abs(A),2)>=0);
 
%2nd check of Diagonal Dominance after 1st exchange

if test==0
    A([2 1],:) = A([1 2],:); %if test=0 then switch back to original A
    b([2 1]) = b([1 2]);
    A([1 3],:) = A([3 1],:); %Elementary Row Transformation of A   
    b([1 3]) = b([3 1]); %Elementary Row Transformation of b
end
d=diag(A);
%Extracting the diagonal Element of A
D=diag(d); %Making the diagonal matrix D from extracted diagonal elements
D_inv=1/D; %inverse of a diagonal matrix is simply 1/D
E=A-D; %Calculating Remainder Matrix
x=rand(3,1); %Initializing the vector 'x' with a random initial value
T=-D_inv*E;
C=D_inv*b;

for j=1:N
%performing Jacobi Iterations for 'N' times
    x=T*x+C;
end
disp(x)
%Display of result

NOTE: The above code is developed for a 3x3 matrix A & 3x1 vector b, but applicable to any dimension matrix, but in order this method to converge you need to provide A such that it is diagonally dominant.

NOTE: We can also save the content of x to see what happens in each step, just by making the statement x=T*x+C; as x(j+1)=T*x(j)+C;

Test Cases:
Case 1: (When Matrix is diagonally Dominant)
>> A=[5 3 2;1 6 2;2 1 7]

A =

     5     3     2
     1     6     2
     2     1     7

>> b=[14;13;24]

b =

    14
    13
    24

>> jocobi(A,b,40)

x =

    1.0000
    1.0000
    3.0000



Case 2: (When Matrix is not diagonally Dominant, thus not converging)
 A=[5 1 5;6 2 3;7 8 1]

A =

     5     1     5
     6     2     3
     7     8     1

>> b=[11;11;18]

b =

    11
    11
    18

>> jocobi(A,b,40)

%this result is obtained only when the elementary row transformation part was not done to make The matrix A diagonally dominant.

x =

   1.0e+22 *

   -0.3795
   -1.1291
   -0.4246


Although by doing simple mldivide we will get the result
x=b\A

x =   0.4364    0.3127    0.1873 






0 comments: