Sparse Matrices
In economic applications it's not uncommon to have matrices with a large number of zero-valued elements and, because C/C++ stores zeros in the same way it stores any other numeric value, these elements can use memory space unnecessarily and can sometimes require extra computing time. Examples of applications that use sparse matrices: graph stored as adjacency matrices, optimization problems using linear algebra and sparse linear equation systems.
Sparse matrices provide a way to store data that has a large percentage of zero elements more efficiently. While full matrices internally store every element in memory regardless of value, sparse matrices store only the nonzero elements and their position. Using sparse matrices can significantly reduce the amount of memory required for data storage.
There are two main advantages of using sparse matrices:
· greatly reduced memory requirements
·
faster computations by eliminating operations on
zero elements
Because sparse matrices store both element values and position for non-zero elements, they require more storage space if number of such elements isn’t much smaller than zero elements. Another drawback of using a sparse matrix is that it doesn’t provide direct access to elements. Basic operations are also more difficult to implement compared to normal matrices.
These tradeoffs must be well considered when deciding whether to use a normal or a sparse matrix. Usually sparse matrices are used when dealing with large volumes of data that contain between 0.15% and 3% non-zero elements.
When storing a sparse matrix we need to consider two kinds of information:
There are many ways to store a sparse matrix in memory:
a) Using three arrays for line indices, column indices and values. The first two arrays can be replaced by a single array containing the index of the elements in the linearized matrix.
b) By grouping the indices and value triplets in a structure and storing them into an array or a linked list.
c) Using a bit array for marking the positioning of the non-zero elements in the linearized matrix and a separate array for values.
In the following sections we’ll use the first form (we’ll use three dynamically allocated arrays for storing line indices, column indices and values). For easier manipulation all the information related to a matrix is grouped into a C structure:
struct MatriceRara
{
//
vectorii pentru stocarea liniei/coloanei si valorii
int
*Linie, *Coloana, *Valoare;
//
numarul de elemente nenule din matrice
int
nrElemente, nrLinii, nrColoane;
};
The conversions between normal matrices and sparse matrices and vice versa are straightforward:
// Conversie matrice normala -> matrice rara
MatriceRara Conversie(int **mat, int m, int n)
{
MatriceRara
matR;
//
determinarea numarului de elemente
int
i, j, nrElemente = 0, indiceMatR = 0;
for
(i = 0; i < m; i++)
for
(j = 0; j < n; j++)
if
(mat[i][j] != 0)
nrElemente++;
//
initializare matrice rara
matR.Linie
= new int[nrElemente];
matR.Coloana
= new int[nrElemente];
matR.Valoare
= new int[nrElemente];
matR.nrElemente
= nrElemente;
matR.nrLinii
= m;
matR.nrColoane
= n;
//
copierea valorilor nenule
for
(i = 0; i < m; i++)
for
(j = 0; j < n; j++)
if
(mat[i][j] != 0)
{
matR.Linie[indiceMatR]
= i;
matR.Coloana[indiceMatR]
= j;
matR.Valoare[indiceMatR]
= mat[i][j];
indiceMatR++;
}
return
matR;
}
// Conversie matrice rara -> matrice normala
int ** Conversie(MatriceRara& matR)
{
//
declarare si alocare matrice
int
**mat;
mat =
new int*[matR.nrLinii];
for
(int i = 0; i < matR.nrLinii; i++)
{
mat[i]
= new int[matR.nrColoane];
for
(int j = 0; j < matR.nrColoane; j++)
mat[i][j]
= 0;
}
//
copiere elemente nenule
for
(int i = 0; i < matR.nrElemente; i++)
mat[matR.Linie[i]][matR.Coloana[i]]
= matR.Valoare[i];
//
returnare rezultat
return
mat;
}
One way to implement matrix operations on sparse matrices is to simulate direct access to elements using Get/Set functions and use the traditional algorithms. This is the simplest method but has the drawback that is very slow because we must consider all elements (even non-zero ones) and use the sequential access functions that need O(nrZ) time to complete.
A better solution is to develop specialized algorithms that use the particularities of the sparse matrices to provide a fast execution time. In order to speed up the operations we’ll assume that elements inside the sparse matrix are sorted according to (line, column); all operations presented here preserve this property on the matrices.
The algorithm to compute the transpose of a sparse matrix is very simple. All we need to do is swap the columns and line indices and sort the resulted arrays in order to preserve the mentioned property:
// Calculeaza transpusa unei matrici stocate ca
matrice rara
MatriceRara Transpusa(const MatriceRara& mat)
{
//
initializare rezultat
MatriceRara
R;
R.Linie
= new int[mat.nrElemente];
R.Coloana
= new int[mat.nrElemente];
R.Valoare
= new int[mat.nrElemente];
R.nrElemente
= mat.nrElemente;
R.nrLinii
= mat.nrColoane;
R.nrColoane
= mat.nrLinii;
//
copiere date si inversare linie-coloana
for
(int i = 0; i < mat.nrElemente; i++)
{
R.Linie[i]
= mat.Coloana[i];
R.Coloana[i]
= mat.Linie[i];
R.Valoare[i]
= mat.Valoare[i];
}
//
sortare rezultat dupa linie, coloana
int
temp;
for
(int i = 0; i < R.nrElemente; i++)
for
(int j = i + 1; j < R.nrElemente; j++)
if
(R.Linie[i] > R.Linie[j] ||
R.Linie[i] == R.Linie[j] && R.Coloana[i]
> R.Coloana[j])
{
temp = R.Linie[i]; R.Linie[i] = R.Linie[j];
R.Linie[j] = temp;
temp = R.Coloana[i]; R.Coloana[i] = R.Coloana[j];
R.Coloana[j] = temp;
}
return
R;
}
The algorithm for sparse matrix addition has three steps:
Because the matrices are sorted by (line, column) we can do only one pass to determine the number of elements and to compute the values. In the process we take care to eliminate zero elements in the result.
// Aduna doua matrici rare
MatriceRara Adunare(const MatriceRara& M1,
const MatriceRara& M2)
{
MatriceRara
R = {NULL, NULL, NULL, 0, 0, 0};
if
(M1.nrLinii != M2.nrLinii && M1.nrColoane != M2.nrColoane)
{
cout << "Eroare: Matricele trebuie sa
aiba dimensiuni egale." << endl;
return
R;
}
//
determinarea numarului de elemente al matricei rezultat
int i
= 0, j = 0, nr = 0;
while
(i < M1.nrElemente && j < M2.nrElemente)
{
if
(M1.Linie[i] < M2.Linie[j])
i++;
else
if (M1.Linie[i] > M2.Linie[j])
j++;
else
{
if
(M1.Coloana[i] < M2.Coloana[j])
i++;
else
if (M1.Coloana[i] > M2.Coloana[j])
j++;
else
{
if
(M1.Valoare[i] + M2.Valoare[j] == 0)
nr
= nr + 2; // element comun nul
else
nr
= nr + 1; // element comun nenul
i++;
j++;
}
}
}
// numarul de elemente al rezultatului
nr =
M1.nrElemente + M2.nrElemente - nr ;
//
initializare rezultat
R.Linie
= new int[nr];
R.Coloana
= new int[nr];
R.Valoare
= new int[nr];
R.nrElemente
= nr;
R.nrLinii
= M1.nrLinii;
R.nrColoane
= M1.nrColoane;
//
calculare rezultat
int k
= 0; // indice in matricea
rezultat
i =
0; j = 0;
while
(k < R.nrElemente)
{
if
(M1.Linie[i] < M2.Linie[j])
{
R.Linie[k]
= M2.Linie[i];
R.Coloana[k]
= M2.Coloana[i];
R.Valoare[k]
= M2.Valoare[i];
i++;
k++;
}
else
if (M1.Linie[i] > M2.Linie[j])
{
R.Linie[k]
= M1.Linie[j];
R.Coloana[k]
= M1.Coloana[j];
R.Valoare[k]
= M1.Valoare[j];
j++;
k++;
}
else
{
if
(M1.Coloana[i] < M2.Coloana[j])
{
R.Linie[k]
= M2.Linie[i];
R.Coloana[k]
= M2.Coloana[i];
R.Valoare[k]
= M2.Valoare[i];
i++;
k++;
}
else
if (M1.Coloana[i] > M2.Coloana[j])
{
R.Linie[k]
= M1.Linie[j];
R.Coloana[k]
= M1.Coloana[j];
R.Valoare[k]
= M1.Valoare[j];
j++;
k++;
}
else
{
if
(M1.Valoare[i] + M2.Valoare[j] == 0)
{
i++;
j++; // element comun nul
}
else
{
R.Linie[k]
= M1.Linie[i];;
R.Coloana[k]
= M1.Coloana[i];
R.Valoare[k] =
M1.Valoare[i] + M2.Valoare[j];
i++;
j++; k++; // element comun nenul
}
}
}
}
return
R;
}
Matrix multiplication is done in two phases. In the first phase we simulate the operation in order to determine the number of elements in the result set and in the second one we compute the result.
The operation is executed only if the matrices have compatible sizes.
In order to compute the value of one element in the result set we do the following:
· we initialize the current value with 0
· for every element in the current column from the first matrix
a. we search for a corresponding element in the second matrix
b. if we find such element we add the product to the current value
· if the current value is non-zero we store it in the result matrix
// Inmultirea a doua matrici rare
MatriceRara Inmultire(MatriceRara& M1,
MatriceRara& M2)
{
MatriceRara
R = {NULL, NULL, NULL, 0, 0, 0};
if
(M1.nrColoane != M2.nrLinii)
{
cout
<< "Eroare: Dimensiunile matricelor nu corespund." <<
endl;
return
R;
}
//
determinarea numarului de elemente al rezultatului (simulare inmultire)
int
nr = 0, suma;
for
(int i = 0; i < M1.nrLinii; i++)
for
(int j = 0; j < M2.nrLinii; j++)
{
suma
= 0;
//
calcul valoare element
for
(int col = 0; col < M1.nrColoane; col++)
for
(int i1 = 0; i1 < M1.nrElemente; i1++)
// daca exista in M1
if
(M1.Linie[i1] == i && M1.Coloana[i1] == col)
for
(int i2 = 0; i2 < M2.nrElemente; i2++)
if (M2.Linie[i2] == col && M2.Coloana[i2]
== j) // si in M2
suma
+= M1.Valoare[i1] * M2.Valoare[i2];
if
(suma) nr++;
}
//
initializare rezultat
R.Linie
= new int[nr];
R.Coloana
= new int[nr];
R.Valoare
= new int[nr];
R.nrElemente
= nr;
R.nrLinii
= M1.nrLinii;
R.nrColoane
= M2.nrColoane;
//
calcul rezultat
int
indice = 0;
for
(int i = 0; i < M1.nrLinii; i++)
for
(int j = 0; j < M2.nrLinii; j++)
{
suma
= 0;
//
calcul valoare element
for
(int col = 0; col < M1.nrColoane; col++)
for
(int i1 = 0; i1 < M1.nrElemente; i1++)
// daca exista in M1
if
(M1.Linie[i1] == i && M1.Coloana[i1] == col) for
(int i2 = 0; i2 < M2.nrElemente; i2++)
if
(M2.Linie[i2] == col && M2.Coloana[i2] == j) // si in M2
suma
+= M1.Valoare[i1] * M2.Valoare[i2];
if
(suma)
{
R.Linie[indice]
= i;
R.Coloana[indice]
= j;
R.Valoare[indice]
= suma;
indice++;
}
}
return
R;
}
1. Write a function to eliminate all elements with values less than a threshold from a sparse matrix.
2. Write functions to read a sparse matrix from the console and to write the result in both forms (sparse and normal).
3. Write a function to multiply a sparse matrix with a constant. Use this function in conjunction with the addition to simulate sparse matrix subtraction.
4. Write a function to multiply a sparse matrix by a line vector.
5. Give a solution to eliminate the necessity of the first pass in the matrix multiplication algorithm.
6. An airline company operates direct flights between n cities. Write an application to determine the minimum number of flight segments between any two cities.
Solution:
We model the problem as an undirected graph with cities represented as vertexes and flight segments as edges and we store the graph as an adjacency matrix A (using a sparse matrix). In this case Ak will contain all flight routes of length k between cities.