Re: comments on Matrix

On Thu, Mar 21, 2013 at 11:21 AM, Rik Cabanier <cabanier@gmail.com> wrote:

>
>
> On Thu, Mar 21, 2013 at 6:01 AM, Benoit Jacob <jacob.benoit.1@gmail.com>wrote:
>
>>
>>
>> 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.
>>
>
> I think that's beyond the grasp of most authors. Let's not create an
> interface that requires a degree in mathematics.
>
> Direct2D also refuses to do the inversion:
> http://msdn.microsoft.com/en-us/library/windows/desktop/dd372277(v=vs.85).aspx
> Given that all the libraries refuse to invert and all of Adobe's graphics
> applications (such as Photoshop and Illustrator) do the same under the hood
> with no ill effects, let's not do this by default.
>

Direct2D doesn't kill your app when the inversion fails. It just returns
failure which you are free to ignore.


>
> If you feel very strongly, we can provide another call that never throws.
> We would have to change the IDL of the matrix interface though since it
> could now contain unrestricted doubles.
>
>
>>
>>
>>>
>>>
>>>>
>>>>
>>>>>
>>>>>
>>>>>>
>>>>>>
>>>>>> 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.
>>
>
> Why would a throw kill an application? Just catch it and abort the
> (unstable) calculation.
>

Because most authors will not realize that everywhere they call inverse
they need a catch and when the user does something (like set fov to 0) the
app will crash and the user will lose work.  If the function does not throw
then the app just keeps running, the user, seeing a mess on the screen from
setting the fov to 0 sets it back to a reasonable value and continues
working.




>
>
>>  - 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
>>
>
> This would be up to the implementation to decide. The spec can be silent
> about this.
>
>
>>
>>
>>>
>>>
>>>> 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
>>>>
>>>
>>>
>>
>