[scilab-Users] Problem with Floating point and determinant of a matrix

[michael.baudin at scilab.org]

 Hi Reinaldo,

 As Gilbert Strang writes, "The way to understand this subject is by 
 example". Here is an example where the matrix A is numerically well 
 conditioned, but has a zero determinant :

 -->n=200;
 -->a=1.e-2;
 -->A=a*eye(n,n);
 -->det(A)
  ans  =
     0.
 -->rcond(A)
  ans  =
     1.

 Hence, we have just found a 200-by-200 matrix A for which the 
 determinant is zero, even if the matrix A is perfectly conditioned. In 
 practice, this particular matrix A behaves much more as a diagonal 
 matrix than as a singular matrix. You can find a similar example in 
 "Computer methods for mathematical computations", by Forsythe, Malcolm 
 and Moler, see chapter 3 "Linear systems of equations", p. 44.

 The kind of conditioning which is computed by rcond is the conditioning 
 with respect with the linear equation A*x=b problem.

 The determinant is a bad measure of the numerical singularity of the 
 matrix A. In some sense, it is "absolute" and does not take into account 
 for scale factors.

 The condition number defined by K(A)=||A|| ||A^-1||, for a given norm, 
 is a measure of the closeness to the set of the singular matrices. For 
 the A*x=b problem, it is a relative error magnification factor.

 All in all, we may write our algorithms as :

 if ( rcond(A) > threshold ) then
   // A is well conditioned
 else
   // A is ill-conditioned
 end

 A reasonable threshold may be threshold = %eps ~ 10^-16. But this 
 condition is only valid for the A*x=b problem. For example, the 
 eigenvalue problem has a different condition number.

 An intuitive presentation of this topic is given in "Numerical 
 computing with Matlab" by Moler (see the chapter 2 "Linear equations"). 
 A more detailed presentation is done in Golub and Van Loan's "Matrix 
 computations". Other good authors are Demmel "Applied numerical linear 
 algebra", Higham "Accuracy and stability of numerical algorithms" and 
 the pioneering work of Wilkinson "Rounding errors in algebraic 
 processes".

 Best regards,

 Michaël Baudin

 On Mon, 27 Dec 2010 13:15:36 -0800 (PST), "Prof. Dr. Reinaldo Golmia 
 Dante" <tiraduvidascefet at yahoo.com> wrote:
> Hi all,
>
> I understand that Scilab stores the real numbers with foating point 
> numbers,
> that is, with limited precision, and the computed value (answer) is
> not exactly
> equal to 0 (page 23 - manual "Introduction to Scilab").
>  
> Take at look for those examples:
>  
> A = [1 2 3; 4 5 6; 7 8 9]
> A  =
>  
>     1.    2.    3. 
>     4.    5.    6. 
>     7.    8.    9. 
>  
> -->det(A)
>  ans  =
>  
>     6.661D-16  <---------- it should be 0
>  
> -->inv(A)
>  ans  =
>  
>  10^15 *
>  
>   - 4.5035996    9.0071993  - 4.5035996 
>     9.0071993  - 18.014399    9.0071993 
>   - 4.5035996    9.0071993  - 4.5035996            <---------- it
> should appear
> an error message because the matrix A is not invertible (or 
> singular).
>
>  
> -->det(inv(A))
>  ans  =
>  
>     9.007D+15     <-------------- The determinant of invertible
> matrix A^(-1)
> does not exist.
>  
> Other example:
>
> -->B = [1 1; 1 1]
> B  =
>  
>     1.    1. 
>     1.    1. 
>  
> -->det(B)
>  ans  =
>  
>     0.  <-------- it is correct !!
>  
> -->inv(B)
>        !--error 19    <-------- it is correct !!
>
>  
> The previously examples show two integer matrices A and B. The
> determinant of
> matrix A is quite zero, but not,
> and this can propagate an error in case the Scilab developer uses
> that result
> into other future calculations or algorithms.
> The determinant of matrix B is equal to 0 and the answer is correct.
> In case the
> Scilab developer uses that value,
> he or she can use the simple statement for testing like to:
>  if ( det(matrix) <> 0 ) then
> <action 1>                     // The Scilab developer knows that the
> matrix is
> invertible (or nonsingular)
> else
> <action 2>                     // The Scilab developer knows that the
> matrix is
> not invertible (or singular)
> end
>  
> My doubt: "How can I proceed to design any algorithm, which uses
> matrix, if the
> determinant of
>
> the matrix could not be zero and, as the same time, that matrix is 
> not
> invertible ?".
> How can I manage this uncertainty ?
>  
> I appreciate to hearing from you some hints to solve this 
> uncertainty.
>  
> Thank you in advance.
>  
> All best,
> Reinaldo.
>  
> PS: I am not MATLAB user, but has MATLAB got this limited precision 
> for
> determinants ?