R語言里的矩陣處理學(xué)習(xí)筆記

關(guān)于矩陣，通常都會使用matlab來做處理。其實使用R也可以對矩陣做出一些簡單的處理。而R語言中提供的matrix，matlab包也提供了不少關(guān)于矩陣處理的東西(可以通過??matlab來查看具體函數(shù)）。

一、矩陣的輸入

通常我們使用函數(shù)matrix來創(chuàng)建矩陣，函數(shù)的介紹如下：

matrix(data = NA,nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)

如果想將數(shù)據(jù)按行輸入矩陣，參數(shù)byrow須改為T。

由于矩陣也是特殊的數(shù)組，我們也可以用生成數(shù)組的函數(shù)array（）。具體格式如下：

array(data = NA, dim= length(data), dimnames = NULL)

這里的dim是一個二維數(shù)組，生成的就是矩陣了。

當(dāng)然dim（）函數(shù)，attr()（這個是一個格式轉(zhuǎn)換的函數(shù)）也是可以生成矩陣的。

還有如diag（）可以產(chǎn)生對角矩陣這樣特殊矩陣的函數(shù)。

例如生成下面的這個矩陣（為了便于下面的敘述，這個矩陣記為x，生成命令x<-matrix（1:16,2,8））:

[,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7]  [,8]
[1,]    1    3      5    7    9   11   13   15
[2,]    2    4      6    8  10   12   14   16

我們可以使用如下幾種命令：

matrix（1:16,2,8）

x<-1:16  ;dim（x）<-c(2,8)

array(1:16,c(2,8))

x<-1:16;attr(x,"dim")<-c(2,8)

二、矩陣的基本操作

這里主要的有：矩陣的加法與乘法，矩陣求秩，矩陣的轉(zhuǎn)置，方陣求行列式，矩陣求逆，解線性方程組

1、矩陣的加法與乘法

加法直接使用加號，實現(xiàn)對應(yīng)元素相加。但是要注意兩個矩陣必須可加

矩陣的乘法：和matlab類似，R也給出了兩種乘法：

“*”：對應(yīng)位置元素做乘法，如x*x得到結(jié)果：

> x*x
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    9   25   49   81  121 169  225
[2,]    4   16   36   64  100  144  196 256

"%*%":這個是通常意義下的矩陣乘法，如x%*%t(x)得到結(jié)果：

> x%*%t(x)
     [,1] [,2]
[1,]  680  744
[2,]  744  816

這里乘法也必須是有意義的才行。

通常我們也使用crossprod()函數(shù)來做乘法，crossprod(x,x)效果與x%*%t(x)相同

2、矩陣求秩

這里可以利用qr分解來解決求秩的問題

qr(x)$rank

可以得到x的秩

3、矩陣的轉(zhuǎn)置

常用的命令為t().

R中還可以使用命令aperm()來實現(xiàn)矩陣的廣義轉(zhuǎn)置，函數(shù)用法如下：

aperm(a, perm, ...)

## Default S3 method:

aperm(a, perm = NULL, resize =TRUE, ...)

## S3 method for class 'table'

aperm(a, perm = NULL, resize =TRUE, keep.class = TRUE, ...)

4、方陣求行列式

命令為det()，無須贅述

5、矩陣求逆與解線性方程組

使用函數(shù)solve（）

對于線性方程組b<-A%*%y

的解使用函數(shù)solve（A,b)即可

從而我們知道取b為單位陣時即可求逆，通常簡化為solve（A)

但是值得注意的是，用直接求逆的辦法解線性方程組y<-solve(A)%*%b不僅穩(wěn)定性低，效率也不咋地。

三、矩陣的分解

這里主要介紹矩陣的lu分解，Cholesky分解，特征值與特征向量，QR分解，奇異值分解

1、LU分解

將系數(shù)矩陣A轉(zhuǎn)變成等價兩個矩陣L和U的乘積 ，其中L和U分別是下三角和上三角矩陣。當(dāng)A的所有順序主子式都不為0時，矩陣A可以分解為A=LU，且分解唯一。其中L是單位下三角矩陣，U是上三角矩陣。

library(Matrix) #這里引入函數(shù)包Matrix

> m
     [,1] [,2] [,3]
[1,]    2   -1    3
[2,]    1    2    1
[3,]    2    4    2

> l <- lu(m)

> l
'MatrixFactorization' of Formal class 'denseLU' [package "Matrix"]with 3 slots
  ..@ x   : num [1:9] 2 1 0.5 -1 5 0.5 3 -1 0
  ..@ perm: int [1:3] 1 3 3
  ..@ Dim : int [1:2] 3 3

> LU <- expand(l) #生成P,L,U

> LU
$L
3 x 3 Matrix of class "dtrMatrix" (unitriangular)
     [,1] [,2] [,3]
[1,]  1.0    .    .
[2,]  1.0  1.0    .
[3,]  0.5  0.5  1.0
$U
3 x 3 Matrix of class "dtrMatrix"
     [,1] [,2] [,3]
[1,]    2   -1    3
[2,]    .    5   -1
[3,]    .    .    0
$P（這個矩陣的意思是保持被變換的矩陣第一行不變，二三行對調(diào)）
3 x 3 sparse Matrix of class "pMatrix"
          
[1,] | . .
[2,] . . |
[3,] . | .

可以驗證A = LU$P%*%LU$L%*%LU$U

P為置換矩陣，L為下單位三角矩陣，U為上三角矩陣；

2、Cholesky分解

如果矩陣A為n階對稱正定矩陣，則存在一個非奇異（滿秩）的下三角實矩陣L，使得：A=L%*%t(L)當(dāng)限定的L的對角元素為正時，分解唯一，成為Cholesky分解

> A
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> chol(A)
        [,1]     [,2]     [,3]    [,4]
[1,] 1.414214 0.7071068 0.7071068 0.7071068
[2,] 0.000000 1.2247449 0.4082483 0.4082483
[3,] 0.000000 0.0000000 1.1547005 0.2886751
[4,] 0.000000 0.0000000 0.0000000 1.1180340
> t(chol(A))%*%chol(A)
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> crossprod(chol(A),chol(A))
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
若矩陣為對稱正定矩陣，可以利用Cholesky分解求行列式的值，如：
> prod(diag(chol(A))^2)
[1] 5
> det(A)
[1] 5
若矩陣為對稱正定矩陣，可以利用Cholesky分解求矩陣的逆，這時用函數(shù)chol2inv()，這種用法更有效。

3、特征值與特征向量

函數(shù)eigen（A）用來計算方陣A的特征值與特征向量，返回一個含有特征值與特征向量的列表

> A
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> eigen(A)
$values
[1]  3.620937e+01 -2.209373e+00 -1.050249e-15  8.203417e-16
$vectors
           [,1]        [,2]      [,3]       [,4]
[1,] -0.4140028 -0.82289268  0.4422036 -0.1001707
[2,] -0.4688206 -0.42193991 -0.3487083  0.5349238
[3,] -0.5236384 -0.02098714 -0.6291942 -0.7693354
[4,] -0.5784562  0.37996563  0.5356989  0.3345823
有時只需特征值時，我們使用eigen（A,only.value=T)$value可以快速得到結(jié)果

4、QR分解

A為m×n矩陣可以進行QR分解，A=QR，其中：Q'Q＝I，在R中可以用函數(shù)qr()進行QR分解，例如：
> A=matrix(1:16,4,4) 數(shù)據(jù)分析培訓(xùn)
> qr(A)
$qr
      [,1]     [,2]       [,3]      [,4]
[1,] -5.4772256 -12.7801930 -2.008316e+01 -2.738613e+01
[2,] 0.3651484 -3.2659863 -6.531973e+00 -9.797959e+00
[3,] 0.5477226 -0.3781696 2.641083e-15 2.056562e-15
[4,] 0.7302967 -0.9124744 8.583032e-01 -2.111449e-16
$rank
[1] 2
$qraux
[1] 1.182574e+00 1.156135e+00 1.513143e+00 2.111449e-16
$pivot
[1] 1 2 3 4
attr(,"class")
[1] "qr"
rank項返回矩陣的秩，qr項包含了矩陣Q和R的信息，要得到矩陣Q和R，可以用函數(shù)qr.Q()和qr.R()作用qr()的返回結(jié)果，例如：
> qr.R(qr(A))
      [,1]     [,2]       [,3]      [,4]
[1,] -5.477226 -12.780193 -2.008316e+01 -2.738613e+01
[2,] 0.000000 -3.265986 -6.531973e+00 -9.797959e+00
[3,] 0.000000   0.000000 2.641083e-15 2.056562e-15
[4,] 0.000000   0.000000 0.000000e+00 -2.111449e-16
> qr.Q(qr(A))
      [,1]       [,2]     [,3]    [,4]
[1,] -0.1825742 -8.164966e-01 -0.4000874 -0.37407225
[2,] -0.3651484 -4.082483e-01 0.2546329 0.79697056
[3,] -0.5477226 -8.131516e-19 0.6909965 -0.47172438
[4,] -0.7302967 4.082483e-01 -0.5455419 0.04882607
> qr.Q(qr(A))%*%qr.R(qr(A))
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16

> t(qr.Q(qr(A)))%*%qr.Q(qr(A))
        [,1]       [,2]      [,3]       [,4]
[1,] 1.000000e+00 -1.457168e-16 -6.760001e-17 -7.659550e-17
[2,] -1.457168e-16 1.000000e+00 -4.269046e-17 7.011739e-17
[3,] -6.760001e-17 -4.269046e-17 1.000000e+00 -1.596437e-16
[4,] -7.659550e-17 7.011739e-17 -1.596437e-16 1.000000e+00
> qr.X(qr(A))
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16

5、svd分解

利用函數(shù)svd（）可以對矩陣做svd分解?？匆粋€R提供的例子：

svd> hilbert <- function(n) { i <- 1:n; 1 /outer(i - 1, i, "+") }
svd> X <- hilbert(9)[,1:6]
svd> (s <- svd(X))
$d
[1] 1.668433e+00 2.773727e-01 2.223722e-02 1.084693e-03 3.243788e-05
[6] 5.234864e-07
$u
            [,1]       [,2]       [,3]        [,4]       [,5]        [,6]
 [1,] -0.7244999  0.6265620  0.27350003 -0.08526902 0.02074121 -0.00402455
 [2,] -0.4281556 -0.1298781 -0.64293597  0.55047428 -0.27253421 0.09281592
 [3,] -0.3121985 -0.2803679 -0.33633240 -0.31418014  0.61632113-0.44090375
 [4,] -0.2478932 -0.3141885 -0.06931246 -0.44667149  0.02945426 0.53011986
 [5,] -0.2063780 -0.3140734  0.10786005 -0.30241655 -0.35566839 0.23703838
 [6,] -0.1771408 -0.3026808  0.22105904 -0.09041508 -0.38878613-0.26044927
 [7,] -0.1553452 -0.2877310  0.29280775  0.11551327 -0.19285565-0.42094482
 [8,] -0.1384280 -0.2721599  0.33783778  0.29312535 0.11633231 -0.16079025
 [9,] -0.1248940 -0.2571250  0.36542543  0.43884649 0.46496714  0.43459954
$v
           [,1]       [,2]      [,3]        [,4]       [,5]         [,6]
[1,] -0.7364928  0.6225002  0.2550021 -0.06976287  0.01328234-0.001588146
[2,] -0.4432826 -0.1818705 -0.6866860  0.50860089 -0.19626669  0.041116974
[3,] -0.3274789 -0.3508553 -0.2611139 -0.50473697  0.61605641 -0.259215626
[4,] -0.2626469 -0.3921783  0.1043599 -0.43747940 -0.40833605 0.638901622
[5,] -0.2204199 -0.3945644  0.3509658  0.01612426 -0.46427916-0.675826789
[6,] -0.1904420 -0.3831871  0.5110654  0.53856351  0.44663632 0.257248908
svd> D <- diag(s$d)
svd> s$u %*% D %*% t(s$v) #  X = U D V'
           [,1]      [,2]      [,3]       [,4]       [,5]      [,6]
 [1,] 1.0000000 0.5000000 0.33333333 0.25000000 0.20000000 0.16666667
 [2,] 0.5000000 0.3333333 0.25000000 0.20000000 0.16666667 0.14285714
 [3,] 0.3333333 0.2500000 0.20000000 0.16666667 0.14285714 0.12500000
 [4,] 0.2500000 0.2000000 0.16666667 0.14285714 0.12500000 0.11111111
 [5,] 0.2000000 0.1666667 0.14285714 0.12500000 0.11111111 0.10000000
 [6,] 0.1666667 0.1428571 0.12500000 0.11111111 0.10000000 0.09090909
 [7,] 0.1428571 0.1250000 0.11111111 0.10000000 0.09090909 0.08333333
 [8,] 0.1250000 0.1111111 0.10000000 0.09090909 0.08333333 0.07692308
 [9,] 0.1111111 0.1000000 0.09090909 0.08333333 0.07692308 0.07142857
svd> t(s$u) %*% X %*% s$v #  D = U' X V
              [,1]         [,2]          [,3]         [,4]          [,5]
[1,]  1.668433e+00  2.009230e-16 -2.333610e-16  1.193300e-16 2.347298e-17
[2,]  1.627828e-17  2.773727e-01  7.318365e-19 3.109966e-17 -5.251265e-17
[3,]  1.828617e-17  1.086828e-17  2.223722e-02 4.511721e-18  1.194020e-17
[4,]  2.420517e-17  1.205777e-17  3.433104e-18 1.084693e-03 -4.584063e-18
[5,] -3.406808e-17 -1.147965e-17 -8.968404e-19  6.405788e-18 3.243788e-05
[6,] -1.591696e-17  2.714931e-18  1.721002e-17 -2.358252e-18 1.170640e-17
              [,6]
[1,]  1.015423e-16
[2,]  2.334931e-17
[3,] -1.562373e-17
[4,]  1.364795e-18
[5,] -2.473743e-18
[6,]  5.234864e-07

6、矩陣廣義逆(Moore-Penrose)
  n×m矩陣A+稱為m×n矩陣A的Moore-Penrose逆，如果它滿足下列條件：
①   A A+A=A；②A+A A+= A+；③(A A+)H=A A+；④(A+A)H= A+A
在R的MASS包中的函數(shù)ginv()可計算矩陣A的Moore-Penrose逆，例如：
library(“MASS”)
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
驗證性質(zhì)1：
> A%*%ginv(A)%*%A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
驗證性質(zhì)2：
> ginv(A)%*%A%*%ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
驗證性質(zhì)3:
> t(A%*%ginv(A))
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7
> A%*%ginv(A)
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7
驗證性質(zhì)4:
> t(ginv(A)%*%A)
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7
> ginv(A)%*%A
  [,1] [,2] [,3] [,4]
[1,] 0.7 0.4 0.1 -0.2
[2,] 0.4 0.3 0.2 0.1
[3,] 0.1 0.2 0.3 0.4
[4,] -0.2 0.1 0.4 0.7

對于矩陣，我們還可以使用cbind(),rbind()構(gòu)造分塊矩陣。