Re: Octave QR factorization

[taltman]

Hi Bart,

Depending on the platform (OS-hardware combo) and the algorithm, you
can set a tolerance for floating-point equality. This would be a good
lesson to your students on the difference between the linear algebraic
ideal, and the approximated realization in hardware utilizing a finite
precision method of calculation.

My guess for Matlab's behavior is:

* Like John's SPARC system, you were just lucky with the
  OS/hardware/software version combination.

* Matlab does some behind-the-scenes tweaking of their definition of
  equality on reals depending on the OS/hardware/software combination.

Their assumption is that this will be "good enough" for most
people. But for any real serious numerical work, you should analyze
how much precision you need, and then how much precision you will lose
due to your calculation method. I don't think any algorithm can be
written to predict how much precision you can rely on due to a given
calculation method ( smacks of the Turing halting problem... ).

~Tomer



On Dec 17, 2003 at 2:16pm, John W. Eaton wrote:

jwe >Date: Wed, 17 Dec 2003 14:16:55 -0600
jwe >From: John W. Eaton <address@hidden>
jwe >To: Bart Vandewoestyne <address@hidden>
jwe >Cc: address@hidden
jwe >Subject: Octave QR factorization
jwe >Resent-Date: Wed, 17 Dec 2003 14:16:57 -0600
jwe >Resent-From: address@hidden
jwe >
jwe >On 17-Dec-2003, Bart Vandewoestyne <address@hidden> wrote:
jwe >
jwe >| While I was checking the correctness of an exercise given to my 
jwe >| students during one of my algebra exercise sessions, i found the 
jwe >| following curiosity:
jwe >| 
jwe >| The purpose of the exercise was for students to calculate the QR 
jwe >| factorization of a given matrix A.  Now all I wanted to do is enter 
jwe >| the factorization we found by hand into Octave and see if the product 
jwe >| of Q and R is again the original matrix A.
jwe >| 
jwe >| What I found is that when I multiply the matrices Q and R we've found 
jwe >| by hand and then check for Q*R==A, Octave tells me that Q*R is not 
jwe >| equal to A, most probably due to rounding errors I guess...
jwe >| 
jwe >| If I run the same script in Matlab, Matlab *does* give me that 
jwe >| Q*R=A... so it doesn't experience the rounding errors as much as 
jwe >| Octave does???
jwe >| 
jwe >| So now my question is: why does Matlab gives different results than 
jwe >| Octave?  Is there a big difference in the implementation of QR between 
jwe >| Matlab and Octave?  Is there a way to decrease the rounding errors 
jwe >| made by Octave so we get the same accuracy as Matlab?
jwe >| 
jwe >| The script I used can be found at 
jwe >| http://www.cs.kuleuven.ac.be/~bartv/p370_oef1_6_17.m
jwe >
jwe >You don't provide details about your installation of Octave, but I'm
jwe >going to guess that it is on x86 hardware.  I ran your test on an AMD
jwe >Athlon system and it "failed".  I also ran it on a SPARC system and it
jwe >"succeeded".  I'm fairly certain that the difference is in the way
jwe >that floating point operations are implemented on different machines
jwe >(unfortunately, there are differences even for systems that mostly
jwe >conform to the IEEE standard for floating point arithmetic).
jwe >
jwe >You might find this message
jwe >
jwe >  http://www.octave.org/mailing-lists/help-octave/2003/2185
jwe >
jwe >helpful.
jwe >
jwe >
jwe >jwe
jwe >
jwe >
jwe >
jwe >-------------------------------------------------------------
jwe >Octave is freely available under the terms of the GNU GPL.
jwe >
jwe >Octave's home on the web:  http://www.octave.org
jwe >How to fund new projects:  http://www.octave.org/funding.html
jwe >Subscription information:  http://www.octave.org/archive.html
jwe >-------------------------------------------------------------
jwe >
jwe >



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------