Как решить неквадратную линейную систему с 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
Если в матрице 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:
http://www.ams.org/samplings/feature-column/fcarc-svd
asvd$v %*% adiag %*% t(asvd$u), которую вы строите, является псевдообратной матрицей A, не так ли?
- person Basj; 26.10.2018
Цель состоит в том, чтобы решить Ax = b
где A равно p на q, x равно q на 1 и b равно p на 1 для x при заданном A и б.
Подход 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
MASS::ginv демонстрационный подход. Когда-нибудь я отредактирую ответ с другой реализацией. Ваше здоровье.
- person Suren; 03.09.2019
solveупоминается, что qr.solve может обрабатывать неквадратные системы. Поэтому, когда я использую приведенные здесь A и B и пытаюсь выполнитьqr.solve(A,B), я получаю сообщение об ошибкеError in qr.solve(A, B) : singular matrix 'a' in solve. Есть идеи? - person Gimelist schedule 21.08.2014