#1147 Norm of NaN-matrix

Here's a patch that addresses the aforementioned problem under the assumption that dim > 0 (is this legitimate?).

Attachments

0001-infinity_norm-Make-fold1-into-a-proper-fold.patch

Now with fewer pointless comparisons.

Attachments

0001-infinity_norm-Make-fold1-into-a-proper-fold.patch

And now a patch that is actually different from the first one.

Attachments

0001-infinity_norm-Make-fold1-into-a-proper-fold.patch

Thanks. I wouldn't want to rule out the case that somebody calls this for a 0x0 matrix. In that case the old code returned 0, and your new code will crash. Can you adapt your patch such that it returns 0 on a 0x0 matrix?

Sure, here's a new version.

Attachments

0001-infinity_norm-Make-fold1-into-a-proper-fold.patch

Or rather, here.

Attachments

0001-infinity_norm-Make-fold1-into-a-proper-fold.patch

Patch has been applied in revision 6819. Many thanks! Not surprisingly, the same bug exists in densevector. Could you be so kind and provide a patch for that too?

Sure.

Attachments

0001-Fix-FieldVector-infinity_norm-in-analogy-to-r6819.patch

Here's a test for the changes to fvector and fmatrix.

Attachments

0001-Add-a-test-for-FS-1147.patch

Forgot an include (for the test).

Attachments

0001-Add-a-test-for-FS-1147.patch

First things first: I introduced a subtle bug in FieldMatrix::infinity_norm_real() that's already been pushed. Sorry about that. Also, the list of patches above is getting a bit confusing (darn you, flyspray!).

So I'll attach three patches to this comment:

0001-Fix-FieldVector-infinity_norm-in-analogy-to-r6819.patch: Addresses the Nan-problem for infinity_norm() and infinity_norm_real() of the FieldVector class

0002-Add-a-test-for-FS-1147.patch: Adds tests for the NaN-handling of FieldVector and FieldMatrix

0003-Fix-up-r6819-Add-a-test-case.patch: Fixes the bug I mentioned in FieldMatrix::infinity_norm_real(). Adds corresponding tests, too.

Attachments

To be clear: The three patches in my previous comment obsolete all the previous patches that have not been pushed.

infinity_norm() and infinity_norm_real() share quite a bit of code. I've factored that out into a private infinity_norm_generic.

So here's a fourth patch.

Attachments

0001-Make-infinity_norm-generic.patch

Same patch, without a typo:

Attachments

0001-Make-infinity_norm-generic.patch

I'm not 100% happy with the current patch (the one Oli applied...).

Comparing the iterators will be rather difficult to statically evaluate the "if". When using size(), it should me easier and not less readable.

@christi: Did that and rebased the other patches on top of that. That leaves five patches (please ignore all the others):

0001-Do-not-compare-iterators-Let-s-hope-the-optimiser-ca.patch: Amendment to the patch that's already been pushed in response to your (Christian's) critique 0002-Fix-FieldVector-infinity_norm-in-analogy-to-r6819.patch: Addresses the NaN-Problem for fieldvectors 0003-Add-a-test-for-FS-1147.patch: Adds a test for both vectors and matrices. Catches the NaN problem and the following bug. 0004-Fix-up-r6819-Add-a-test-case.patch: Fixes a bug I introduced in FMatrix::infinity_norm_real() through the NaN-Fix that's already been pushed 0005-Make-infinity_norm-generic.patch: An attempt to factor out code and keep bugs like the one mentioned above from happening. Not exactly pretty, though.

Attachments

gcc 4.7 doesn't seem to like patch 0005.

Is anything holding back patches 1-4?

I just committed the first two, plus the fix-up part of 0004. Thanks a lot! I would push 0005, but indeed it doesn't compile with gcc-4.7. Any ideas?

It is nice to have the tests from 0003-0004. However, traditionally we have grouped all tests for a single class into a single file. Hence can I ask you to merge you new tests into the files that test FieldMatrix and FieldVector.

Patches attached.

Attachments

I don't understand why gcc 4.7 doesn't like patch 0005 either btw.

Just goes to show that patch is evil, I guess -- hardly any gain and it makes error diagnosis a lot harder.

So let's forget about patch 0005. The rest is committed now. Thanks!

LAPACK returns NaN for the matrix norm if any entry is NaN(*). I think that's sane. Shouldn't we do the same?

(*) http://www.netlib.org/lapack/lapack-3.4.2.html under "More details"

I think returning NaN is a good idea.

@christi: Just to be clear: You agree that a matrix/vector, for which (at least) a single entry is NaN, should have norm NaN?

It seems this situation is actually quite a bit more complex than I thought it was.

Leaving matrices aside for now, what this bug was initially about is essentially this scenario

assert(std::numeric_limits<double>::has_quiet_NaN);
double mynan = std::numeric_limits<double>::quiet_NaN();
using V = Dune::FieldVector<double, 3>;
V v = { mynan, mynan, mynan };
v.infinity_norm();

where the last line would previously yield 0 and now NaN (which I think is an improvement).

The way I came across this issue was that I computed an error, my code told me the error had infinity norm 0 yet I could clearly see that my code was misbehaving. It turned out that the error computation had failed completely and the error vector consisted of NaNs only. From my point of view, a norm of 0 did not make any sense. For someone who uses NaN to signify incomplete (or in this case entirely missing) data, I take it the 0 norm actually would.

Now here's the issue: While the infinity norm computation for vectors with no or exclusively NaN entries is now in some sense sane, anything inbetween is highly surprising to me.

My expectation was that NaN entries would simply be ignored (again, undesirable given the current all-NaN-behaviour), since I assumed the following to hold true: std::max(x,NaN) and std::max(NaN, x) both return x. Yet somehow e.g. std::max<double> behaves quite differently from std::fmax, creating the following bizarre situation:

V v = { mynan, mynan, mynan }; v.infinity_norm(); // nan
V v = { 1, mynan, mynan }; v.infinity_norm(); // 1(!)
V v = { mynan, 1, mynan }; v.infinity_norm(); // nan
V v = { mynan, mynan, 1 }; v.infinity_norm(); // nan

I've reduced this a bit to the following behaviour:

    assert(std::numeric_limits<double>::has_quiet_NaN);
    double mynan = std::numeric_limits<double>::quiet_NaN();

    std::max<double>(mynan, 2.0); // nan(!)
    std::max<double>(2.0, mynan); // 2

    std::fmax(mynan, 2.0); // 2
    std::fmax(2.0, mynan); // 2

Since std::max_element uses std::max internally, it, too, returns different values for the std::vectors { 2, mynan } and { mynan, 2 } by the way.

I'm not sure if there's a sane way to handle all this without explicit calls to isnan for every entry, which would be far too costly... Maybe we should actually leave NaN handling entirely to the user (again)?

The issue of max/min and nan is nicely explained here: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=47706

(Deleted my previous comment because it was entirely wrong)

So the upshot is then: std::max makes assumptions of comparability that special values of double violate, hence std::max can only be used with finite double values.

If the user knows her field type to be double, then using std::fmax(x,y) instead of std::max(x,y) has the advantage of being invariant under x-y-swapping. In contrast to what I believe we want here, std::fmax treats NaN values as missing data, so that it, too, cannot be used here, though...

So maybe the only workable solution is to assume the user will sanitise her data, so that vectors will never contain NaN values and continue to use std::max (plus revert my changes, which would make the code a lot nicer again).

This bug has gotten rather messy over time and I feel that I should state more clearly what it's really about and why you might want to care.

For a vector, there are currently two types of norms: norms that sum entries somehow, like the one_norm() and the two_norm(), and a norm that doesn't, namely infinity_norm. Because any sum that contains NaN somewhere will be NaN, the summing norms will always return NaN as soon as a single vector entry is NaN. This is not so for the infinity_norm(), which I find inconsistent.

I made the situation infinitesimally better by supplying a fix to infinity_norm() for the case where the vector consists solely of NaN entries (thinking back then that this was the only situation that caused problems). At the same time, I made things a bit worse because now the order of the vector entries actually matters. The vector { NaN, 2.0 } will have norm NaN whereas { 2.0, NaN } will have norm 2.0. This is because instead of std::fmax, the current implementation uses std::max. Now one might suggest a switch to std::fmax for floating point types, which would remove the order dependence, but it would still be so that any vector that contains even a single proper number would not have infinity_norm() return NaN.

One could of course role own's own mymax() that applies isnan() to both of its arguments. I don't know how costly such calls would be, but they would actually be a bit wasteful here because twice as many NaN checks as necessary would be carried out in the infinity_norm() loop this way. One could instead check for NaN in the loop itself but would have to specialise this function for floating point types then, and finally one would have to answer if all of this is really worth it or if the user should really make sure her vectors are free of NaNs (this kind of added security check could also be enabled through a define).

I did some experiments and measured some timings. If you do it correctly, you will have no performance penalty compared to the "standard" implementation (below named version 0).

The code can be found at https://gitlab.dune-project.org/snippets/11

After compiling with g++-5.2 -O3 I get the following timings, these do not change when using sse4.2. They might differ somewhat for other architectures:

> ./test_inf_norm
Verify Implementations 
  Version	0 failed
  Version	1 works
  Version	2 works
  Version	3 works
  Version	4 works
  Version	5 works
  Version	6 works
----
v = double[1000000000]
v[999999999] = NaN
  Version 0	1.65497 ms
  Version 1	3.54432 ms
  Version 2	2.23769 ms
  Version 3	2.30446 ms
  Version 4	0.965323 ms
  Version 5	1.28024 ms
  Version 6	1.27563 ms
----

The funny thing is that with clang-3.7 the timings are significantly different. Apperently gcc is much better at optimizing the code:

v = double[1000000000]
v[999999999] = NaN
  Version 0	5.18355 ms
  Version 1	3.81165 ms
  Version 2	4.44596 ms
  Version 3	4.46749 ms
  Version 4	4.52507 ms
  Version 5	4.76469 ms
  Version 6	4.75915 ms

With g++-4.9 I get roughly the same timings as with g++-5.2. Looking at the numbers, version 1 is never really fast. On the other hand, g++ manages to get improved performance for a couple of implementations. Version 4 is the fastest for g++, while version 5 and 6 will have a better performance, if the NAN is not the last value.

Oh! And the numbers do not change if I have no NAN entry. I even tested a version of the code, where every value in v is initialized with a random value, but again this gave to the same timings, just that 90% of the time was consumed by the std::rand.

In summary I suggest to implement version 4. I tracks the NAN by computing the sum of all values and at the end corrects the max with the NAN.

@christi: Is version 4 correct? On first glance, the norm of the null vector will be nan (because 0 / 0 should be nan).

@martin.nolte Oh you are right, this looks like a corner case. I'll update the test accordingly.

I updated the test and added a 7th version, which is basically the the same version as 4, but instead of computing a normalization, it explicitly tests for NAN. I also updated version 4 to start the sum with 1 ans sums abs(v[i]), so that it doesn't matter if the whole vector is 0.

I further added a second test, which splits the vector into tiny chunks of size 10.

In the first timing test version 7 and version 4 perform comparably, but in the second test (which is much harder) version 4 takes ~1.5 times longer than version 0 and version 7 take twice as long. All other versions are considerably more expensive.

The funny thing is that with clang-3.7 the timings are significantly different. Apperently gcc is much better at optimizing the code:

This might be related to https://llvm.org/bugs/show_bug.cgi?id=25566 which suggests that one should benchmark with clang 3.5 rather than with 3.7 (I cannot comment on 3.6).

@pipping Interesting ... and yes, with clang-3.6 everything has roughly the same timings as with g++. What still puzzles me is that in the 10^9 entries test run, I get better timing for version 4 than for the (wrong) standard implementation.

The other funny thing is that for clang version 7 is slightly faster than version 4, while for g++ version 4 was always faster, sometime even a factor 1.5.

I think that's related to caching. If I reverse the order, in timing_tests and use

for (int version = 7; version >= 0; version--) {

with with g++ 4.9, I get something like

  Version 7	nan 	1.47424 ms
  Version 6	nan 	1.08488 ms
  Version 5	nan 	1.08359 ms
  Version 4	nan 	0.813823 ms
  Version 3	nan 	1.00453 ms
  Version 2	nan 	1.00366 ms
  Version 1	nan 	2.98307 ms
  Version 0	nan 	0.815933 ms

so that version 4 is longer faster than version 0.

I've changed the snippet such that it loops through version 0 to 7 and then 7 to 0 and uploaded it as https://gitlab.dune-project.org/snippets/12

The second iteration appears to give more reasonable timings.

Finally, it might be worth mentioning that the differences between version 4 and version 7 essentially disappear for me (both with gcc and clang) if I compile with -O2 -march=native.

@pipping

you are absolutely right, the caching is an issue here, I'd assume running version 0 once to avoid cold caches and then doing the benchmark should be sufficient.
Today, when explaining the issue to Jö, I had an other idea. I'll update the test try it out.

I just performed a test for boost::rational, which does not have a NaN value. Version 4 obviously goes through, but the performance loss is really big. Perhaps we will have to specialize for field types with NaN and for types without. If we do so, we could instead of version 4 also use version 7.

Hm. I'm a tad surprised now because the latest Boost rational documentation seems to state

This library does not support representations for infinity or NaN

under Exceptions.

Update: Never mind, I misread entirely what you wrote. What should we specialise to then? Certainly anything for which std::is_floating_point<T>::value is true. But also mpfr_float from boost/multiprecision/mpfr.hpp?

The problem with version 7 is that there is no isnan(std::complex<double>) (which we don't need for version 4).

I've now removed version 7 for said reason and extended the snippet to cover std::complex<double> (Please see https://gitlab.dune-project.org/snippets/13). Here, it really doesn't seem to matter all that much which version one goes with.

Combined with the double performance, this suggests that we might want to specialise infinity_norm() and use version 4 for floating point types.

mentioned in merge request core/dune-common!15 (merged)

mentioned in commit core/dune-common@1f6422da

This issue is solved with the latest patches by Elias.

Status changed to closed

mentioned in merge request core/dune-istl!10 (merged)

As suggested in core/dune-common!15 (merged), continuing here with an extension to dune-istl classes...

Oops... and I've broken the compilation of dune-common for everyone. Please kindly push the following fix, too

0001-Add-missing-include.patch

mentioned in commit core/dune-istl@c8ad7196

mentioned in commit pipping/dune-istl@c8ad7196

mentioned in commit pipping/dune-common@1f6422da

mentioned in merge request core/dune-istl!15 (merged)

mentioned in commit core/dune-istl@c9ee3bdf

mentioned in commit felix.gruber/dune-common@1f6422da

Mentioned in commit core/dune-common@1f6422da

mentioned in issue core/dune-common#72

Property	Value
Reported by	Elias Pipping (elias.pipping@fu-berlin.de)
Reported at	Jul 6, 2012 09:42
Type	Bug Report
Version	Git (pre2.4) [autotools]
Operating System	Unspecified / All
Last edited by	Elias Pipping (elias.pipping@fu-berlin.de)
Last edited at	Oct 20, 2012 12:05

#1147 Norm of NaN-matrix

Metadata

Description

Designs

Child items ...

Activity

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments

Attachments