R中的矩阵运算

最新推荐文章于 2023-04-10 16:22:34 发布

狗狗诗人061

最新推荐文章于 2023-04-10 16:22:34 发布

阅读量1.8k

点赞数 1

文章标签： R

本文链接：https://blog.csdn.net/qq_42412550/article/details/100902927

版权

R中的矩阵运算

创建一个向量

> x=c(1,2,3,4)
> x
[1] 1 2 3 4

向量化算子

vec<-function (x){
t(t(as.vector(x)))
}
vech<-function (x){
t(x[lower.tri(x,diag=T)])
}

> x=matrix(1:12,3,4)
> x
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> vec(x)
      [,1]
 [1,]    1
 [2,]    2
 [3,]    3
 [4,]    4
 [5,]    5
 [6,]    6
 [7,]    7
 [8,]    8
 [9,]    9
[10,]   10
[11,]   11
[12,]   12

> vech(x)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    5    6    9

创建一个矩阵

在R中可以用函数**matrix()**来创建一个矩阵，应用该函数时需要输入必要的参数值。

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

data为必要的矩阵元素，是一个向量数据。nrow为行数，ncol为列数，注意nrow与ncol的乘积应为矩阵元素个数。byrow项控制排列元素时是否按行进行，当byrow=T时，生成矩阵的数据按行放置，缺省时相当于byrow=F，数据按列放置，dimnames是数组维的名字，即给定行和列的名称，缺省时为空。

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

> matrix(1:12,nrow=3,ncol=4)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> matrix(1:12,nrow=4,ncol=3)
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

> matrix(1:12,nrow=4,ncol=3,byrow=T)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

> rowname=c("r1","r2","r3")
> rowname
[1] "r1" "r2" "r3"

> colname=c("c1","c2","c3","c4")
> colname
[1] "c1" "c2" "c3" "c4"

> matrix(1:12,nrow=3,ncol=4,dimnames=list(rowname,colname))
   c1 c2 c3 c4
r1  1  4  7 10
r2  2  5  8 11
r3  3  6  9 12

创建一个服从正态分布的随机数矩阵：

> matrix(rnorm(16),nrow=4)
           [,1]        [,2]       [,3]       [,4]
[1,] -1.0404372  0.12465269  0.4497447 -0.7273449
[2,]  0.1113022  0.51174249  0.4619477  1.1711118
[3,] -0.8835299  0.03437571 -1.4223462  1.4187844
[4,]  0.3888612 -1.69345992 -0.3732066  2.2947479

创建数字相同的n列m行矩阵：

matrix(2,ncol=n,nrow=m)

创建对角矩阵：

diag(x,ncol=n,nrow=m)

> diag(c(2,5,1),3,3)
     [,1] [,2] [,3]
[1,]    2    0    0
[2,]    0    5    0
[3,]    0    0    1
> diag(1,4,4)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

矩阵的运算

矩阵的维数

在R中很容易得到一个矩阵的维数，函数**dim()**将返回一个矩阵的维数，**nrow()**返回行数，**ncol()**返回列数。

> A=matrix(1:12,3,4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> nrow(A)
[1] 3
> ncol(A)
[1] 4

向量只有定义了维数向量后才能被看作是数组。

> A=1:12
> dim(A)=c(3,4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
#注意生成矩阵按列排

矩阵转置

A为m×n矩阵，求A’在R中可用函数t()。

> A=matrix(1:12,nrow=3,ncol=4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> t(A)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
[4,]   10   11   12

若将函数t()作用于一个向量x，则R默认x为列向量，返回结果为一个行向量。

> x
[1] 1 2 3 4
> t(x)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
> class(x)
[1] "numeric"
> class(t(x))
[1] "matrix"

若想得到一个列向量，可用t(t(x))

> x
[1] 1 2 3 4
> t(t(x))
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4

> y=t(t(x))
> t(t(y))
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4

矩阵求逆

矩阵求逆可用函数solve()，应用solve(a, b)运算结果是解线性方程组ax = b，若b缺省，则系统默认为单位矩阵，因此可用其进行矩阵求逆。

> A=matrix(rnorm(16),4,4)
> A
           [,1]       [,2]        [,3]       [,4]
[1,] -0.8712401  0.1113022  0.51174249  0.4619477
[2,] -1.0227890 -0.8835299  0.03437571 -1.4223462
[3,] -1.9354118  0.3888612 -1.69345992 -0.3732066
[4,] -1.0404372  0.1246527  0.44974466 -0.7273449

> solve(A) #矩阵的逆
           [,1]        [,2]        [,3]        [,4]
[1,] -0.5895417 -0.15091808 -0.17801790  0.01204078
[2,] -0.7554458 -0.98079749  0.11771324  1.37778621
[3,]  0.3018511 -0.05568050 -0.37015008  0.49052214
[4,]  0.9004925  0.01336331  0.04594339 -0.85265354

> solve(A) %*%A
              [,1]          [,2]          [,3]          [,4]
[1,]  1.000000e+00  7.372575e-18  5.724587e-17 -1.561251e-17
[2,] -2.220446e-16  1.000000e+00  1.110223e-16  2.220446e-16
[3,]  0.000000e+00  1.387779e-17  1.000000e+00  0.000000e+00
[4,]  0.000000e+00 -1.387779e-17 -1.110223e-16  1.000000e+00

矩阵加减乘

在R中对同行同列矩阵相加减，可用符号：“＋”、“－”。

> A=B=matrix(1:12,nrow=3,ncol=4)
> A+B
     [,1] [,2] [,3] [,4]
[1,]    2    8   14   20
[2,]    4   10   16   22
[3,]    6   12   18   24
> A-B
     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
> A*B #对应位置的元素进行运算
     [,1] [,2] [,3] [,4]
[1,]    1   16   49  100
[2,]    4   25   64  121
[3,]    9   36   81  144
> round(A/B,2)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    1    1
[3,]    1    1    1    1

> A <- matrix(1:9,3,3)
> B <- matrix(10:18,3,3,byrow=T)
> round(A/B,2)
     [,1] [,2] [,3]
[1,] 0.10 0.36 0.58
[2,] 0.15 0.36 0.53
[3,] 0.19 0.35 0.50
> round(A/B,3)
      [,1]  [,2]  [,3]
[1,] 0.100 0.364 0.583
[2,] 0.154 0.357 0.533
[3,] 0.188 0.353 0.500

数与矩阵相乘

A为m×n矩阵，c>0，在R中求cA可用符号：“*”。

> A=matrix(rnorm(16),4,4)
> c=2
> c*A
          [,1]       [,2]       [,3]       [,4]
[1,]  2.342224 -0.1160524 -0.1329442 -4.1713005
[2,]  2.837569  1.0332378 -0.4191544 -0.7643557
[3,]  4.589496 -0.9605139  2.1827769  3.1975189
[4,] -2.600122  2.8490885 -0.2917054 -0.1194255

矩阵相乘

A为m×n矩阵，B为n×k矩阵，在R中求AB可用符号：“％*％”，要求第一个矩阵的列数与行数相同。

> A=matrix(1:12,nrow=3,ncol=4)
> B=matrix(1:12,nrow=4,ncol=3)
> A%*%B
     [,1] [,2] [,3]
[1,]   70  158  246
[2,]   80  184  288
[3,]   90  210  330

两个矩阵的内积：若A为n×m矩阵，要得到A’B，可用函数crossprod()，该函数计算结果与t(A)%*%B相同，但是效率更高。

> A=matrix(1:12,nrow=4,ncol=3)
> B=matrix(1:12,nrow=4,ncol=3)

> t(A)%*%B
[,1] [,2] [,3]
[1,] 30 70 110
[2,] 70 174 278
[3,] 110 278 446

> crossprod(A,B)
[,1] [,2] [,3]
[1,] 30 70 110
[2,] 70 174 278
[3,] 110 278 446

两个矩阵的外积：outer()，又叫叉积tcrossprod(A,B)，表示A%*%t(B)。

> tcrossprod(A,B)
     [,1] [,2] [,3] [,4]
[1,]  107  122  137  152
[2,]  122  140  158  176
[3,]  137  158  179  200
[4,]  152  176  200  224

矩阵Hadamard积：若A={aij}m×n, B={bij}m×n, 则矩阵的Hadamard积定义为：A⊙B={aij bij }m×n。
R中Hadamard积可以直接运用运算符“*”。

> A=matrix(1:16,4,4)
> 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
> B=A
> A*B
     [,1] [,2] [,3] [,4]
[1,]    1   25   81  169
[2,]    4   36  100  196
[3,]    9   49  121  225
[4,]   16   64  144  256

矩阵的行和、列和、行平均与列平均

> A=matrix(1:16,4,4)
> 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

> rowSums(A) #行和
[1] 28 32 36 40
> rowMeans(A) #行平均
[1]  7  8  9 10
> colSums(A) #列和
[1] 10 26 42 58
> colMeans(A) #列平均
[1]  2.5  6.5 10.5 14.5

上述关于矩阵行和列的操作，还可以使用**apply()**函数实现。
apply()函数一般有三个参数，第一个参数代表矩阵对象；第二个参数代表要操作矩阵的维度，1表示对行进行处理，2表示对列进行处理；第三个参数就是处理数据的函数。

> A=matrix(1:16,4,4)
> 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
> args(apply)

function (X, MARGIN, FUN, …)
#X为矩阵，MARGIN用来指定是对行运算还是对列运算，MARGIN＝1表示对行运算，MARGIN＝2表示对列运算，FUN用来指定运算函数, …用来给定FUN中需要的其它的参数。

> apply(A,1,sum) #行和
[1] 28 32 36 40
> apply(A,1,mean) #行平均
[1]  7  8  9 10
> apply(A,2,sum) #列和
[1] 10 26 42 58
> apply(A,2,mean) 列平均
[1]  2.5  6.5 10.5 14.5

**apply()**函数功能强大，我们可以对矩阵的行或者列进行其它运算，例如：计算每一列的方差。

> A=matrix(rnorm(100),20,5)
> apply(A,2,var)
[1] 0.5381173 0.9386903 0.8095553 1.0751036 0.4523867

> apply(A,2,function(x,a)x*a,a=2)
             [,1]       [,2]         [,3]       [,4]        [,5]
 [1,]  0.16316464  1.1447812 -0.009093235  2.0000152 -0.33264978
 [2,] -0.01178833 -2.6486865  2.015007204 -1.2134291 -1.58910117
 [3,]  1.26282645 -1.0950278  1.560019080  0.7063075 -1.14669047
 [4,] -1.54228344 -2.6213897 -0.677652808 -0.3721412  2.47238161
 [5,] -0.29886324 -2.2479107  0.451237883  0.3697244 -0.03865527
 [6,] -1.21685700 -0.5023253 -1.280240346  0.7858224 -0.47129852
 [7,] -2.10971059  4.0672129 -4.494801955 -1.4999128 -0.61827955
 [8,]  3.57077083  0.2402861  1.229849894 -1.7856335 -0.83637219
 [9,]  1.62288005  0.7930775 -0.402749107  1.5328548  2.65776496
[10,]  0.57121073 -1.6156684  1.019192208  1.1743835  0.25518786
[11,] -1.83587538 -2.1246697  1.143074435 -5.8006707  0.58895132
[12,] -0.05865423 -1.6899474  0.428914557 -0.3648873 -2.01370231
[13,]  0.71927161 -2.5710846 -0.304246618 -3.2018336 -2.03747487
[14,]  0.92100703 -1.9311210 -1.463643998 -2.8150015 -0.08794078
[15,] -0.49440385  0.9075725  1.852034686 -2.4689345 -1.35004016
[16,]  1.59260273  1.1453355  2.254886784  2.4937487 -1.86244455
[17,] -1.03195424  1.4415105 -1.620298692  0.2978265 -0.51564602
[18,] -0.69378657 -3.9347814 -1.797394644 -0.6810874 -0.01785887
[19,] -2.46964212 -1.7145394  3.287047265 -2.4820883  1.61570536
[20,]  0.24441263 -1.9243323 -1.067046947  1.3601310  0.54606481

注意：apply(A,2,function(x,a)x*a,a=2)与A×2效果相同，此处旨在说明如何应用alpply函数。

矩阵对角元素相关运算

> A=matrix(1:16,nrow=4,ncol=4)
> 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
> diag(A)
[1]  1  6 11 16

对一个向量应用**diag()**函数将产生以这个向量为对角元素的对角矩阵，例如：

> diag(diag(A))
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    6    0    0
[3,]    0    0   11    0
[4,]    0    0    0   16

对一个正整数z应用**diag()**函数将产生以z维单位矩阵，例如：

> diag(3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

取矩阵的上、下三角部分

在R中，我们可以很方便的取到一个矩阵的上、下三角部分的元素，函数**lower.tri()和函数upper.tri()**提供了有效的方法。

> args(lower.tri)
function (x, diag = FALSE)

函数将返回一个逻辑值矩阵，其中下三角部分为真，上三角部分为假，选项diag为真时包含对角元素，为假时不包含对角元素。upper.tri()的效果与之相反。

> A=matrix(1:12,3,4)
> A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> lower.tri(A)
      [,1]  [,2]  [,3]  [,4]
[1,] FALSE FALSE FALSE FALSE
[2,]  TRUE FALSE FALSE FALSE
[3,]  TRUE  TRUE FALSE FALSE

> lower.tri(A,diag=T)
     [,1]  [,2]  [,3]  [,4]
[1,] TRUE FALSE FALSE FALSE
[2,] TRUE  TRUE FALSE FALSE
[3,] TRUE  TRUE  TRUE FALSE

> upper.tri(A)
      [,1]  [,2]  [,3] [,4]
[1,] FALSE  TRUE  TRUE TRUE
[2,] FALSE FALSE  TRUE TRUE
[3,] FALSE FALSE FALSE TRUE

> upper.tri(A,diag=T)
      [,1]  [,2] [,3] [,4]
[1,]  TRUE  TRUE TRUE TRUE
[2,] FALSE  TRUE TRUE TRUE
[3,] FALSE FALSE TRUE TRUE

> A[lower.tri(A)]=0
> A
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    0    5    8   11
[3,]    0    0    9   12
> A[upper.tri(A)]=0
> A
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    5    0    0
[3,]    0    0    9    0

矩阵的特征值与特征向量

矩阵A的谱分解为A=UΛU’,其中Λ是由A的特征值组成的对角矩阵，U的列为A的特征值对应的特征向量，在R中可以用函数**eigen()**得到U和Λ，

> args(eigen)
function (x, symmetric, only.values = FALSE, EISPACK = FALSE)

其中：x为矩阵，symmetric项指定矩阵x是否为对称矩阵，若不指定，系统将自动检测x是否为对称矩阵。

> B=diag(4)+1
> B
     [,1] [,2] [,3] [,4]
[1,]    2    1    1    1
[2,]    1    2    1    1
[3,]    1    1    2    1
[4,]    1    1    1    2
> B.eigen=eigen(B,symmetric=T)
> B.eigen
eigen() decomposition
$values
[1] 5 1 1 1
$vectors
     [,1]       [,2]       [,3]       [,4]
[1,] -0.5  0.8660254  0.0000000  0.0000000
[2,] -0.5 -0.2886751 -0.5773503 -0.5773503
[3,] -0.5 -0.2886751 -0.2113249  0.7886751
[4,] -0.5 -0.2886751  0.7886751 -0.2113249
> B.eigen$vectors%*%diag(B.eigen$vectors)
              [,1]
[1,] -5.551115e-17
[2,]  5.773503e-01
[3,]  2.113249e-01
[4,]  2.113249e-01
> t(B.eigen$vectors)%*%B.eigen$vectors
              [,1]          [,2]          [,3]          [,4]
[1,]  1.000000e+00 -5.551115e-17 -1.110223e-16 -9.714451e-17
[2,] -5.551115e-17  1.000000e+00 -5.551115e-17 -5.551115e-17
[3,] -1.110223e-16 -5.551115e-17  1.000000e+00  0.000000e+00
[4,] -9.714451e-17 -5.551115e-17  0.000000e+00  1.000000e+00

矩阵广义逆(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包中的函数**ginv()**可计算矩阵A的Moore-Penrose逆。

> library("MASS")
> A=matrix(1:16,4,4)
> 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

验证性质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

验证性质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

验证性质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

验证性质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

矩阵X’X的逆

在统计计算中，我们常常需要计算这样矩阵的逆，如OLS估计中求系数矩阵。R中的包“strucchange”提供了有效的计算方法。

> library("strucchange")
> args(solveCrossprod)
function (X, method = c("qr", "chol", "solve"))

method指定求逆方法，选用“qr”效率最高，选用“chol”精度最高，选用“slove”与slove(crossprod(x,x))效果相同

> A=matrix(rnorm(16),4,4)
> solveCrossprod(A,method="qr")
         [,1]     [,2]      [,3]     [,4]
[1,] 1.855991 1.132912  5.190787 1.874650
[2,] 1.132912 1.067508  3.881481 1.283942
[3,] 5.190787 3.881481 17.581683 6.265738
[4,] 1.874650 1.283942  6.265738 2.450694

> solveCrossprod(A,method="chol")
         [,1]     [,2]      [,3]     [,4]
[1,] 1.855991 1.132912  5.190787 1.874650
[2,] 1.132912 1.067508  3.881481 1.283942
[3,] 5.190787 3.881481 17.581683 6.265738
[4,] 1.874650 1.283942  6.265738 2.450694

> solveCrossprod(A,method="solve")
         [,1]     [,2]      [,3]     [,4]
[1,] 1.855991 1.132912  5.190787 1.874650
[2,] 1.132912 1.067508  3.881481 1.283942
[3,] 5.190787 3.881481 17.581683 6.265738
[4,] 1.874650 1.283942  6.265738 2.450694

> solve(crossprod(A,A))
         [,1]     [,2]      [,3]     [,4]
[1,] 1.855991 1.132912  5.190787 1.874650
[2,] 1.132912 1.067508  3.881481 1.283942
[3,] 5.190787 3.881481 17.581683 6.265738
[4,] 1.874650 1.283942  6.265738 2.450694

矩阵Kronecker积

n×m矩阵A与h×k矩阵B的kronecker积为一个nh×mk维矩阵，在R中kronecker积可以用函数**kronecker()**来计算。

> A=matrix(1:4,2,2)
> B=matrix(rep(1,4),2,2)
> A
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> B
     [,1] [,2]
[1,]    1    1
[2,]    1    1

> kronecker(A,B)
     [,1] [,2] [,3] [,4]
[1,]    1    1    3    3
[2,]    1    1    3    3
[3,]    2    2    4    4
[4,]    2    2    4    4

矩阵的Choleskey分解

对于正定矩阵A，可对其进行Choleskey分解，即：A=P’P，其中P为上三角矩阵，在R中可以用函数**chol()**进行Choleskey分解。

> 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)  #A=P'P，求上三角阵P
[,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

若矩阵为对称正定矩阵，可以利用Choleskey分解求行列式的值。

> prod(diag(chol(A))^2)
[1] 5
> det(A)
[1] 5

若矩阵为对称正定矩阵，可以利用Choleskey分解求矩阵的逆，这时用函数chol2inv()，这种用法更有效。

> chol2inv(chol(A))
[,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8

> solve(A)
[,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8

矩阵奇异值分解

A为m×n矩阵，rank(A)= r, 可以分解为：A=UDV’,其中U’U=V’V=I。在R中可以用函数**scd()**进行奇异值分解。

> A=matrix(1:18,3,6)
> A
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    4    7   10   13   16
[2,]    2    5    8   11   14   17
[3,]    3    6    9   12   15   18

> svd(A)
$d
[1] 4.589453e+01 1.640705e+00 1.366522e-15
$u
           [,1]        [,2]       [,3]
[1,] -0.5290354  0.74394551  0.4082483
[2,] -0.5760715  0.03840487 -0.8164966
[3,] -0.6231077 -0.66713577  0.4082483
$v
            [,1]        [,2]       [,3]
[1,] -0.07736219 -0.71960032 -0.4076688
[2,] -0.19033085 -0.50893247  0.5745647
[3,] -0.30329950 -0.29826463 -0.0280114
[4,] -0.41626816 -0.08759679  0.2226621
[5,] -0.52923682  0.12307105 -0.6212052
[6,] -0.64220548  0.33373889  0.2596585

> A.svd=svd(A)
> A.svd$u%*%diag(A.svd$d)%*%t(A.svd$v)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    4    7   10   13   16
[2,]    2    5    8   11   14   17
[3,]    3    6    9   12   15   18

> t(A.svd$u)%*%A.svd$u
             [,1]         [,2]         [,3]
[1,] 1.000000e+00 3.330669e-16 1.665335e-16
[2,] 3.330669e-16 1.000000e+00 5.551115e-17
[3,] 1.665335e-16 5.551115e-17 1.000000e+00

> t(A.svd$v)%*%A.svd$v
             [,1]          [,2]          [,3]
[1,] 1.000000e+00  2.775558e-17  2.775558e-17
[2,] 2.775558e-17  1.000000e+00 -2.081668e-16
[3,] 2.775558e-17 -2.081668e-16  1.000000e+00

矩阵QR分解

A为m×n矩阵可以进行QR分解，A=QR，其中：Q’Q＝I，在R中可以用函数**qr()**进行QR分解。

> A=matrix(1:16,4,4)
> 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  1.601186e-15  2.217027e-15
[4,]  0.7302967  -0.9124744 -5.547002e-01 -1.478018e-15

$rank
[1] 2
$qraux
[1] 1.182574e+00 1.156135e+00 1.832050e+00 1.478018e-15
$pivot
[1] 1 2 3 4

attr(,"class")
[1] "qr"

rank项返回矩阵的秩，qr项包含了矩阵Q和R的信息，要得到矩阵Q和R，可以用函数qr.Q()和qr.R()作用qr()的返回结果。

> 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  1.601186e-15  2.217027e-15
[4,]  0.000000   0.000000  0.000000e+00 -1.478018e-15

> 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 -1.665335e-16  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 -5.551115e-17  0.000000e+00  2.081668e-17
[2,] -5.551115e-17  1.000000e+00 -2.775558e-17 -6.938894e-17
[3,]  0.000000e+00 -2.775558e-17  1.000000e+00  2.775558e-17
[4,]  2.081668e-17 -6.938894e-17  2.775558e-17  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

backsolve&fowardsolve函数

这两个函数用于解特殊线性方程组，其特殊之处在于系数矩阵为上或下三角。

> args(backsolve)
function (r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)
> args(forwardsolve)
function (l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)

r或者l为n×n维三角矩阵，x为n×1维向量，对给定不同的upper.tri和transpose的值，方程的形式不同。
对于函数backsolve()而言：

> A=matrix(1:9,3,3)
> A
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> x=c(1,2,3)
> x
[1] 1 2 3

> B=A
> B[upper.tri(B)]=0
> B
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    2    5    0
[3,]    3    6    9

> C=A
> C[lower.tri(C)]=0
> C
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    0    5    8
[3,]    0    0    9

> backsolve(A,x,upper.tri=T,transpose=T)
[1]  1.00000000 -0.40000000 -0.08888889

> solve(t(C),x)
[1]  1.00000000 -0.40000000 -0.08888889

> backsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333  0.3333333

> solve(C,x)
[1] -0.8000000 -0.1333333  0.3333333

> backsolve(A,x,upper.tri=F,transpose=T)
[1] 0.0000000 0.0000000 0.3333333

> solve(t(B),x)
[1] 0.0000000 0.0000000 0.3333333

> backsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0

> solve(B,x)
[1] 1 0 0

对于函数forwardsolve()而言：

> A=matrix(1:9,3,3)
> A
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> x=c(1,2,3)
> x
[1] 1 2 3

> B=A
> B[upper.tri(B)]=0
> B
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    2    5    0
[3,]    3    6    9

> C=A
> C[lower.tri(C)]=0
> C
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    0    5    8
[3,]    0    0    9

> forwardsolve(A,x,upper.tri=T,transpose=T)
[1]  1.00000000 -0.40000000 -0.08888889

> solve(t(C),x)
[1]  1.00000000 -0.40000000 -0.08888889

> forwardsolve(A,x,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333  0.3333333

> solve(C,x)
[1] -0.8000000 -0.1333333  0.3333333

> forwardsolve(A,x,upper.tri=F,transpose=T)
[1] 0.0000000 0.0000000 0.3333333

> solve(t(B),x)
[1] 0.0000000 0.0000000 0.3333333

> forwardsolve(A,x,upper.tri=F,transpose=F)
[1] 1 0 0

> solve(B,x)
[1] 1 0 0

row()与col()函数

在R中定义了的这两个函数用于取矩阵元素的行或列下标矩阵，例如矩阵A={aij}m×n，row()函数将返回一个与矩阵A有相同维数的矩阵，该矩阵的第i行第j列元素为i，函数col()类似。

> x=matrix(1:12,3,4)
> row(x)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    2    2    2    2
[3,]    3    3    3    3
> col(x)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
[3,]    1    2    3    4

这两个函数同样可以用于取一个矩阵的上下三角矩阵：

> x
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> x[row(x)<col(x)]=0
> x
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    2    5    0    0
[3,]    3    6    9    0

> x=matrix(1:12,3,4)

> x[row(x)>col(x)]=0
> x
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    0    5    8   11
[3,]    0    0    9   12

行列式的值

在R中，函数det(x)将计算方阵x的行列式的值。
注意：x必必须是正方形矩阵。

> x=matrix(rnorm(16),4,4)
> x
           [,1]       [,2]      [,3]       [,4]
[1,]  0.1848622  0.7664274 -1.600917  0.1489132
[2,]  0.3929112  0.5871918 -1.407501 -0.3405437
[3,] -0.7499564 -2.9003354 -1.234467 -1.2410442
[4,] -0.8928168 -0.1824437  1.246874  0.6800655

> det(x)
[1] -1.195379

狗狗诗人061

关注

1
点赞
踩
15

收藏

觉得还不错? 一键收藏
0
评论
R中的矩阵运算

R中的矩阵运算创建一个向量> x=c(1,2,3,4)> x[1] 1 2 3 4创建一个矩阵在R中可以用函数**matrix()**来创建一个矩阵，应用该函数时需要输入必要的参数值。function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)data项为必要的矩阵元素，nro...
复制链接

扫一扫