Линейная алгебра в R

Матричные операции в 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-разложение позволяет обойти это ограничение.