Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

CompressedColumnStorage<T>.Transpose yields wrong result leading to sequential issues in Solvers #43

Closed
ChristoWolf opened this issue Mar 27, 2024 · 9 comments · Fixed by #44

Comments

@ChristoWolf
Copy link

Hi!

Issue

I have several issues, most of which stem from the first one here (let SparseMatrix be the one from Math.NET and CSparseMatrix the one from Csparse.NET):

  • Transposing a CSparseMatrix yields incorrect results.
    E.g. when transposing sparse matrix a
    var a = (SparseMatrix)SparseMatrix.Build.SparseOfRowArrays([[1, 2], [3, 4], [5, 6]]);
    using
    var aStorage = (SparseCompressedRowMatrixStorage<double>)a.Storage;
    var aSparse = new CSparseMatrix(a.RowCount, a.ColumnCount)
    {
          // Based on https://github.com/wo80/CSparse.NET/wiki/Math.NET-Numerics-and-CSparse.
      RowIndices = aStorage.ColumnIndices,
      ColumnPointers = aStorage.RowPointers,
      Values = aStorage.Values
    };
    var aDash = aSparse.Transpose();
    one would expected aDash to be
    1 3 5
    2 4 6
    
    in matrix notation.
    However, I instead get
    1 2 0
    3 4 0
    
    or
    %%MatrixMarket matrix coordinate real general
    2 3 4
    1 1 1
    2 1 3
    1 2 2
    2 2 4
    
    in Matrix Market format (obtain via MatrixMarketWriter.WriteMatrix).
  • aSparse is incorrectly dumped as
    %%MatrixMarket matrix coordinate real general
    3 2 4
    1 1 1
    2 1 2
    1 2 3
    2 2 4
    
    although its Values property is [1,2,3,4,5,6] (with 3 rows and 2 columns).
  • Multiplying matrices is unclear to me.
    Is supposed A = B*C equivalent to var A = B.ParallelMultiply(C) or var A = C.ParallelMultiply(B)?
    When I execute var aTilde = aDash.ParallelMultiply(aSparse); I get
    5 11
    11 25
    
    and with var aTilde = aSparse.ParallelMultiply(aDash); I get
    10 14 0
    14 20 0
    0 0 0
    
    .
    Mathematically, aDash*aSparse and aSparse*aDash should be
    10 14
    14 20
    
    and
    5 11 0
    11 25 0
    0 0 0
    
    respectively.
    Here, I used the values for aDash and aSparse that were computed by the library (i.e. the unexpected ones).
  • These issues of course propagate to or also happen in the Solvers.
    E.g. var solver = SparseQR.Create(aTilde, ColumnOrdering.Natural);
    yields incorrect factorization matrices and applying Solve (with ColumnOrdering.Natural) yields the wrong result.

My original mathematical problem

I want to numerically solve the regularized l_2-minimization problem|| Ax -b ||_2 + λ|| x ||_2, where A is some noise operator, b a noisy signal, x the denoised signal that I want to find, and λ > 0 is an arbitrary regularization parameter.
´x´ and ´b´ have the same number of entries.

This is equivalent to an l_2-minimization problem of the form || Cx_c - y ||_2, with C in |R^mxn being a sparse matrix dependent on λ with m > n.
According to the documentation of SparseQR.Solve,

If m >= n a least-squares problem (min |Ax-b|) is solved.

which is exactly what I want there.

Unfortunately, this returns wrong values for my actual problem that I want to solve, where m is about 8000 and n about 4000.
So I heavily simplified the problem, which yields the results above.

Versions

  • .NET SDK: 8.0.202
  • CSparse 3.8.1
  • MathNet.Numerics 5.0.0
  • MathNet.Numerics.MKL.Win-x64 3.0.0
@wo80
Copy link
Owner

wo80 commented Mar 27, 2024

I haven't tested any of the code you provided, but the first issue I see is

var aSparse = new CSparseMatrix(a.RowCount, a.ColumnCount)

That would only be ok for a square matrix. If a.RowCount != a.ColumnCount, you should write

var aSparse = new CSparseMatrix(a.ColumnCount, a.RowCount)

Can you please test this and report back. If the problems reside, I'll have a closer look.

@ChristoWolf
Copy link
Author

ChristoWolf commented Mar 28, 2024

Hi @wo80!

Thanks for the quick help, that was it!
But I do not understand at all why this is needed.
Also, it gives the transposed matrix, so I would need to transpose my minimization problem.
No biggy, but still would like to understand the thought process behind this limitation.

@wo80
Copy link
Owner

wo80 commented Mar 28, 2024

You need to understand how the matrix entries are stored (see for example https://www.smcm.iqfr.csic.es/docs/intel/mkl/mkl_manual/appendices/mkl_appA_SMSF.htm).

MathNet uses CSR storage (matrix entries are stored row-wise) while CSparse uses CSC storarge (entries are stored column-wise). So if you take the CSR storage and just interpret it as CSC, it will represent the transposed matrix. This means, if your original matrix stored as CSR is m x n, then it will represent a n x m matrix in CSC.

Also, it gives the transposed matrix, so I would need to transpose my minimization problem. No biggy, but still would like to understand the thought process behind this limitation.

Yes, if you want to solve a least-squares problem, make sure that the CSparse m x n matrix has m >= n, so the MathNet matrix should be built as n x m. You can also feed matrices with m < n to SparseQR, but it will then compute the transposed matrix internally:

else
{
// Ax=b is underdetermined
var AT = A.Transpose();
var p = AMD.Generate(AT, order);
// Ordering and symbolic analysis
C.SymbolicAnalysis(AT, p, order == ColumnOrdering.Natural);
// Numeric QR factorization of A'
C.Factorize(AT, progress);
}

@wo80
Copy link
Owner

wo80 commented Mar 28, 2024

I just looked up the documentation of the Multiply method and I agree it's a bit sparse:

/// <summary>
/// Sparse matrix multiplication, C = A*B
/// </summary>
/// <param name="other">The sparse matrix multiplied to this instance.</param>
/// <returns>C = A*B</returns>
public CompressedColumnStorage<T> Multiply(CompressedColumnStorage<T> other)

Do you think changing the summary would be enough:

/// <summary>
/// Sparse matrix multiplication, C = A*B, where A is current instance.
/// </summary>

@ChristoWolf
Copy link
Author

Thanks a ton for clarifying, seems that I really misunderstood a few things.

I am actually thinking about rewriting most of my code to use CSparse directly, the performance is incredible!

Do you think changing the summary would be enough:

Absolutely, that makes it very explicit.
The issue that I have with

The sparse matrix multiplied to this instance.

is that one could interpret it as left multiplication, as right multiplication is quite uncommon.

@epsi1on
Copy link
Collaborator

epsi1on commented Mar 29, 2024

I am actually thinking about rewriting most of my code to use CSparse directly, the performance is incredible!

I think you are right. The CSparse.NET is a port of CSparse by Tim. Davis. If i remember right, CSparse.NET had reasonably good performance compared to original CSparse although the original one is in native C or C++.

@wo80
Copy link
Owner

wo80 commented Mar 29, 2024

@ChristoWolf

I am actually thinking about rewriting most of my code to use CSparse directly, the performance is incredible!

In case you were using the dense QR from MathNet (or even native MKL), the difference will be considerable. You can set up your matrix using the CSparse.NET CoordinateStorage class, so in case you aren't using any other parts from MathNet, switching should be easy.

The issue that I have with [...] is that one could interpret it as left multiplication, as right multiplication is quite uncommon.

Alright, I will make sure to add a note that the matrix A is multiplied by the other matrix from the right.

@wo80
Copy link
Owner

wo80 commented Mar 29, 2024

Please leave this issue open! I will close it with the next pull request.

@ChristoWolf
Copy link
Author

Thanks a ton, @wo80!

@ChristoWolf ChristoWolf changed the title CompressedColumnStorage<T>.Transpose yields wrong result leading to sequential issues in Solvers CompressedColumnStorage<T>.Transpose yields wrong result leading to sequential issues in Solvers Apr 4, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants