Fast transpose of (sparse) matrices C++

Fast transposed sparse matrix algorithm

In many applications, especially deep learning, we often encounter matrices in which most of the elements are 0, which are called sparse matrices. Using traditional matrix storage methods wastes a lot of storage space. In order to solve this problem, special storage methods are usually used in data structures to store sparse matrices.

Pre-environment

WSL2 Ubuntu-20.04
gcc-11.4
C++20 standard

1. Data structure: triplet

First, in order to represent non-zero elements in a sparse matrix, we introduce a structure called Triplet. This structure contains three elements: row index, column index and value (the default type here is int). Note that although it is said to be a matrix, Triplet is actually stored in the actual storage space according to the sequential structure of the array, and in row-major order, that is, in a continuous storage space middle.

struct Triplet
{<!-- -->
    int row, col, value;
Triplet( int r, int c, int v ) : row( r ), col( c ), value( v ) {<!-- -->}
    Triplet() = default;
};

2. Read sparse matrix

In general programming languages, the subscripts ofarrays start from0while the matrix subscripts in mathematics start from1for adaptation. The matrix in mathematics uses the subscript here and starts with 1. In order not to waste the first storage space, the first storage space is used to save the rows, columns and the number of non-zero elements of the matrix (useful for subsequent algorithms) .
First, the data of the sparse matrix needs to be read from the standard input. The size of the matrix and the number of non-zero elements are read first, followed by the row index, column index and value of each non-zero element.
Test input: The first row has three integers, representing the row n, column m and the number of non-zero elements k of the matrix, followed by k rows, each row is a triplet. like:
7,8,6
1,4,-5
2,2,16
2,7,1
4,2,2
4,7,10
5,5,-10

int n, m, k;
    char ch; //ignore commas
    cin >> n >> ch >> m >> ch >> k;
    vector<Triplet> matrix(k + 1);//The first element is used to save the relevant information of the array
    matrix[0].row = n;
    matrix[0].col = m;
    matrix[0].value = k;
    for (int i = 1; i <= k; i + + )
    {<!-- -->
        cin >> matrix[i].row >> ch >> matrix[i].col >> ch >> matrix[i].value;
    }

3. Fast transposition algorithm

The core idea of fast transpose is to pre-compute the transposed position instead of directly swapping rows and columns.

  • Calculate the number of non-zero elements in each column: We first build an array named num to store the number of non-zero elements in each column.
 //Construct num array, the number of non-zero elements in the previous column of each column
    vector<int> num(m + 1, 0);//Construct an array of m + 1 zeros
    for (int i = 1; i <= k; i + + )
    {<!-- -->
        num[matrix[i].col] + + ;
    }
  • Determine the position of the first non-zero element of each row after transposition: Next, we construct the cpot array, each element of which represents the first element of each row after transposition The position of the non-zero element. For cpot, there is a mathematical formula

    {

    cpot

    ?

    [

    1

    ]

    =

    1

    ;

    cpot

    ?

    [

    col

    ?

    ]

    =

    cpot

    ?

    [

    col

    ?

    ?

    1

    ]

    +

    num

    ?

    [

    col

    ?

    ?

    1

    ]

    ;

    (

    2

    c

    o

    l

    m

    )

    \left\{\begin{array}{ll}\operatorname{cpot}[1]=1;\\operatorname{cpot}[\operatorname{col}]=\operatorname{cpot}[\operatorname{col}- 1] + \operatorname{num}[\operatorname{col}-1]; & amp;\quad\mathrm{(2\leq col\leq m)}\end{array}\right.

    {cpot[1]=1;cpot[col]=cpot[col?1] + num[col?1];?(2≤col≤m)?
    cpot[1]=1 is because the first row stores the entire matrix size and the number of non-zero elements, and cannot be empty.
    When col>=2, since the actual storage of the matrix is sequential storage, the i after transposition >The first non-zero element of the row is in column i-1 before transposition (that is, after transposition All non-zero elements of row i-1), and before transposition after the first non-zero element of the previous row, hence this calculation formula.

Based on this formula, we can construct the cpot array

//Construct the cpot array, that is, the first non-zero element subscript of each row after transposition
    vector<int> cpot( m + 1, {<!-- --> 0 });
    cpot[1] = 1;
    for (int i = 2; i <= m; i + + )
    {<!-- -->
        cpot[i] = cpot[i - 1] + num[i - 1];
    }
  • Transpose: With the cpot array, we can quickly determine the position of each element after transposition. In this way, we can directly place each element into its transposed position without the need to swap rows and columns.
 //result array
    vector<Triplet> result( k + 1 );
    result[0] = matrix[0];
    for (int i = 1; i <= k; i + + )
    {<!-- -->
        //First fill in the first element position, in + +
        result[cpot[matrix[i].col] + + ] = matrix[i];
    }
    for ( auto & amp; & amp;i : result )
    {<!-- -->
        //Exchange rows and columns
        swap(i.row, i.col);
    }

This cpot matrix is not unchanged. When the position of the first non-zero element in each row is confirmed, in the sequential storage space, the non-first non-zero element of this row will be vacated between each non-zero element. position. At this time, you only need to add one to the value of cpot[col] to get the position of the second non-zero element in this row, and so on.

4. Output the transposed sparse matrix

After completing the above steps, directly output the transposed sparse matrix.

void out( const vector<Triplet> & amp;v, int pos)
{<!-- -->
    for ( int i = pos; i < v.size(); i + + )
    {<!-- -->
        cout << v[i].row << "," << v[i].col << "," << v[i].value << endl;
    }
}

5. Time complexity

Since the array is traversed twice during the construction process and once when calculating the position after transposition, the time complexity of this algorithm is linear and is only related to the number of non-zero elements k and the size of the matrix n, m irrelevant, that is

T

(

n

,

m

,

k

)

=

O

(

k

)

\quad T(n,m,k)=O(k)

T(n,m,k)=O(k)

Complete code

#include 
using namespace std;
struct Triplet
{
    int row, col, value;
    Triplet( int r, int c, int v ) : row( r ), col( c ), value( v ) {}
    Triplet() = default;
};
void out( const vector<Triplet> & amp;v, int pos)
{<!-- -->
    for ( int i = pos; i < v.size(); i + + )
    {<!-- -->
        cout << v[i].row << "," << v[i].col << "," << v[i].value << endl;
    }
}
int main()
{
    int n, m, k;
    char ch; //ignore commas
    cin >> n >> ch >> m >> ch >> k;
    vector matrix( k + 1 );
    matrix[0].row = n;
    matrix[0].col = m;
    matrix[0].value = k;
    for (int i = 1; i <= k; i + + )
    {
        cin >> matrix[i].row >> ch >> matrix[i].col >> ch >> matrix[i].value;
    }
    //Construct num array, the number of non-zero elements in the previous column of each column
    vector num( m + 1, 0);
    for (int i = 1; i <= k; i + + )
    {
        num[matrix[i].col] + + ;
    }
    //Construct the cpot array, that is, the first non-zero element subscript of each row after transposition
    vector<int> cpot( m + 1, {<!-- --> 0 });
    cpot[1] = 1;
    for (int i = 2; i <= m; i + + )
    {<!-- -->
        cpot[i] = cpot[i - 1] + num[i - 1];
    }
    //result array
    vector<Triplet> result( k + 1 );
    result[0] = matrix[0];
    for (int i = 1; i <= k; i + + )
    {<!-- -->
        //First fill in the first element position, in + +
        result[cpot[matrix[i].col] + + ] = matrix[i];
    }
    for ( auto & amp; & amp;i : result )
    {<!-- -->
        //Exchange rows and columns
        swap(i.row, i.col);
    }
    out(result, 0);
    return 0;
}

Test output:
8,7,6
2,2,16
2,4,2
4,1,-5
5,5,-10
7,2,1
7,4,10

syntaxbug.com © 2021 All Rights Reserved.