Матричные операции в R
Для начала, создадим 2 матрицы:
A <- matrix(1:4, 2,2)
A
## [,1] [,2] ## [1,] 1 3 ## [2,] 2 4
B <- matrix(1:6, 3,2)
B
## [,1] [,2] ## [1,] 1 4 ## [2,] 2 5 ## [3,] 3 6
Транспонирование матрицы:
t(A)
## [,1] [,2] ## [1,] 1 2 ## [2,] 3 4
t(B)
## [,1] [,2] [,3] ## [1,] 1 2 3 ## [2,] 4 5 6
Инвертировать можно только квадратные матрицы:
solve(A)
## [,1] [,2] ## [1,] -2 1.5 ## [2,] 1 -0.5
solve(B)
## Error in solve.default(B): 'a' (3 x 2) must be square
Умножение транспонированной матрицы на себя дает квадратную матрицу:
t(B) %*% B
## [,1] [,2] ## [1,] 14 32 ## [2,] 32 77
Но не наоборот:
B %*% t(B)
## [,1] [,2] [,3] ## [1,] 17 22 27 ## [2,] 22 29 36 ## [3,] 27 36 45
Тот же результат что и t(B) %*% B
можно получить проще с помощью crossprod
:
crossprod(B)
## [,1] [,2] ## [1,] 14 32 ## [2,] 32 77
Метод наименьших квадратов в R
Допустим, реальная зависимость между двумя переменными описывается уравнением y = 2 + 3*x
В матричном виде зависимость будет выражать уравнением Y = X * β + ε, где ε это шум или ошибки в измерении y.
true_beta <- c(2,3)
X <- cbind(1, 1:100)
y <- X %*% true_beta + rnorm(100, 0, 10)
plot(X[,2],y)
abline( true_beta[1], true_beta[2], col = 'red')
Матричное решение системы уравнений, минимизируещее сумму квадратов ошибок:
solve( t(X) %*% X ) %*% t(X) %*% y
## [,1] ## [1,] 4.318327 ## [2,] 2.975728
Сравним полученные значения с “истинными” значениями β, которые использовались для генерации данных:
true_beta
## [1] 2 3
Или более компактный способ при помощи crossprod()
:
solve(crossprod(X)) %*% crossprod(X,y)
## [,1] ## [1,] 4.318327 ## [2,] 2.975728
Данный способ является неустойчивым для данных, миниммальные и максимальные значения которых сильно отличаются. Смоделируем такую ситуацию при помощи полинома 2 порядка:
nmin <- 1
nmax <- 1e6
x <- nmin : nmax
X <- cbind(1, x, x^2)
true_beta <- c(1,1,1)
y <- X %*% true_beta + rnorm(nmax - nmin + 1, mean = 0, sd = 100)
plot(x,y, 'l')
Попробуем найти решение LS “обычным” способом:
solve(crossprod(X)) %*% crossprod(X,y)
## Error in solve.default(crossprod(X)): system is computationally singular: reciprocal condition number = 5.55549e-25
Более стабильным будет являться решение с помощью т.н. QR-декомпозиции. QR-декомпозиция разбивает матрицу X на две матрицы Q и R:
QR <- qr(X)
Q <- qr.Q(QR)
R <- qr.R(QR)
Проверим, действительно ли произведение Q и R матриц дает матрицу X
(Q %*% R)[10:15,]
## x ## [1,] 1 10 99.99203 ## [2,] 1 11 120.99200 ## [3,] 1 12 143.99200 ## [4,] 1 13 168.99203 ## [5,] 1 14 195.99200 ## [6,] 1 15 224.99200
X[10:15,]
## x ## [1,] 1 10 100 ## [2,] 1 11 121 ## [3,] 1 12 144 ## [4,] 1 13 169 ## [5,] 1 14 196 ## [6,] 1 15 225
Матрица Q обладает особым свойством: t(Q)*Q = I
t(Q) %*% Q
## [,1] [,2] [,3] ## [1,] 1.000000e+00 1.149254e-16 4.599619e-15 ## [2,] 1.149254e-16 1.000000e+00 -1.307418e-14 ## [3,] 4.599619e-15 -1.307418e-14 1.000000e+00
Матрица R также обладает особым свойством: это треугольная матрица c нулями ниже диагонали:
R
## x ## [1,] -1000 -500000500 -3.333338e+14 ## [2,] 0 288675135 2.886754e+14 ## [3,] 0 0 -7.453560e+13
Решение LS c помощью QR-разложения:
solve.qr(QR, y)
## [,1] ## 1.103111 ## x 0.999999 ## 1.000000
Напомним, что данная система уравнений не имела решения при помощи стандартного подхода нахождения обратной матрицы с помощью solve()
, т.к. численные методы реализованные для данной функции не позволяют найти обратную матрицу в данном случае:
solve( t(X) %*% X)
## Error in solve.default(t(X) %*% X): system is computationally singular: reciprocal condition number = 5.55549e-25
QR-разложение позволяет обойти это ограничение.