This project has moved. For the latest updates, please go here.

Sparse matrix manipulation / Solver iteration count

Nov 14, 2012 at 1:14 PM

Hello,

given a sparse matrix and an index i, what would be the most efficent way to set row i and column i to zero? Is there a better way than doing

int i = someIndex;
for (int j = 0; j < matrix.ColumnCount; j++)
{
    matrix[j, i] = 0;
    matrix[i, j] = 0;
}
And is it possible to get the iteration count when using iterative solvers?

Regards,

Chris

Nov 17, 2012 at 10:53 AM

So I tried to be clever and did

 

var enumerator = sparse.IndexedEnumerator();
int i = someIndex;
foreach (var entry in enumerator)
{
    // Set any entry in the ith row to zero.
    if (entry.Item1 == i)
    {
        sparse[i, entry.Item2] = 0;
    }

    // Set any entry in the ith column to zero.
    if (entry.Item2 == i)
    {
        sparse[entry.Item1, i] = 0;
    }
}

 

but despite the fact that it isn't faster, the results differ from the first method (maybe I just did some stupid mistake, but I don't see). Here's the code:

private static void TestSetZero(Matrix<double> matrix)
{
    int n = matrix.RowCount;
    int step = 2; // Set every n-th row/col to zero

    List<int> indices = new List<int>(n / step);
    for (int i = 0; i < n; i += step)
    {
        indices.Add(i);
    }

    SparseMatrix sparse1 = matrix.Clone() as SparseMatrix;
    SparseMatrix sparse2 = matrix.Clone() as SparseMatrix;

    if (sparse1 != null && sparse2 != null)
    {
        Console.WriteLine("Input matrix has " + sparse1.NonZerosCount + " non-zeros.");
        Console.WriteLine("Setting every " + step + "-th row/column to zero.");

        // Begin test full
        Stopwatch sw = new Stopwatch();
        sw.Start();

        foreach (var i in indices)
        {
            for (int j = 0; j < n; j++)
            {
                sparse1[j, i] = 0;
                sparse1[i, j] = 0;
            }
        }

        sw.Stop();
        Console.WriteLine("Timing full: " + sw.ElapsedMilliseconds + "");
        Console.WriteLine("Matrix has " + sparse1.NonZerosCount + " non-zeros.");

        // Begin test sparse
        sw.Restart();

        var enumerator = sparse2.IndexedEnumerator();

        foreach (var i in indices)
        {
            foreach (var entry in enumerator)
            {
                if (entry.Item1 == i)
                {
                    sparse2[i, entry.Item2] = 0;
                }

                if (entry.Item2 == i)
                {
                    sparse2[entry.Item1, i] = 0;
                }
            }
        }

        sw.Stop();
        Console.WriteLine("Timing sparse: " + sw.ElapsedMilliseconds + "");
        Console.WriteLine("Matrix has " + sparse2.NonZerosCount + " non-zeros.");

        Console.WriteLine("Matrices equal? " + sparse1.Equals(sparse2) + "");
        //Console.WriteLine("Matrices deep equal? " + DeepEqual(sparse1, sparse2) + "");

        //DumpMatrix(sparse1, "sparse1.txt");
        //DumpMatrix(sparse2, "sparse2.txt");
    }
    else
    {
        Console.WriteLine("Please load a sparse matrix.");
    }
}
I used the sparse-small.mat file for testing.

 

Nov 17, 2012 at 10:58 AM
Edited Nov 17, 2012 at 11:03 AM

BTW, I think

SparseMatrix sparse = (SparseMatrix)matrix;

should not throw an exception for a dense matrix, but instead return null. In the above code I use

SparseMatrix sparse = matrix as SparseMatrix;

which works fine.

EDIT: Though, of course, an excption while casting is more instructive than a null reference exception at a later point... don't know...

Nov 17, 2012 at 11:33 AM
Edited Nov 17, 2012 at 11:35 AM

I think the problem is that 

if (entry.Item1 == i){sparse2[i, entry.Item2] = 0;}
if (entry.Item2 == i){sparse2[entry.Item1, i] = 0;}

modify the data the enumerator is created over (i.e. the column indices and row pointers). Either MN should throw an exception here or the enumerator should be rewritten to handle modifying the internal storage while iterating over it (not sure that is possible).

Regards
Marcus 

Coordinator
Nov 17, 2012 at 11:46 AM

Yes, you cannot modify any enumerable while enumerating over it, matrices are no difference.

To set a submatrix to zero, simply use matrix.Storage.Clear(rowIndex,rowCount,columnIndex,columnCount)

Btw, you can always directly access the CRS data from a spare matrix by casting it's .Storage to SparseCompressedRowMatrixStorage.

Thanks,
Christoph

Nov 17, 2012 at 12:25 PM

Right. This works, but it's slow:

List<Tuple<int, int>> toZero = new List<Tuple<int, int>>(sparse.NonZerosCount);

var enumerator = sparse.IndexedEnumerator();

foreach (var i in indices)
{
    foreach (var entry in enumerator)
    {
        if (entry.Item1 == i)
        {
            toZero.Add(new Tuple<int, int>(i, entry.Item2));
        }
        else if (entry.Item2 == i)
        {
            toZero.Add(new Tuple<int, int>(entry.Item1, i));
        }
    }
}

foreach (var item in toZero)
{
    sparse[item.Item1, item.Item2] = 0;
}

So how do I use matrix.Storage.Clear? Tried the following, which doesn't work:

foreach (var i in indices)
{
    // Set any entry in the ith row to zero.
    sparse.Storage.Clear(i, 1, 0, n);

    // Set any entry in the ith column to zero.
    sparse.Storage.Clear(0, n, i, 1);
}

Coordinator
Nov 17, 2012 at 12:41 PM

Actually, scratch Matrix.Storage.Clear on sparse matrices until the next release, I've noticed an issue that was missed by the unit tests.

Coordinator
Nov 24, 2012 at 9:45 AM

FYI:

The issue I mentioned above has been fixed in the meantime. The next release of Math.NET Numerics will also add the following members to Matrix:

  • ClearColumn(index)
  • ClearRow(index)
  • ClearSubMatrix(rowIndex,rowCount,columnIndex,columnCount)

Thanks,
Christoph

Nov 24, 2012 at 11:37 AM

Thanks for your effort, this sounds great. I checked out the latest stuff, and there still seems to be a problem. Using direct indexing and the new clear methods give different results:

// List<int> indices = list of rows/columns to clear

foreach (var i in indices)
{
    for (int j = 0; j < n; j++)
    {
        sparse1[j, i] = 0;
        sparse1[i, j] = 0;
    }
}

foreach (var i in indices)
{
    sparse2.ClearRow(i);
    sparse2.ClearColumn(i);
}

Console.WriteLine("Equal? " + sparse1.Equals(sparse2)); // Nope, no luck!
Again using the sparse-small.mat file for testing.

Coordinator
Nov 24, 2012 at 1:29 PM

Indeed, thanks for testing. I've converted the case to a unit test and took the chance to finally implement sparse sub-matrix clearing properly. As a side effect it is now also more efficient, especially for clearing rows.

Thanks,
Christoph