Definició i notacions

Calcular la inversa d’una matriu quadrada en R és ben fàcil, fem servir la funció solve(). Més endavant entendreu perquè calcular la inversa és el mateix que “resoldre”.

A <- matrix(c(2, 0, 3, 4, 7, 5,  5, 4, 9), ncol = 3)
A
     [,1] [,2] [,3]
[1,]    2    4    5
[2,]    0    7    4
[3,]    3    5    9
solve(A)
           [,1]        [,2]       [,3]
[1,]  1.4827586 -0.37931034 -0.6551724
[2,]  0.4137931  0.10344828 -0.2758621
[3,] -0.7241379  0.06896552  0.4827586

R fa sempre els càlculs amb la versió decimal (se’n diu de coma flotant, vegeu el final de l’anterior guia) dels nombres implicats. Entre altres coses, això té l’aventatge de ser molt més eficient. I no importa massa la dimensió de la matriu.

A <- matrix(sample(-24:24), nrow = 7) # són els nombres des de -24 fins a 24 reordenats aleatoriament
iA <- solve(A)
round(iA,4)
        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]
[1,]  0.0210  0.0293 -0.0539 -0.0018 -0.0005  0.0395 -0.0518
[2,]  0.0098 -0.0158  0.0178  0.0023 -0.0047 -0.0765  0.0192
[3,] -0.0115 -0.0186 -0.0254  0.0019 -0.0253 -0.0372  0.0147
[4,]  0.0213  0.0115 -0.0369  0.0184 -0.0162  0.0279 -0.0310
[5,]  0.0317  0.0666  0.0806  0.0455  0.0124 -0.0362  0.0286
[6,]  0.0147 -0.0055  0.0523  0.0468  0.0179 -0.0858  0.0675
[7,]  0.0110  0.0442  0.0007  0.0253  0.0115 -0.0207 -0.0379
A %*% iA
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,]  1.000000e+00  2.220446e-16 -9.367507e-17  1.110223e-16  2.775558e-17
[2,]  3.469447e-18  1.000000e+00 -1.279359e-16  1.249001e-16 -5.204170e-17
[3,] -4.857226e-17 -5.551115e-17  1.000000e+00 -2.775558e-17  6.938894e-18
[4,]  1.110223e-16 -1.110223e-16  1.353084e-16  1.000000e+00  5.551115e-17
[5,]  5.551115e-17  1.110223e-16  1.058181e-16 -1.665335e-16  1.000000e+00
[6,]  8.326673e-17  2.220446e-16  5.204170e-18  1.665335e-16 -1.387779e-17
[7,]  1.110223e-16  2.220446e-16 -2.255141e-17  1.110223e-16  0.000000e+00
              [,6]          [,7]
[1,]  0.000000e+00 -4.440892e-16
[2,]  3.122502e-16 -1.387779e-16
[3,] -6.938894e-17 -5.551115e-17
[4,] -5.551115e-17  6.661338e-16
[5,]  4.440892e-16 -8.326673e-17
[6,]  1.000000e+00  1.665335e-16
[7,]  1.665335e-16  1.000000e+00

Com veieu, aixó té l’inconvenient que de vegades no s’entén a la primera el què hem obtingut, aquesta matriu identitat no ho sembla gaire degut als petitissims errors que s’acumulen en els càlculs amb decimals.

Podem demanar el resultat arrodonit i ho veiem clar:

round(A %*% iA,4)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    0    0    0    0    0    0
[2,]    0    1    0    0    0    0    0
[3,]    0    0    1    0    0    0    0
[4,]    0    0    0    1    0    0    0
[5,]    0    0    0    0    1    0    0
[6,]    0    0    0    0    0    1    0
[7,]    0    0    0    0    0    0    1

Resoldre un sistema d’equacions lineals és multiplicar per la inversa.

Com que un sistema d’equacions lineals es pot escriure matricialment com a \({\mathbf A}{\mathbf X}={\mathbf b}\), la solució serà \({\mathbf X}={\mathbf A}^{-1}{\mathbf b}\). Per tant, per resoldre el sistema \[\begin{array}{lcc} x+2y-z&=&-5\\2x-y+z&=&6\\x-y-3z&=&-3 \end{array}\] farem:

A <- matrix(c(1,2,1,2,-1,-1,-1,1,-3), nrow = 3)
A
     [,1] [,2] [,3]
[1,]    1    2   -1
[2,]    2   -1    1
[3,]    1   -1   -3
b <- c(-5, 6, -3)
solve(A,b)
[1]  1 -2  2

Observa que la funció solve() que hem fet servir per calcular la inversa, si li donem un segon argument, ens torna el producte de la inversa per aquest segon argument.

solve(A) %*% b- solve(A,b)  # pràcticament zero, és el mateix
              [,1]
[1,] -2.220446e-16
[2,]  0.000000e+00
[3,]  0.000000e+00

Matrius ortogonals

La matriu de l’exemple de la guia:

P<- matrix(c(1,-1,1,1),nrow = 2) / sqrt(2)
P
           [,1]      [,2]
[1,]  0.7071068 0.7071068
[2,] -0.7071068 0.7071068
solve(P) # la inversa és la transposada
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068
P %*% t(P)
     [,1] [,2]
[1,]    1    0
[2,]    0    1
t(P) %*% P
     [,1] [,2]
[1,]    1    0
[2,]    0    1

Per ortogonalitzar una matriu, R fa servir un mètode anomenat de Gram-Schmidt, disponible al paquet pracma. No hi entrem en aquest curs, però el podeu trobar força ben explicat a (https://ca.wikipedia.org/wiki/Proc%C3%A9s_d%27ortogonalitzaci%C3%B3_de_Gram-Schmidt).

L’apliquem a un exemple.

library(pracma) # potser us cal fer install.packages("pracma") si no el teniu instal·lat

A <- matrix(sample(-12:12), nrow = 5) # són els nombres des de -12 fins a 12 reordenats aleatoriament
A
     [,1] [,2] [,3] [,4] [,5]
[1,]   -7    8   -1  -12    2
[2,]    9    5   -3    6    4
[3,]   -2   -5   -9  -10  -11
[4,]   10    1   -4    3    0
[5,]   11   -6   12   -8    7
nA <- gramSchmidt(A)$Q
nA
           [,1]       [,2]       [,3]        [,4]        [,5]
[1,] -0.3715214  0.5773351  0.1826652 -0.69583768 -0.10535329
[2,]  0.4776704  0.5411482 -0.2364724  0.03352659  0.64957400
[3,] -0.1061490 -0.4467784 -0.6770099 -0.52671119  0.23098560
[4,]  0.5307449  0.2187771 -0.3941799 -0.09768329 -0.71100393
[5,]  0.5838194 -0.3554833  0.5449716 -0.47719932  0.08985195
round(nA %*% t(nA), 4)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1