[racket] matrix-solve and approximation errors

2014-04-16 Thread Laurent
I've just been bitten by a bad case of floating-point error with `matrix-solve` (and a bad CPU that has some floating-point issues): (let ([n 0.8201 #;(+ (* .9 .9)(* .1 .1))]) (matrix-solve (matrix [[ 1 0 .9 1] [ 0 1 .1 1] [.9 .1 n 1] [ 1 1

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Laurent
Forgot to mention that with the true value of n=0.82, this of course returns the correct solution: (let ([n 0.82 #;(+ (* .9 .9)(* .1 .1))]) (matrix-solve (matrix [[ 1 0 .9 1] [ 0 1 .1 1] [.9 .1 n 1] [ 1 1 1 0]]) (col-matrix [0 0 0 1]))) ; -> (array #

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Jens Axel Søgaard
I forgot to note, that I get the same results as you, so I don't think the CPU is to blame. 2014-04-16 12:18 GMT+02:00 Jens Axel Søgaard : > Hi Laurent, > > I think the underlying problem is that the matrix is *very* close to > an invertible one: > >> (matrix-determinant > (matrix [[ 1 0

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Jens Axel Søgaard
Hi Laurent, I think the underlying problem is that the matrix is *very* close to an invertible one: > (matrix-determinant (matrix [[ 1 0 9/10 1] [ 0 1 1/10 1] [ 9/10 1/10 82/100 1] [ 1 1 1 0]])) 0 /Jens Axel

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Jens Axel Søgaard
2014-04-16 12:18 GMT+02:00 Jens Axel Søgaard : > Hi Laurent, > > I think the underlying problem is that the matrix is *very* close to > an invertible one: non-invertible >> (matrix-determinant > (matrix [[ 1 0 9/10 1] > [ 0 1 1/10 1] > [ 9

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Laurent
Interesting; even more when using rounded numbers returns a correct solution! To protect against such errors, maybe `matrix-solve` could run a post-check to verify that M×X=B up to some epsilon, depending on an optional argument? Or maybe just mention to do this check systematically in the docs?

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Neil Toronto
FWIW, Jens Axel's "very close" needs quantifying: the entries are less than 2^-52 away from comprising a non-invertible matrix. :) But approximation error is a big deal, and quantifying and responding to it currently has little support in `math/matrix'. We need two things to start with. One re

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Konrad Hinsen
Neil Toronto writes: > The design decision is to determine how much the solvers should automate > detecting and responding to floating-point error. I have no idea what > extent popular linear algebra libraries go to, nor whether we should try > to do better. Any linear-algebra library mean

Re: [racket] matrix-solve and approximation errors

2014-04-16 Thread Laurent
Sounds like a good idea, but possibly cumbersome to program. Or maybe just define a error-checking-level parameter that can be adjusted from `unsafe` to `noob`? On Wed, Apr 16, 2014 at 9:18 PM, Konrad Hinsen wrote: > Neil Toronto writes: > > > The design decision is to determine how much the s

Re: [racket] matrix-solve and approximation errors

2014-04-17 Thread Konrad Hinsen
Laurent writes: > Sounds like a good idea, but possibly cumbersome to program. > Or maybe just define a error-checking-level parameter that can be adjusted > from > `unsafe` to `noob`? I'd expect this more difficult to design than to program. My first idea would be to base everything on sin

Re: [racket] matrix-solve and approximation errors

2014-04-17 Thread Matthias Felleisen
Neil, a two-layer design, using option contracts, could be a solution here. On Apr 16, 2014, at 3:18 PM, Konrad Hinsen wrote: > Any linear-algebra library meant for "serious" work does not do any > non-essential work

Re: [racket] matrix-solve and approximation errors

2014-04-20 Thread Neil Toronto
On 04/17/2014 10:07 AM, Matthias Felleisen wrote: Neil, a two-layer design, using option contracts, could be a solution here. Interesting! So we could have both math/matrix and (say) math/unchecked-matrix, where the former enables the option contracts (e.g. that check matrix condition number

Re: [racket] matrix-solve and approximation errors

2014-04-20 Thread Matthias Felleisen
This makes a lot of sense. I checked this out a few years back when I wanted to formulate a contract for adaptive integration as an example for a paper. I couldn't find much and gave up. I should have realized that this is a research problem, possibly a dissertation problem. -- Matthias O