I've been working on an algorithm for removing non-zero elements from a large sparse (ish) matrix as a way of decreasing the computation time of a numerical solution to a large system of coupled ODEs of the form of the logistic equation.
The basic idea is that every T timesteps you assess which elements in the state matrix, B, are greater than some threshold. FOr those deemed insignificant (i.e. below the threshold) you remove (set to zero) the associated row/column in the sparse interaction matrix, A, so that the computation of AxB will require fewer flops.(Note the diagonal remains non-zero so that the associated elements of the state matrix do not blow up)
The algorithm is as follows:
Scan through B vector for elements less that thresh and record locations
Iterate through non-zero elements of A matrix corresponding to locations recorded in 1. and record locations
- Set elements of A corresponding to locations recorded in 2. to zero.
I need to separate steps 2. and 3. because if you set non-zero elements of a sparse matrix (arma::sp_mat) to zero, you destroy the iterators and therefore can't achieve the result in a single pass.
Can anyone offer any suggestions as to how to improve this algorithm?
(Note I'm not a computer scientist and so there may be an obvious improvement I'm not aware of. MMMMMany thanks for trying!!)
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
int main() {
mat A = {{1, 0.2, 0, 0.2, 0, 0, 0.2, 0, 0.2, 0},
{0, 1, 0.2, 0, 0.2, 0, 0, 0.2, 0 , 0},
{0, 0, 1, 0.2, 0, 0, 0.2, 0, 0, 0.2},
{0.2, 0, 0, 1, 0.2, 0, 0.2, 0.2, 0, 0.2},
{0, 0.2, 0, 0.2, 1, 0.2, 0, 0, 0, 0.2},
{0, 0, 0.2, 0, 0.2, 1, 0, 0.2, 0, 0},
{0.2, 0, 0, 0.2, 0.2, 0, 1, 0, 0, 0.2},
{0, 0.2, 0, 0.2, 0, 0, 0, 1, 0.2, 0},
{0, 0, 0.2, 0, 0.2, 0.2, 0, 0, 1, 0},
{0, 0, 0, 0.2, 0, 0.2, 0, 0, 0.2, 1}};
sp_mat Asp(A); // Interaction matrix in sparse representation
vec B = {1, 1, 0.1, 1, 0.1, 1, 1, 0.1, 1, 1};
double thresh = 0.2;
vec index(Asp.n_rows); // obj for storing locations of insignificant elements b
int count = 0;
for (int i=0; i<B.n_rows; i++) {
if (B(i) < thresh) {
index(count) = i;
count ++;
}
}
index.resize(count);
sp_mat tmp = Asp; // temp copy of interaction matrix
for (int i=0; i<index.n_rows; i++) {
// iterate through non-zero elements in col corresponding to each element of index
sp_mat::iterator bc = tmp.begin_col(index(i));
sp_mat::iterator ec = tmp.end_col(index(i));
// iterate through non-zero elements in row corresponding to each element of index
sp_mat::row_iterator br = tmp.begin_row(index(i));
sp_mat::row_iterator er = tmp.end_row(index(i));
// mat for storing locations of non-zero elements to be removed
mat loc(2, 2 * tmp.n_cols);
int count2 = 0;
for (auto j = bc; j != ec; j++) {
loc(0, count2) = j.row();
loc(1, count2) = j.col();
count2++;
}
for (auto j = br; j != er; j++) {
loc(0, count2) = j.row();
loc(1, count2) = j.col();
count2++;
}
loc.resize(2, count2);
for (int j = 0; j < loc.n_cols; j++) {
tmp(loc(0, j), loc(1, j)) = 0.0;
}
}
tmp.diag().ones();
cout << "Number of non-zero terms interaction matrix = " << Asp.n_nonzero << endl;
cout << "Number of non-zero terms reduced matrix = " << tmp.n_nonzero << endl;
}
Output:
Number of non-zero terms interaction matrix = 45
Number of non-zero terms reduced matrix = 24