Kratos: Matrix Comparison

[roigcarlo]

Dear all,

I've recently noticed that the number of test that need to compare matrices has been steadily increasing over the last few months.

We don't have any natural way to perform this operation, in fact I've run some performance tools and I suspect that the way we currently do it, cell per cell, has a negative impact in the compile times.
With the inclusion of the AMatrix library I think it is a good moment to think about implementing one check specific for matrices and vectors.

I am not a heavy user of matrices in Kratos so I would like some feedback on how it would be the most natural way for you to compare this and how you would expect it to behave. I pretend to start with non sparse matrices for the moment and extend it in case is useful, so only dynamic and static for now.
I mainly have three concerns.

• How the reference matrix would be created?
• Which is the interface you would expect such a comparison have?
• How do you want the error message to be, should we print only the first difference, should we print them all? is the value important, only the position, both?

In order to start to collect ideas consider this:

// The calculated matrix
Matrix<3,3> M = CalculateSomething();

// Comparing to an std::vector
std::vector<double> Rv = {1, 2, 3, 4, 5, 6, 7 ,8, 9};

// Comparing to a Matrix with an list initializer 
Matrix<3,3> Rm = {1, 2, 3, 4, 5, 6, 7 ,8, 9};                // Is this possible to do?

// Comparing to a Matrix without initializer 
Matrix<3,3> Re;
Re(0,0) = 1;  Re(0,1) = 2;  Re(0,2) = 3;
Re(1,0) = 4;  Re(1,1) = 5;  Re(1,2) = 6;
Re(2,0) = 7;  Re(2,1) = 8;  Re(2,2) = 9;

KRATOS_CHECK_MATRIX_EQUAL(M, R);
KRATOS_CHECK_MATRIX_EQUAL(M, R, 3, 3);
KRATOS_CHECK_MATRIX_EQUAL(M, R, 2, 2);

• What would be the most natural way for you to generate the reference matrix?
• Would you expect to be able to check std::vector<std::vector<double>> against Matrix?
• Would you expect to be able to check std::vector<double> against Matrix?
• Would you expect the check to validate the size of both matrices, and in this case, would it be useful/impactful to perform as much of the check as possible (min of the sizes)?
• Would it be interesting to specify the size of the check, ex [2,2] in a [3,3] matrix?
• If so, would it be interesting to be able to compare sub-matrices?

Please @KratosMultiphysics/team-maintainers ping anyone interested in the discussion.

@pooyan-dadvand since you are implementing AMatrix.
@loumalouomega , @rubenzorrilla I know you are intensive users of this, so you may be interested.

[bodhinandach]

This is a good idea @roigcarlo. I just created the tests in ParticleMechanicsApplication and think that will be good to have this feature to check matrices and vectors. I have a few opinions:

• I would think that KRATOS_CHECK_MATRIX_NEAR is also needed as sometime you have different results due to precision. Therefore, KRATOS_CHECK_MATRIX_NEAR(matrixA, matrixB, tolerance) will also be usefull.
• If this is just for testings, I think type should be the same, throw an error otherwise.
• Also if the size of both matrices are not the same, an error is expected.
• For the error difference, I think having all of the differences printed out is better, so that we can clearly see the overall problems. Otherwise you need to compile and recompile to check it one by one and therefore time consuming.

[loumalouomega]

In #3236 I compare the nonzero values of a sparse matrix, which to compare is not quite straightforward. We can do the difference and compute the norm, but that is a costly operation.

For the dense small matrices I am OK with computing the norm in a macro, as it is not so expensive

[loumalouomega]

I agree with the proposals of @bodhinandach

[roigcarlo]

In #3236 I compare the nonzero values of a sparse matrix

That's why I explicitly stated that I do not want this operation to be valid for sparse matrix for now :smile:

[pooyan-dadvand]

The exact equality should work with KRATOS_CHECK_EQUAL as both ublas and AMatrix are defining the == operator.
Nevertheless the KRATOS_CHECK_MATRIX_NEAR is not defined as would be a very useful macro.
About the generations AMatrix allows the initialization with array of double.

[GuillermoCasas]

About matrix proximity, there is not a single way to do it. However, I think that comparing the norm of the difference to the tolerance, as @loumalouomega suggested, is one consistent way to do it that probably covers much of what researchers want to do.

To start off with, I would include the -norm family https://en.wikipedia.org/wiki/Matrix_norm#L2,1_and_Lp,q_norms, perhaps only the , and (the latter is the max norm).

The call proposed by @bodhinandach could be generalized by adding an optional last parameter (the norm tag) lie this: KRATOS_CHECK_MATRIX_NEAR(matrixA, matrixB, tolerance, '2') . If this is complicated, one could just define several macros: KRATOS_CHECK_MATRIX_NEAR_L1 and so on.

Note that one can implement the check in a loop, so that as soon as the accumulated norm is bigger than the tolerance it stops calculating.

I agree with @bodhinandach that it is ok to require same sized, equal-typed matrices. I would even go further, not printing anything. Why should we overload so much this macro? I would just make it do what it says it does: assert if the matrices are equal.

I agree that many of the capabilities mentioned above would be nice to have, but perhaps they can be implemented elsewhere. For instance, about comparing sub-matrices, why not just do
KRATOS_CHECK_MATRIX_EQUAL(M(0,0,3,3), R(0,0,3,3))? (sorry, have to check the syntax of Amatrix)

[GuillermoCasas]

I am having second-thoughts about printing or not. Now I think maybe it is best to print the norm of the error, along with the difference matrix in case of error (only compute these things in case of error).

[pablobecker]

but two matrices might have the same norm while being completely different..
for matrices to be defined as almost equal, IMO the check should be element wise.

[GuillermoCasas]

The following points summarize the current consensus of the @KratosMultiphysics/implementation-committee about this issue:

1. Support for submatrices should be provided from the outside.
2. The macro should check matrix size and print the corresponding specific error for that case
3. Any type of matrix should allowed (no need for equal types)
4. There should be a macro for exact equality that applies the existing macro for numbers entry by entry.
5. There should be another macro (or set of macros) that implement the tolerance-based approximate equality in a specified norm (the norm of the difference matrix, that is). It is suggested that support is given for the norms mentioned in the comment by @GuillermoCasas above.

Passing the question to the @KratosMultiphysics/technical-committee for the final approval.

[jcotela]

@philbucher made me notice that PR #5430 partially implements this issue. I will try to summarize what is in that PR in relation to the discussion here:

• Macros to check for equality (and "closeness" based on a tolerance) are added for dense matrices. They check first that both matrices have the same dimensions and then check that the difference term-by-term is not larger than the provided tolerance.

• In case of failure, the error message specifies the cause of the error, prints the entire matrices and indicates the first term were the tolerance was not respected.

• The matrix/vector types are not relevant for this check, they don't even need to be the same for the two arguments. The only requirements (which both amatrix and ublas satisfy) is that matrices have size1(), size2() methods and operator() based indexed access.

• The matrix equality is for now an overload of the matrix near macro, instead of being based on operator==. This may not be as efficient, but I believe it is more informative (operator== should not be expected to give a detailed report on the differences between its arguments).

• There is no support for sparse matrices. As @roigcarlo mentioned, that would be a completely different beast.

• There is no support for checking only submatrices (but I agree with @GuillermoCasas and the @KratosMultiphysics/implementation-committee that this should come from the outside, not be part of the macro).

• There is no support for norm-based checks.

[roigcarlo]

@KratosMultiphysics/technical-committee closes this one as the basic macro is implemented and other tolerance based ones can be don in a separate one when AMatrix is merged.