Решение неквадратной линейной системы с помощью R
Как решить неквадратную линейную систему с R: A X = B?
(в случае, если система не имеет решения или бесконечно много решений)
Пример:
A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
B=matrix(c(-17,28,11),3,1,T)
A
[,1] [,2] [,3] [,4]
[1,] 0 1 -2 3
[2,] 5 -3 1 -2
[3,] 5 -2 -1 1
B
[,1]
[1,] -17
[2,] 28
[3,] 11
2 ответа
Если матрица A имеет больше строк, чем столбцов, то вам следует использовать метод наименьших квадратов.
Если матрица A имеет меньше строк, чем столбцов, то вам следует выполнить разложение по сингулярным значениям. Каждый алгоритм делает все возможное, чтобы дать вам решение, используя предположения.
Вот ссылка, которая показывает, как использовать SVD в качестве решателя:
http://www.ecse.rpi.edu/~qji/CV/svd_review.pdf
Давайте применим его к вашей проблеме и посмотрим, работает ли он:
Ваша матрица ввода A и известный вектор RHS B:
> A=matrix(c(0,1,-2,3,5,-3,1,-2,5,-2,-1,1),3,4,T)
> B=matrix(c(-17,28,11),3,1,T)
> A
[,1] [,2] [,3] [,4]
[1,] 0 1 -2 3
[2,] 5 -3 1 -2
[3,] 5 -2 -1 1
> B
[,1]
[1,] -17
[2,] 28
[3,] 11
Давайте разложим ваш A матрица:
> asvd = svd(A)
> asvd
$d
[1] 8.007081e+00 4.459446e+00 4.022656e-16
$u
[,1] [,2] [,3]
[1,] -0.1295469 -0.8061540 0.5773503
[2,] 0.7629233 0.2908861 0.5773503
[3,] 0.6333764 -0.5152679 -0.5773503
$v
[,1] [,2] [,3]
[1,] 0.87191556 -0.2515803 -0.1764323
[2,] -0.46022634 -0.1453716 -0.4694190
[3,] 0.04853711 0.5423235 0.6394484
[4,] -0.15999723 -0.7883272 0.5827720
> adiag = diag(1/asvd$d)
> adiag
[,1] [,2] [,3]
[1,] 0.1248895 0.0000000 0.00000e+00
[2,] 0.0000000 0.2242431 0.00000e+00
[3,] 0.0000000 0.0000000 2.48592e+15
Вот ключ: третье собственное значение в d очень маленький; наоборот, диагональный элемент в adiag очень большой Перед решением установите его равным нулю:
> adiag[3,3] = 0
> adiag
[,1] [,2] [,3]
[1,] 0.1248895 0.0000000 0
[2,] 0.0000000 0.2242431 0
[3,] 0.0000000 0.0000000 0
Теперь давайте вычислим решение (см. Слайд 16 по ссылке, которую я дал вам выше):
> solution = asvd$v %*% adiag %*% t(asvd$u) %*% B
> solution
[,1]
[1,] 2.411765
[2,] -2.282353
[3,] 2.152941
[4,] -3.470588
Теперь, когда у нас есть решение, давайте заменим его обратно, чтобы увидеть, дает ли оно то же самое B:
> check = A %*% solution
> check
[,1]
[1,] -17
[2,] 28
[3,] 11
Это B сторона, с которой вы начали, так что я думаю, что мы хороши.
Вот еще одно приятное обсуждение SVD от AMS:
Цель состоит в том, чтобы решить Ax = b
где A - это p по q, x - это q по 1, а b - это p по 1 для x, заданного для A и b.
Подход 1: Обобщенный обратный: Мур-Пенроуз https://en.wikipedia.org/wiki/Generalized_inverse
Умножая обе части уравнения, получим
A'Ax = A' b
где A ' - транспонирование A. Обратите внимание, что теперь A'A - это матрица q на q. Один из способов решить эту проблему теперь умножить обе части уравнения на обратную величину A'A. Который дает,
x = (A'A) ^ {- 1} A 'b
Это теория за обобщенным обратным. Здесь G = (A'A)^{-1} A' является псевдообратным к A.
library(MASS)
ginv(A) %*% B
# [,1]
#[1,] 2.411765
#[2,] -2.282353
#[3,] 2.152941
#[4,] -3.470588
Подход 2: Обобщенный обратный с использованием SVD.
@duffymo использовала SVD для получения псевдообратного значения A.
Тем не менее, последние элементы svd(A)$d может быть не такой маленький, как в этом примере. Таким образом, вероятно, не следует использовать этот метод как есть. Вот пример, где ни один из последних элементов A не близок к нулю.
A <- as.matrix(iris[11:13, -5])
A
# Sepal.Length Sepal.Width Petal.Length Petal.Width
# 11 5.4 3.7 1.5 0.2
# 12 4.8 3.4 1.6 0.2
# 13 4.8 3.0 1.4 0.1
svd(A)$d
# [1] 10.7820526 0.2630862 0.1677126
Одним из вариантов будет выглядеть как единственные значения в cor(A)
svd(cor(A))$d
# [1] 2.904194e+00 1.095806e+00 1.876146e-16 1.155796e-17
Теперь ясно, что присутствуют только два больших сингулярных значения. Итак, теперь можно применить svd к A, чтобы получить псевдообратную, как показано ниже.
svda <- svd(A)
G = svda$v[, 1:2] %*% diag(1/svda$d[1:2]) %*% t(svda$u[, 1:2])
# to get x
G %*% B