Re: comments on Matrix

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.

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
>>
>
>