0

Suppose the following system of equations Ax = b with:

> A <- matrix(c(2,0,-1,0,0,2,2,1,-1,2,0,0,0,1,0,0), ncol = 4)
> A
     [,1] [,2] [,3] [,4]
[1,]    2    0   -1    0
[2,]    0    2    2    1
[3,]   -1    2    0    0
[4,]    0    1    0    0

> b <- c(-2,5,0,0)

Solving these equations with solve() yields:

> x <- solve(A,b)
> x
[1] 6.66e-16 4.44e-16 2.00e+00 1.00e+00

This is just an example, but A and b can be of any form.

I need to detect whether any component of x is 0. Now, the first two components should actually be 0, but they are both higher than the machine epsilon .Machine$double.eps = 2.22e-16 which makes them very small, but not equal to zero.

I think I understand that this is caused by rounding errors in floating point arithmetic inside solve(). What I need to know is whether it is possible (from a practical point of view) to determine the higher bound of these errors, so 0s can be detected. For example, instead of 

> x == 0
[1] FALSE FALSE FALSE FALSE

one would use something like this:

> x > -1e-15 & x < 1e-15
[1]  TRUE  TRUE FALSE FALSE

Giving more insight into this problem would be appreciated.

Davor Josipovic
  • 5,296
  • 1
  • 39
  • 57
  • can't you just use `round()` ? – mtoto Jan 20 '17 at 10:25
  • Up to what decimal? I could use so many quick & dirty things, but before I use any, I would like some (statistical) justification for it. So this question might have many answers, one better than the other... You can certainly start with that one :) – Davor Josipovic Jan 20 '17 at 10:29
  • 1
    Did someone say statistics? If that's the case, then think about what you would consider zero. Anything below/above that is effectively zero. – Roman Luštrik Jan 20 '17 at 10:40
  • @RomanLuštrik, my own 'statistical - misnomer here' approach was under assumption of `A %*% x == b` to calculate the worst case boundaries of `x`. If 0 falls within, then it is 0. Problem is that in this case `A %*% x == b` yields `[1] TRUE TRUE FALSE FALSE`, and thus can not be relied upon... That is why I was hoping someone could shed more light into workings of `solve()` and what the reasonable assumptions could be. – Davor Josipovic Jan 20 '17 at 10:49

1 Answers1

1

One way to approach this is to check if we can find a better solution to the linear system if we assume the components to be zero. For that we would want to solve A[3:4]%*%y=b since A%*%c(0,0,x[3],x[4])=A[3:4]%*%c(x[3],x[4]). This is an overdetermined system so we can't use solve to find a solution. We can however use qr.solve:

> x.new = c(0,0,qr.solve(A[,3:4],b))

It remains to check if this solution is really better:

> norm(A%*%x.new - b) < norm(A%*%x - b)
[1] TRUE

Thus we have a good reason to suspect that x[1]==x[2]==0.

In this simple example it is obviously possible to guess the true solution by looking at the approximate solution:

> x.true = c(0,0,2,1)
> norm(A%*%x.true - b)
[1] 0

This is however not very helpful in the general case.

Jonathan von Schroeder
  • 1,683
  • 11
  • 14