Re: comments on Matrix

2013/3/21 Benoit Jacob <jacob.benoit.1@gmail.com>

>
>
> 2013/3/21 Rik Cabanier <cabanier@gmail.com>
>
>>
>>
>> On Wed, Mar 20, 2013 at 8:38 PM, Benoit Jacob <jacob.benoit.1@gmail.com>wrote:
>>
>>>
>>>
>>> 2013/3/20 Rik Cabanier <cabanier@gmail.com>
>>>
>>>>
>>>>
>>>> On Wed, Mar 20, 2013 at 6:29 PM, Benoit Jacob <jacob.benoit.1@gmail.com
>>>> > wrote:
>>>>
>>>>>
>>>>>
>>>>> 2013/3/20 Rik Cabanier <cabanier@gmail.com>
>>>>>
>>>>>>
>>>>>>
>>>>>> On Tue, Mar 19, 2013 at 6:38 PM, Benoit Jacob <
>>>>>> jacob.benoit.1@gmail.com> wrote:
>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> Seeing that a matrix API was being discussed (
>>>>>>> https://dvcs.w3.org/hg/FXTF/raw-file/default/matrix/index.html), I
>>>>>>> thought I'd take a look and chime in.
>>>>>>>
>>>>>>> Here are the first few things that scare me, glancing at the current
>>>>>>> draft:
>>>>>>>
>>>>>>> 1. Some functions like inverse() throw on singular matrices. The
>>>>>>> problem is it's impossible to define very firmly what singular means, in
>>>>>>> floating-point arithmetic --- except perhaps by prescribing a mandatory
>>>>>>> order in which the arithmetic operations inside of inverse() are performed,
>>>>>>> which would be insane --- so saying that inverse() throws on singular
>>>>>>> matrices means in effect that there are realistic matrices on which it is
>>>>>>> implementation-defined whether inverse() throws or not. The consensus in
>>>>>>> all the serious matrix libraries that I've seen, for closed-form matrix
>>>>>>> inversion, is to just blindly perform the computation --- in the worst case
>>>>>>> you'll get Inf or NaN values in the result, which you will have anyway on
>>>>>>> some input matrices unless you mandate unduly expensive checks. More
>>>>>>> generally, I would never throw on singular-ness in any function, and in the
>>>>>>> case of 4x4 matrices and closed-form computations I wouldn't attempt to
>>>>>>> report on singular-ness otherwise than by inf/nan values.
>>>>>>>
>>>>>>
>>>>>> Are you suggesting that the user should check all the values to make
>>>>>> sure that the matrix inversion didn't fail? That seems very expensive.
>>>>>> Can you point us to an algorithm for inversion that you think the
>>>>>> spec should include?
>>>>>>
>>>>>
>>>>> No, I'm saying that we should _not_ check anything, or if we do (say
>>>>> in a separate "checked inverse" method), it should be in a graceful way
>>>>> e.g. an output boolean parameter, not throwing an exception which by
>>>>> default halts the program. That's what I meant by "The consensus[...] is to
>>>>> just blindly perform the computation --- in the worst case you'll get Inf
>>>>> or NaN values in the result" above.
>>>>>
>>>>
>>>> That's certainly tempting! I can see how that's easier (and less
>>>> disrupting).
>>>>
>>>> I checked 2 of our internal libraries, as well as Cairo, Skia and Flash
>>>> and they all refuse to calculate the matrix if it's not invertible. So, not
>>>> all matrix libraries do as you suggest.
>>>> SVGMatrix is also documented to throw.
>>>> These are all matrix classes designed for graphics so I *think* it
>>>> makes sense to follow what they did.
>>>>
>>>
>>> I had dedicated matrix libraries in mind, rather than the casual
>>> incidental matrix computation in something like Cairo.
>>>
>>
>> So, are you saying the matrix class should not throw?
>>
>
> Yes. Matrices are just a generalization of numbers (which are 1x1
> matrices). Numbers don't throw --- not even 1/0 or 0/0. Matrices should
> not, either.
>
>
>>
>>
>>>
>>>
>>>>
>>>>
>>>>>
>>>>>
>>>>> The algorithm for inversion is going to be plain closed-form matrix
>>>>> inversion (i.e. by computing the cofactors) as in
>>>>> http://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3.C3.973_matricesbut instead of mandating implementers to use particular formulas, a good
>>>>> matrix API would rather carefully avoid being too sensitive to the exact
>>>>> floating-point values (e.g. by not throwing on singular matrices) because
>>>>> even mandating formulas won't ensure that the results are the same
>>>>> everywhere up to the last bit of precision.
>>>>>
>>>>
>>>> I don't know how important this is as it seems like this would only
>>>> apply to edge cases.
>>>> Are there any w3c specs that define how precise a calculation should
>>>> happen? I can see how this would be very important for mathematical
>>>> programs but less so for graphics.
>>>>
>>>
>>> That's my point: the last bits of mantissa are irrelevant to graphics,
>>> so we should just make sure that our API doesn't accidentally make them
>>> unduly important (e.g. by throwing on singular matrix or trying to avoid
>>> methods like is2D() that can return true or false depending on them). Once
>>> we've ensured that, we really don't need no worry about different browsers
>>> giving slightly different values.
>>>
>>
>> I'm a bit confused.
>> It seems that if the last bits are unimportant and the determinant
>> returns something that is in the range of those small bits, 'invert' should
>> fail and throw.
>> If not, you might generate an inverted matrix with very large values and
>> things will go very wrong...
>>
>
> If the determinant is only a few times machine epsilon then the result is
> indeed unreliable, i.e. it will sometimes be very imprecise.
>

I need to qualify this a bit:

If the determinant is only a few times machine epsilon TIMES a quantity
that depends on the order of magnitude of the matrix coefficients then the
result is unreliable.

Let me explain this:

Let's use floats. Machine epsilon is roughly 1e-7 and the range of absolute
values is roughly from 1e-30 to 1e+30. The following matrix:

  1e-10 0
  0     1e-10

is perfectly fine (not even close to approaching the boundaries of the
range of float values) and has determinant 1e-20 which is far below machine
epsilon. Yet it is perfectly invertible, and its inverse is

  1e10 0
  0    1e10

So just comparing the determinant to machine epsilon doesn't work in full
generality. It might work if the matrix coefficient were on the order of
magnitude of 1 but that's not even the case of graphics matrices, were
coefficients can typically be of the order of magnitude of the size of the
page in CSS pixels, I suppose, which would be about 1e3.

Another issue is that it should be clear that any determinant check is at
most a necessary condition, not a sufficient condition for the inverse to
be reliable. For example, take this matrix:

1 0     1e4
1 1     1e4
1 1e4 1

If we compute its determinant by developing it along the first column, as
it is most usually computed, we get

  1 - 1e8 + 1e8 + 1e4

which is mathematically equal to 1e4 + 1, but a with floats, 1 - 1e8 == -
1e8 so we get 1e4 instead.

That's a relative imprecision of 1e-4 only, but that's already over 1e3
time machine epsilon, and that was only a very naive example, and the
determinant of 1e4 would have passed any check guarding against very small
(or very large) determinants.

Just to make it clear that simple checks like checking against small
determinants only buy you so much.

Benoit


> But if on that basis one decides to do something like
>
> function inverse() {
>   if (determinant() < THRESHOLD) {
>     // return an error
>   }
>   // compute and return the inverse
> }
>
> Then one creates a huge discontinuity at the threshold point, which will
> disrupt a majority of application much more than this risk of imprecision.
>
> It is legitimate for an application to say that they really want this
> check to be done. In this case, if you want (to maximize efficiency) to do
> this without computing the determinant twice, you can add a
> checkedInverse() method returning not only the inverse but also a boolean
> indicating whether it is to be trusted.
>
> But:
>  - please don't throw because that defaults to killing the application,
> which is not something that should happen as a result of floating point
> imprecision.
>  - be very wary of hardcoding an arbitrary fixed THRESHOLD. There's no
> one-size-fits-all there. That's what we learnt the hard way in Eigen and we
> eventually decided to let the only inverse-with-check function we have,
> take a 'threshold' argument (though it has a vaguely reasonable default
> value).
> http://eigen.tuxfamily.org/dox/classEigen_1_1MatrixBase.html#abcd1eb30af2979b147e1f65477caa209
>
> Benoit
>
>
>
>>
>>
>>> That's a good thing, because it would be completely unrealistic to
>>> expect exact values. Keep in mind that in floating point, a*(b*c) and
>>> (a*b)*c give different results i.e. associativity doesn't even hold. Same
>>> for distributivity: a*b+a*c != a*(b+c). So for example just a single
>>> browser having a SSE and a non-SSE codepath would almost certainly get
>>> different results from them.
>>>
>>> Benoit
>>>
>>
>>
>