Обнаружение нуля в решении линейной системы уравнений (Ах = б)

Question

Предположим следующую систему уравнений Ax = b с:

> 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) 

Решение этих уравнений с solve() выходами:

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

Это просто пример, но A и b могут быть любой формы.

Мне нужно определить, равен ли какой-либо компонент x 0. Теперь первые два компонента должны быть равны 0, но они оба выше, чем машина epsilon .Machine$double.eps = 2.22e-16, что делает их очень маленькими, но не равными нулю.

Я думаю, я понимаю, что это вызвано ошибками округления в арифметике с плавающей запятой внутри solve(). Мне нужно знать, можно ли (с практической точки зрения) определить более высокую оценку этих ошибок, поэтому можно обнаружить 0s. Например, вместо

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

можно было бы использовать что-то вроде этого:

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

Давать больше понимания этой проблемы будут оценены.

Answer 1

Один из способов приблизиться к этому - проверить, можно ли найти лучшее решение для линейной системы, если мы предположим, что компоненты равны нулю. Для этого мы хотели бы решить A[3:4]%*%y=b с A%*%c(0,0,x[3],x[4])=A[3:4]%*%c(x[3],x[4]). Это переопределенная система, поэтому мы не можем использовать solve, чтобы найти решение. Однако мы можем использовать qr.solve:

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

Остается проверить, если это решение действительно лучше:

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

Таким образом, у нас есть веские основания подозревать, что x[1]==x[2]==0.

---

В этом простом примере, очевидно, можно угадать истинное решение, глядя на приближенном решении:

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

Это, однако, не очень полезно в общем случае.