R语言系统教程（三）：多维数组和矩阵

最新推荐文章于 2025-06-19 10:56:06 发布

原创最新推荐文章于 2025-06-19 10:56:06 发布 · 1.2w 阅读

39 ·

CC 4.0 BY-SA版权

文章标签：

#r语言 #统计学

R语言系统教程专栏收录该内容

10 篇文章

订阅专栏

本文详细介绍了R语言中数组和矩阵的操作，包括如何生成数组和矩阵，数组下标的操作，矩阵的转置、行列式、内积、外积、乘法等运算，以及解决线性方程组、求矩阵的逆、特征值和奇异值分解。此外，还讨论了与矩阵运算相关的函数，如apply()、cbind()和rbind()等。

3.1 生成数组或矩阵

数组可以看成是带多个下标的类型相同的元素的集合，但是注意，尽管数学上向量可以看成一维数组，但是R中不行，因为数组一定是有维度属性的，但向量没有。

3.1.1 将向量定义为数组

向量只有定义了维数(dim)属性后才能看做数组。

z = 1:12
dim(z)
## NULL
dim(z) = c(3,4)
z
##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12
dim(z) = 12
z
## [1] 1 2 3 4 5 6 7 8 9 10 11 12

3.1.2 用array()函数构造多维数组

array()函数声明如下：
array(data = NA, dim = length(data), dimnames = NULL)
其中data为向量数据，dim是维度向量，默认为原向量的长度，dimnames为各维度的名字。注意如果data数据个数不足或多于构造多维数组所需的数量，则采用循环或只取较前部分的方式。因此使用单0或单1数据构造全0或全1矩阵的方法非常常见。

X = array(1:20, dim = c(4, 5))
X
##     [,1] [,2] [,3] [,4] [,5]
## [1,]   1    5    9   13   17
## [2,]   2    6   10   14   18
## [3,]   3    7   11   15   19
## [4,]   4    8   12   16   20
Z = array(0, dim = c(3, 4, 2))
Z
## , , 1
##
##     [,1] [,2] [,3] [,4]
## [1,]   0    0    0    0
## [2,]   0    0    0    0
## [3,]   0    0    0    0
##
## , , 2
##
##     [,1] [,2] [,3] [,4]
## [1,]   0    0    0    0
## [2,]   0    0    0    0
## [3,]   0    0    0    0

上面这种全0的定义方式常用来对数组进行初始化。

3.1.3 用matrix()函数构造矩阵

矩阵即是二维数组，是R中用的最多的数据格式之一，函数matrix()是构造矩阵的函数，其构造形式如下：
matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
其中data是向量数据，nrow是矩阵行数，ncol是矩阵列数，当byrow=TRUE时，数据按行存放，否则按列存放，dimnames为维数名字。
注意：若仅指定nrow或ncol时，如果data数据个数不够，则另一个取1，循环补齐，若为nrow或ncol的倍数时，则未指定的自动取其另一个因数，比其大却不是倍数，warning！（没看懂没关系，看例子你就懂了(╹▽╹)）

A = matrix(1:15, nrow = 3, ncol = 5, byrow = TRUE)
A
##     [,1] [,2] [,3] [,4] [,5]
## [1,]   1    2    3    4    5
## [2,]   6    7    8    9   10
## [3,]  11   12   13   14   15
matrix(1:15, nrow=3,byrow=TRUE)
##     [,1] [,2] [,3] [,4] [,5]
## [1,]   1    2    3    4    5
## [2,]   6    7    8    9   10
## [3,]  11   12   13   14   15
matrix(1:15, ncol=5,byrow=TRUE)
##     [,1] [,2] [,3] [,4] [,5]
## [1,]   1    2    3    4    5
## [2,]   6    7    8    9   10
## [3,]  11   12   13   14   15

3.2 数组下标

与向量相同，数组同样可以按照下标访问或修改元素。

3.2.1 数组下标

使用方括号形式，并用逗号分隔开不同维度的下标即可

a = 1:24
dim(a) = c(2,3,4)
a[2,1,2]
## [1] 8

下标同样为向量形式，如果不写则表示此维全取

a[1, 2:3, 2:3]
##     [,1] [,2]
## [1,]   9   15
## [2,]  11   17
a[1, , ]
##     [,1] [,2] [,3] [,4]
## [1,]   1    7   13   19
## [2,]   3    9   15   21
## [3,]   5   11   17   23
a[, 2, ]
##     [,1] [,2] [,3] [,4]
## [1,]   3    9   15   21
## [2,]   4   10   16   22
a[1, 1, ]
## [1] 1 7 13 19

还有一种是通过向量下标方式取值，但要注意这样取出来的是向量，而不是数组。
注：由于格式问题，从这里开始作者更换了代码输出样式，>后代表的是代码，其他的是输出。

> b = a[3:6]
> b
[1] 3 4 5 6

3.2.2 不规则的数组下标

如果想把数组中的任意位置元素作为数组访问，其方法是使用一个二维数组作为数组下标，二维数组的每一行是一个元素的下标，列数为数组维数。

> b = matrix(c(1,1,1,2,2,3,1,3,4,2,1,4), ncol = 3, byrow = T)
> b
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    3
[3,]    1    3    4
[4,]    2    1    4
> a[b]
[1]  1 16 23 20
> dim(a[b])
NULL

但是注意取出来的仍然是向量。

3.3 数组的四则运算

可以对数组之间进行四则运算，与向量相似，这时进行的是数组对应元素的四则运算，一般要求数组规模相同，例如

> A = matrix(1:6, nrow = 2, byrow = T)
> A
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
> B = matrix(1:6, nrow = 2)
> B
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> C = matrix(c(1,2,2,3,3,4), nrow = 2)
> C
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    4
> D = 2*C + A/B
> D
     [,1]     [,2] [,3]
[1,]    3 4.666667  6.6
[2,]    6 7.250000  9.0

对于形状不一致的向量或数组也可以进行四则运算，一般规则是将向量或数组中的数据与对应向量或数组中的数据进行运算，把短向量或数组中的数据循环使用，从而可以与长向量进行匹配，并尽可能的保留共同属性，如

> x1 = c(100, 200)
> x2 = 1:6
> x1 + x2
[1] 101 202 103 204 105 206
> x3 = matrix(1:6, nrow = 3)
> x1 + x3
     [,1] [,2]
[1,]  101  204
[2,]  202  105
[3,]  103  206
> x2 = 1:5
> x1 + x2
[1] 101 202 103 204 105
Warning message:
In x1 + x2 :
  longer object length is not a multiple of shorter object length

可以看到当向量与数组匹配时，向量按列匹配，当两个数组不匹配时，则会提出警告。
但对于不匹配的数组的四则运算，则会报错

> x1 = matrix(1:6, nrow = 2)
> x2 = matrix(1:6, nrow = 3)
> x1+x2
Error in x1 + x2 : non-conformable arrays

3.4 矩阵的运算

3.4.1 转置运算

> A = matrix(1:6, nrow = 2)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> t(A)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6

3.4.2 行列式

> det(matrix(1:4, ncol = 2))
[1] -2

3.4.3 向量内积

%*%应用于两向量时表示求二者内积

> x = 1:5
> y = 2*x
> x %*% y
     [,1]
[1,]  110

当%*%应用于向量和矩阵时，是将向量与矩阵按列求内积，并不是数学直观上的将向量看做列矩阵运算，因此为了避免带来一些没有必要的问题，建议在涉及到向量与矩阵间运算时，一律事先使用matrix函数将向量转为矩阵再进行运算。
还有一种函数写法，即crossprod()函数，另外tcrossprod()可以用来求外积。

> crossprod(x,y)
     [,1]
[1,]  110
> crossprod(x)
     [,1]
[1,]   55

3.4.4 向量外积

%o%可以用来求向量外积

> x %o% y
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    4    6    8   10
[2,]    4    8   12   16   20
[3,]    6   12   18   24   30
[4,]    8   16   24   32   40
[5,]   10   20   30   40   50

还有outer()函数也可以实现，声明如下
outer(X, Y, fun = “*”,…)
X和Y求外积的向量，fun则是外积运算函数，默认为乘法，outer()函数可以用来很方便的生成自定义网格。

> x = 1:5
> y = 1:5
> outer(x,y)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    2    4    6    8   10
[3,]    3    6    9   12   15
[4,]    4    8   12   16   20
[5,]    5   10   15   20   25
> outer(x,y,FUN = "+")
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    3    4    5    6
[2,]    3    4    5    6    7
[3,]    4    5    6    7    8
[4,]    5    6    7    8    9
[5,]    6    7    8    9   10
> S = function(X, Y){
+   return(X+2*Y)
+ }
> outer(x,y,FUN = S)
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    5    7    9   11
[2,]    4    6    8   10   12
[3,]    5    7    9   11   13
[4,]    6    8   10   12   14
[5,]    7    9   11   13   15

3.4.5 矩阵乘法

%*%应用于两矩阵时表示矩阵乘法

> A = array(1:9, dim = (c(3,3)))
> A
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> B = array(9:1, dim = (c(3,3)))
> B
     [,1] [,2] [,3]
[1,]    9    6    3
[2,]    8    5    2
[3,]    7    4    1
> C = A*B
> C
     [,1] [,2] [,3]
[1,]    9   24   21
[2,]   16   25   16
[3,]   21   24    9
> D = A %*% B
> D
     [,1] [,2] [,3]
[1,]   90   54   18
[2,]  114   69   24
[3,]  138   84   30

另外，根据其数学意义，可以知道crossprod(A,B)表示的是 t(A) %%B, 而 tcrossprod(A,B)表示的是 A%%t(B)

3.4.6 生成对角阵和矩阵取对角运算

均是使用diag(v)函数，但当v为向量时，diag(v)生成一个以v为对角的方阵，而当v为矩阵时，diag(v)取其对角元素。
下面的例子除了一般情况外，还考察了v为一维数组和二维单列数组的情况，大家可以从中看到区别。

> v = c(1,4,5)
> diag(v)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    4    0
[3,]    0    0    5
> M = array(1:9,dim = c(3,3))
> diag(M)
[1] 1 5 9
> dim(v) = 3
> diag(v)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    4    0
[3,]    0    0    5
> dim(v) = c(1,3)
> diag(v)
[1] 1
> dim(v) = c(1,3,1)
> diag(v)
Error in diag(v) : 'x'阵列但却不是一维的

3.4.7 解线性方程组和求矩阵的逆矩阵

使用solve(A)函数即可求矩阵A的逆，而solve(A，b)可以求线性方程组 $A x = b$ 的解

> A = t(array(c(1:8, 10), dim = c(3,3)))
> b = c(1,1,1)
> x = solve(A,b)
> x
[1] -1.000000e+00  1.000000e+00  3.330669e-16 # =0
> B = solve(A)
> B
           [,1]      [,2] [,3]
[1,] -0.6666667 -1.333333    1
[2,] -0.6666667  3.666667   -2
[3,]  1.0000000 -2.000000    1
> B %*% matrix(b, nrow = 3)
              [,1]
[1,] -1.000000e+00
[2,]  1.000000e+00
[3,]  4.440892e-16 # =0

可以看到两种求法结果是几乎一样的，不完全一样的原因是因为计算机不同的数值解法存在舍入误差，完全可以忽略。另外solve()函数只适用于A为方阵的情况，非方阵会报错。

> A = t(array(c(1:5, 7), dim = c(2,3)))
> b = c(1,1)
> solve(A,b)
Error in solve.default(A, b) : 'a' (3 x 2) must be square

3.4.8 求矩阵的特征值与特征向量

使用eigen()函数求对称矩阵特征值与特征向量，结果由列表形式给出，列表的概念见第四期。

> A = t(array(c(1:8, 10), dim = c(3,3)))
> Sm<-crossprod(A,A)
> ev <- eigen(Sm)
> ev
eigen() decomposition
$values
[1] 303.19533618   0.76590739   0.03875643

$vectors
           [,1]         [,2]       [,3]
[1,] -0.4646675  0.833286355  0.2995295
[2,] -0.5537546 -0.009499485 -0.8326258
[3,] -0.6909703 -0.552759994  0.4658502

3.4.9 矩阵的奇异值分解

使用svd()函数对矩阵进行奇异值分解，即 $A=UDV^T$ ，其中 $U, V$ 为正交阵， $D$ 为对角阵，也就是其奇异值。svd()的返回值也是列表。

> svdA = svd(A)
> svdA
$d
[1] 17.4125052  0.8751614  0.1968665

$u
           [,1]        [,2]       [,3]
[1,] -0.2093373  0.96438514  0.1616762
[2,] -0.5038485  0.03532145 -0.8630696
[3,] -0.8380421 -0.26213299  0.4785099

$v
           [,1]         [,2]       [,3]
[1,] -0.4646675 -0.833286355  0.2995295
[2,] -0.5537546  0.009499485 -0.8326258
[3,] -0.6909703  0.552759994  0.4658502
> svdA$u %*% diag(svdA$d) %*% t(svdA$v)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8   10

3.4.10 QR分解

QR分解是将矩阵分解为一个正交矩阵Q和一个上三角矩阵R的方法。使用qr()函数。
$qr是主要结果，其上对角部分是R，下对角部分则包含Q的部分信息（我也不太懂这里的部分信息是怎么个意思，有知道的可以写在评论区）。一般，可以用qr.Q()和qr.R()函数分别得到其QR矩阵结果。
另外，QR分解在计算最小二乘拟合时，计算误差比一般方法要小。

> qrA = qr(A)
> qrA
$qr
           [,1]       [,2]        [,3]
[1,] -8.1240384 -9.6011363 -11.9398746
[2,]  0.4923660  0.9045340   1.5075567
[3,]  0.8616404  0.9954736   0.4082483

$rank
[1] 3

$qraux
[1] 1.1230915 1.0950385 0.4082483

$pivot
[1] 1 2 3

attr(,"class")
[1] "qr"
> qr.Q(qrA)
           [,1]       [,2]       [,3]
[1,] -0.1230915  0.9045340  0.4082483
[2,] -0.4923660  0.3015113 -0.8164966
[3,] -0.8616404 -0.3015113  0.4082483
> qr.R(qrA)
          [,1]      [,2]        [,3]
[1,] -8.124038 -9.601136 -11.9398746
[2,]  0.000000  0.904534   1.5075567
[3,]  0.000000  0.000000   0.4082483

3.5 与矩阵(数组)运算有关的函数

3.5.1 取矩阵的维数

dim(), nrow() 和 ncol()

> A = matrix(1:6, nrow = 2)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> dim(A)
[1] 2 3
> nrow(A)
[1] 2
> ncol(A)
[1] 3

3.5.2 矩阵的合并

cbind()可以按列合并（横向）矩阵，而rbind()可以按行合并（纵向），均要求矩阵对应维度相同，不同的则循环补足，应用于矩阵时矩阵规模固定，但应用于向量时，cbind()中按列向量处理，rbind()中按行向量处理。

> cbind(c(1,2), c(3,4))
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x1 = rbind(c(1,2), c(3,4))
> x1
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> x2 = 10 + x1
> x3 = cbind(x1, x2)
> x3
     [,1] [,2] [,3] [,4]
[1,]    1    2   11   12
[2,]    3    4   13   14
> x4 = rbind(x1, x2)
> x4
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]   11   12
[4,]   13   14
> cbind(1, x1)
     [,1] [,2] [,3]
[1,]    1    1    2
[2,]    1    3    4

3.5.3 矩阵的拉直

使用as.vector()函数将矩阵转化为向量。

> A = matrix(1:6, nrow = 2)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> as.vector(A)
[1] 1 2 3 4 5 6

3.5.4 数组的维名字

数组有一个属性dimnames保存各个维度下标的名字，默认NULL。
下例是三种定义方法，大家实践中应当灵活选用。

> X = matrix(1:6, ncol = 2,
+            dimnames = list(c("one","two","three"),c("First", "Second")), byrow = T)
> X
      First Second
one       1      2
two       3      4
three     5      6
> X1 = matrix(1:6, ncol = 2, byrow = T)
> dimnames(X1) = list(c("one","two","three"),c("First", "Second"))
> X1
      First Second
one       1      2
two       3      4
three     5      6
> X2 = matrix(1:6, ncol = 2, byrow = T)
> colnames(X2) = c("First", "Second")
> rownames(X2) = c("one", "two", "three")
> X2
      First Second
one       1      2
two       3      4
three     5      6

3.5.5 数组的广义转置

使用aperm(A, perm)函数重排各个维度的次序，例如，perm中的2即代表原来的第二个维度。

> A = array(1:24, dim = c(2,3,4))
> A
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

, , 3

     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14   16   18

, , 4

     [,1] [,2] [,3]
[1,]   19   21   23
[2,]   20   22   24

> B = aperm(A, c(2,3,1))
> B
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    3    9   15   21
[3,]    5   11   17   23

, , 2

     [,1] [,2] [,3] [,4]
[1,]    2    8   14   20
[2,]    4   10   16   22
[3,]    6   12   18   24

3.5.6 apply函数

apply()函数用于对数组某一维度或若干维度进行计算，声明如下：
apply(A, MARGIN, FUN, …)
A是数组，MARGIN是固定不变的维度，FUN是处理用的函数，看下例更容易理解：

> A = matrix(1:6, nrow = 2)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> apply(A, 1, sum)
[1]  9 12
> apply(A, 2, mean)
[1] 1.5 3.5 5.5