Add path for DenseMatrix multiplication with SparseVector #933
Conversation
|
Actually I realise this breaks behaviour around double.NaN :( Closing [Test]
public void DenseVsSparse()
{
Matrix<double> m = Matrix<double>.Build.DenseOfArray(new double[,] { { double.NaN, 0 }, { 0, 0 } });
double[] arr = new double[] { 0, 0 };
Vector<double> sparseResult = m * Vector<double>.Build.SparseOfArray(arr);
Vector<double> denseResult = m * Vector<double>.Build.DenseOfArray(arr);
for (var i = 0; i < denseResult.Count; i++)
{
Assert.AreEqual(denseResult[i], sparseResult[i], $"Failed at index {i}");
/* Message:
Failed at index 0
Expected: NaN
But was: 0.0d
*/
}
} |
|
Failures on master:
Failures on PR in addition to above:
So it is not as cut-and-dried as I thought. I am reopening and will await your decision on what to do, if anything. |
| } | ||
| else | ||
|
|
||
| if (rightSide.Storage is SparseVectorStorage<Complex> sparseRight) |
There was a problem hiding this comment.
This logic appears 3 times in this file, only the matrix indices are swapped and/or its element is conjugated. I would extract the core part in one function, and call it 3 times; each time changing only the logic which manipulates the matrix value
| } | ||
| else | ||
|
|
||
| if (rightSide.Storage is SparseVectorStorage<Complex32> sparseRight) |
There was a problem hiding this comment.
Same comment as in the Complex/DenseMatrix.cs file.
| } | ||
| else | ||
|
|
||
| if (rightSide.Storage is SparseVectorStorage<double> sparseRight) |
There was a problem hiding this comment.
Same comment as in the Complex/DenseMatrix.cs file.
| } | ||
| else | ||
|
|
||
| if (rightSide.Storage is SparseVectorStorage<float> sparseRight) |
There was a problem hiding this comment.
Same comment as in the Complex/DenseMatrix.cs file.
It fails the same 9 tests as the comment above, relating to NaN/Inf behaviour differences between Sparse and Dense multiplication (AFAIR): those which exist already on master (exposed by 8d6317d) and those which this PR introduces more of. There needs to be a decision as to what to do (if anything) about the behavioural differences before the PR goes anywhere. It is easy enough to use the fast logic in consuming code (that's what I ended up doing). |
|
I don't really see where the master branch fails, can you please share a link? |
Currently when multiplying a
DenseMatrixwith aVectorthere is a path to aILinearAlgebraProviderwhen the multiplicand and the result are bothDenseVector, otherwise we fall back to a base implementation:mathnet-numerics/src/Numerics/LinearAlgebra/Double/Matrix.cs
Lines 158 to 169 in 4b13ff4
When
rightSideis aSparseVectorthenrightSide[j]is a binary search on the sparse indices. We can instead just iterate the sparse indices as the column indices, similar to theDiagonalMatrixpath here:mathnet-numerics/src/Numerics/LinearAlgebra/Double/DenseMatrix.cs
Line 611 in 4b13ff4
This PR does that in the second commit. The first commit tries to increase test coverage to make sure this path is tested. I wanted to add some tests to check arithmetic results are equal across different storage types (a generalisation of #912), I didn't achieve that but I think the changes add a little extra robustness in that area.
I ran a rudimentary benchmark in a distinctly unscientific environment (third commit, you may not want to take that) just for fun: