[racket] matrix-solve and approximation errors
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 requires engineering, the other a design decision.
The engineering is implementing a function that computes condition
numbers, which you can use to estimate error in the results of
`matrix-solve' and other inverse-like functions. IIRC, it requires
computing eigendecomposition or singular value decomposition, both of
which are useful for other things (e.g. principal component analysis,
different solving strategies) and are on the radar.
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.
Neil ⊥
On 04/16/2014 05:34 AM, Laurent wrote:
> 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?
>
>
> On Wed, Apr 16, 2014 at 1:31 PM, Jens Axel Søgaard
> <jensaxel at soegaard.net <mailto:jensaxel at soegaard.net>> wrote:
>
> 2014-04-16 12:18 GMT+02:00 Jens Axel Søgaard <jensaxel at soegaard.net
> <mailto:jensaxel at soegaard.net>>:
> > 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/10 1/10 82/100 1]
> > [ 1 1 1 0]]))
> > 0
>
>
>
>
> ____________________
> Racket Users list:
> http://lists.racket-lang.org/users
>