Mathjs: QR decomposition -- zero assert is misleading / useless?

[c-dante]

The assertion here that R is lower triangular and that it "has to be close to zero" doens't really add anything.
https://github.com/josdejong/mathjs/blob/master/lib/function/algebra/decomposition/qr.js#L229

In V8 and chrome, the idea of dividing by a large number to coerce floating point values to be zero is not a great check -- and the if/throw is just useless. Is there a good reason to keep this in? Why feel uneasy about coercing it to zero? Floating point is inherently going to be imprecise at times.

If you REALLY want a "close enough" check, doing an epsilon ball is good enough (especially if it can be configured).

[harrysarson]

I wrote this function :)

I am worried that a large number could appear here by some error or edge case which is why I want to check. I also thought that a self respecting qr decomposition would need to actually be upper triangular.

[c-dante]

I can get that -- what I'm offering is that floating point numbers are almost always going to lead to cases where sufficiently small numbers appear in the lower part of the R matrix. It's incredibly easy to get offset with floating point math (looking at you 0.2 + 0.1), and the "close to zero" check via division isn't a great way to do that kind of assertion.

In general, you're better off with an epsilon ball -- this is because, as an example:

1e-10 / 1e5 // 1e-15
1e-10 / 1e5 === 0 // false

I would absolutely argue 1e-10 is sufficiently close to zero for me to be happy with the results of the QR decomposition (within floating point math error).

The fact that this "close to zero check" throws an exception means using this QR decomposition will just _fail_ for a bunch of matrices. It's also tough to just assert what you think is "close enough" for other people's applications/tolerances. Perhaps something like:

math.qr(matrix [, epsilon = 1e-5]) {

// ...

    for (i = 0; i < rows; ++i) {
      for (j = 0; j < i && j < cols; ++j) {
        if (Math.abs(Rdata[i][j]) > epsilon) {
          throw new Error('math.qr(): Value exceeded allowed epsilon - ' + 
           'R is not lower triangular (element (' + 
            i + ', ' + j + ')  = ' + Rdata[i][j] + ')'
          );
        }
        Rdata[i][j] = multiplyScalar(Rdata[i][j], 0);
      }
    }

// ...

[harrysarson]

The equal function I used uses an epsilon value.

[c-dante]

I'm suspect of the logic behind nearlyEqual -- that is, it may use an epsilon, but it fails hard for a lot of cases that I'd expect it to pass. That is, it's not a "within an epsilon" check -- it's more to do with addition of floating points, so it fails on something like 0.1 and 0.11 (when you get VERY small):

const nearlyEqual = (x, y) => Math.abs(x - y) <= Math.max(Math.abs(x), Math.abs(y)) * Number.EPSILON
nearlyEqual(0, 0) // true
nearlyEqual(0, 1e-100) // false, but I'd expect that to be true (?)
nearlyEqual(1e-200, 0) // false, but I'd expect that to be true (?)
relError(1e-200, 2e-200) // false, but can expect that -- 1 !== 2
relError(1, 1 + Number.EPSILON) // true -- this is the exact expected case
relError(1e-50, 11e-51) // false -- this is strange because it's effectively a sig fig problem

[harrysarson]

I am out of my depth here I think. I might refer up @josdejong

[josdejong]

@c-dante thanks for all your input!

The _nearly equal_ implementation uses in math.js calculates the relative difference, described in more detail in the docs of mathjs here: http://mathjs.org/docs/datatypes/numbers.html#equality. As far as I know this is quite a common way to compare values, but I'm always open to reconsider earlier choices. I don't remember exactly what the pros and cons of various comparison mechanisms where again, but we could list them again.

As for the check against numerical instability in the qr function: If we're certain that a solution totally makes no sense it's indeed a good idea to protect the user against this. But it's really tricky to determine that and there is a large gray field of "precise enough" results (all calculations suffer from some loss of precision). Also, when using the qr function with BigNumbers having a high precision, this check may be too conservative. @HarrySarson do you think this warning is really critical? I think we could remove it if needlessly restricts the function too much - it's hard for me to estimate whether this check does more harm than good, but I agree with @c-dante that it's also partly up to the user to deal/reckon with instabilities when working with matrix operations in general.

[harrysarson]

So a possible solution would be to separate the qr function from a new function that also zeros out elements.

Then we can test the qr function without worrying about missing hidden errors which was my greatest concern.

Missing errors in the tests was my great concern so I would be happy with this change

[josdejong]

👍 would be great if you could give this a try and rework the qr function Harry.

[harrysarson]

@josdejong could you assign this to me? Otherwise I might forget about it again.

[josdejong]

@harrysarson Sure! If you are willing to become a collaborator on the project :smiley: I've send you an invitation. Can you share me your email too privately (you find mine on my profile page)?

[harrysarson]

Yep have done :)

[josdejong]

This issue is fixed now in v6.2.4. Thanks again @harrysarson!