Вычислить обратную неквадратную матрицу в R

Я новичок в языке R и пытаюсь выяснить, как можно вычислить обратную матрицу, которая не является квадратной. (не квадратный? неправильный? Я не уверен в правильной терминологии).

Из моей книги и быстрого поиска в Google (см. источник) я нашел вас может использовать solve(a), чтобы найти обратную матрицу, если a квадратная.

Матрица, которую я создал, насколько я понимаю, не квадратная:

  > matX<-matrix(c(rep(1, 8),2,3,4,0,6,4,3,7,-2,-4,3,5,7,8,9,11), nrow=8, ncol=3);
  > matX

       [,1] [,2] [,3]
  [1,]    1    2   -2
  [2,]    1    3   -4
  [3,]    1    4    3
  [4,]    1    0    5
  [5,]    1    6    7
  [6,]    1    4    8
  [7,]    1    3    9
  [8,]    1    7   11
  >

Есть ли функция для решения матрицы такого размера или нужно что-то делать с каждым элементом? поскольку функция solve() выдает эту ошибку:

  Error in solve.default(matX) : 'a' (8 x 3) must be square

Расчет, который я пытаюсь получить из приведенной выше матрицы: (matX'matX)^-1

Заранее спасибо.


r matrix matrix-inverse
person Reanimation    schedule 26.01.2014    source источник
comment
Ссылаясь на вашу последнюю строку, попробуйте: solve(t(matX) %*% matX)   -  person Mark Heckmann    schedule 26.01.2014

Ответы (2)


arrow_upward
6
arrow_downward

ginv ginv в пакете MASS даст обобщенную обратную матрицу. Предварительное умножение исходной матрицы на нее даст тождество:

library(MASS)
inv <- ginv(matX)

# test it out
inv %*% matX
##               [,1]         [,2]          [,3]
## [1,]  1.000000e+00 6.661338e-16  4.440892e-15
## [2,] -8.326673e-17 1.000000e+00 -1.110223e-15
## [3,]  6.938894e-17 8.326673e-17  1.000000e+00

Как было предложено в комментариях, это можно отобразить лучше, используя zapsmall:

zapsmall(inv %*% matX)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Инверсия matX'matX теперь:

tcrossprod(inv)
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

решить однако, если вы хотите вычислить обратное число matX'matX, вам это не нужно. Это не касается MASS:

solve(crossprod(matX))
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

svd Можно также использовать svd, и он также не требует MASS:

with(svd(matX), v %*% diag(1/d^2) %*% t(v))
##              [,1]         [,2]         [,3]
## [1,]  0.513763935 -0.104219636 -0.002371406
## [2,] -0.104219636  0.038700372 -0.007798748
## [3,] -0.002371406 -0.007798748  0.006625269

ДОБАВЛЕНО немного дополнительной информации.

person G. Grothendieck    schedule 26.01.2014
comment
zapsmall может быть удобно для более четкого отображения результатов inv %*% matX ... - person Ben Bolker; 26.01.2014
comment
О, круто :D Спасибо, что нашли время ответить :D Это здорово. - person Reanimation; 26.01.2014
comment
Для решения svd зачем квадратировать диагональную матрицу? диаг.(1/д^2) - person user1143669; 21.02.2020
comment
@ user1143669, пожалуйста, прочитайте предпоследнюю строку вопроса. - person G. Grothendieck; 21.02.2020

arrow_upward
2
arrow_downward

Вы можете сделать так называемую «псевдоинверсию Мура-Пенроуза». Вот функция exp.mat, которая сделает это за вас. Также есть пример использования здесь.

exp.mat():

#The exp.mat function performs can calculate the pseudoinverse of a matrix (EXP=-1)
#and other exponents of matrices, such as square roots (EXP=0.5) or square root of 
#its inverse (EXP=-0.5). 
#The function arguments are a matrix (MAT), an exponent (EXP), and a tolerance
#level for non-zero singular values.
exp.mat<-function(MAT, EXP, tol=NULL){
    MAT <- as.matrix(MAT)
    matdim <- dim(MAT)
    if(is.null(tol)){
        tol=min(1e-7, .Machine$double.eps*max(matdim)*max(MAT))
    }
    if(matdim[1]>=matdim[2]){ 
        svd1 <- svd(MAT)
        keep <- which(svd1$d > tol)
        res <- t(svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep]))
    }
    if(matdim[1]<matdim[2]){ 
        svd1 <- svd(t(MAT))
        keep <- which(svd1$d > tol)
        res <- svd1$u[,keep]%*%diag(svd1$d[keep]^EXP, nrow=length(keep))%*%t(svd1$v[,keep])
    }
    return(res)
}

пример использования:

source("exp.mat.R")
X <- matrix(c(rep(1, 8),2,3,4,0,6,4,3,7,-2,-4,3,5,7,8,9,11), nrow=8, ncol=3)
iX <- exp.mat(X, -1)
zapsmall(iX %*% X) # results in identity matrix
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
person Marc in the box    schedule 26.01.2014
comment
Привет. Спасибо за ваше время. Есть ли какой-то конкретный пакет, который мне нужно использовать для использования функции, которую вы указали выше? - person Reanimation; 26.01.2014
comment
О, подождите, я думаю, вам нужно создать эту функцию самостоятельно в вашей программе? - person Reanimation; 26.01.2014
comment
@Reanimation - Добро пожаловать. Нет, вы можете просто загрузить функцию как есть (т. е. передать код функции в вашу среду R). Если вы предпочитаете загружать пакет, а не sourceэту функцию, аналогичная функция (pseudoinverse) также может быть найдена в пакете corpcor. Преимущество функции exp.mat в том, что она также может принимать другие экспоненциальные значения (например, EXP=-0.5). - person Marc in the box; 26.01.2014
comment
MASS::ginv также реализует псевдоинверсию и поставляется со стандартной установкой R... - person Ben Bolker; 26.01.2014
comment
@BenBolker - Ура. Не могу поверить, что не знал о ginv — спасибо, что указали на это. - person Marc in the box; 26.01.2014