Re: comments on Matrix

On Thu, Mar 21, 2013 at 7:38 AM, Benoit Jacob <jacob.benoit.1@gmail.com>wrote:

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

I don't question your knowledge of matrices.
However, this does not seem to be a problem in the field.


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