r语言中矩阵QR分解_学习 R 语言:矩阵和数组

本文内容来自《R 语言编程艺术》(The Art of R Programming),有部分修改

矩阵 (Matrix) 是一种特殊的向量,包含两个附加属性:行数和列数。

数组 (Array) 是更一般的对象,可以有多个维度。矩阵是二维数组。

创建矩阵

R 中下标从 1 开始,矩阵按列存储。

使用 matrix() 函数创建矩阵

 matrix(c(1, 2, 3, 4),
nrow=2,
ncol=2,
)print(y)
     [,1] [,2]
[1,] 1 3
[2,] 2 4

在给定所有元素的情况下,nrow 和 ncol 两个参数只需要给出一个

 matrix(1:4,
nrow=2
)print(y)
     [,1] [,2]
[1,] 1 3
[2,] 2 4

另一种方法是创建空的矩阵,再对每一位赋值

 matrix(nrow=2, ncol=2)print(y)
     [,1] [,2]
[1,] NA NA
[2,] NA NA
1] 
     [,1] [,2]
[1,] 1 3
[2,] 2 4

设置参数 byrow=TRUE,可以按行提供的矩阵数据。但矩阵实际还是按列存储

 matrix(1:6,
nrow=2,
byrow=TRUE,
)print(m)
     [,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6

一般矩阵运算

线性代数运算

矩阵乘法

print(y 
     [,1] [,2]
[1,] 7 15
[2,] 10 22

矩阵元素乘法

print(
     [,1] [,2]
[1,] 3 9
[2,] 6 12

矩阵加法

print(y 
     [,1] [,2]
[1,] 2 6
[2,] 4 8

矩阵索引

向量索引同样适用于矩阵索引

提取矩阵列

 matrix(c(1, 2, 3, 4,1, 1, 0, 0,1, 0, 1, 0
),
nrow=4,
)print(z)
     [,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 0
[3,] 3 0 1
[4,] 4 0 0
print(z[, 
     [,1] [,2]
[1,] 1 1
[2,] 1 0
[3,] 0 1
[4,] 0 0

提取矩阵行

 matrix(c(11, 21, 31,12, 22, 32
),
ncol=2,
)print(y)
     [,1] [,2]
[1,] 11 12
[2,] 21 22
[3,] 31 32
print(y[2
     [,1] [,2]
[1,] 21 22
[2,] 31 32
print(y[2
[1] 22 32

对子矩阵赋值

 matrix(1:6,
nrow=3
)print(y)
     [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
[c(
     [,1] [,2]
[1,] 1 8
[2,] 2 5
[3,] 1 12

另一个示例

 matrix(nrow=3, ncol=3)
y matrix(c(4, 5, 2, 3),
nrow=2
)print(y)
     [,1] [,2]
[1,] 4 2
[2,] 5 3
 matrix(1:6,
nrow=3
)print(y)
     [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
print(y[
     [,1] [,2]
[1,] 1 4
[2,] 3 6

扩展案例:图像操作

library(pixmap)
Pixmap image
Type : pixmapGrey
Size : 194x259
Resolution : 1x1
Bounding box : 0 0 259 194
plot(mtrush1)
c530066570ce2a240b8cd47dbb7895fe.png
str(mtrush1)
Formal class 'pixmapGrey' [package "pixmap"] with 6 slots
..@ grey : num [1:194, 1:259] 0.278 0.263 0.239 0.212 0.192 ...
..@ channels: chr "grey"
..@ size : int [1:2] 194 259
..@ cellres : num [1:2] 1 1
..@ bbox : num [1:4] 0 0 259 194
..@ bbcent : logi FALSE
@grey[28, 
0.796078431372549
 mtrush1
mtrush2@grey[84:163, 135:177] 1plot(mtrush2)
bee0513c491965113d5f5f1370d6b539.png
 function(img, rows, cols, q) {
lrows length(rows)
lcols length(cols)
newimg img
randomnoise matrix(nrow=lrows, ncol=lcols, runif(lrows*lcols))
newimg@grey[rows, cols] (1-q) * img@grey[rows, cols] + q * randomnoisereturn(newimg)
}
 blurpart(
mtrush1,84:163,135:177,0.2
)plot(mtrush3)
526d71d17eea68214e0a5551b4ea47a1.png

矩阵元素筛选

 matrix(c(1, 2, 3, 2, 3, 4),
nrow=3
)print(x)
     [,1] [,2]
[1,] 1 2
[2,] 2 3
[3,] 3 4
print(x[x[,
     [,1] [,2]
[1,] 2 3
[2,] 3 4

分析上面的命令

 x[,2] >= 3print(j)
[1] FALSE  TRUE  TRUE
print(x[j,])
     [,1] [,2]
[1,] 2 3
[2,] 3 4

注意:求 j 是向量化计算

另一个示例

 matrix(1:6,
nrow=3
)print(m)
     [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
 c(5, 12, 13)print(m[z %% 2 == 1,])
     [,1] [,2]
[1,] 1 4
[2,] 3 6
print(m[m[,
[1] 3 6

分析

 m[,1] > 1print(a)
[1] FALSE  TRUE  TRUE
 m[,2] > 5print(b)

[1] FALSE FALSE TRUE
 a & bprint(c)
[1] FALSE FALSE  TRUE
print(m[c, ])
[1] 3 6

注意:这里使用 & 表示向量的逻辑与运算。而 && 是用于 if 语句的标量逻辑与运算

向量运算也适用于矩阵

 matrix(c(5, 2, 9, -1, 10, 11),
nrow=3
)print(m)
     [,1] [,2]
[1,] 5 -1
[2,] 2 10
[3,] 9 11
print(
[1] 1 3 5 6

扩展案例:生成协方差矩阵

row() 函数和 col() 函数

协方差矩阵是对称的,对角线元素均为 1

 function(rho, n) {
m matrix(nrow=n, ncol=n)
m ifelse(row(m) == col(m), 1, rho)return(m)
}

示例

print(
     [,1] [,2] [,3]
[1,] 1.0 0.2 0.2
[2,] 0.2 1.0 0.2
[3,] 0.2 0.2 1.0

解释

 matrix(3:8,
nrow=3
)print(z)
     [,1] [,2]
[1,] 3 6
[2,] 4 7
[3,] 5 8
print(
     [,1] [,2]
[1,] 1 1
[2,] 2 2
[3,] 3 3
print(
     [,1] [,2]
[1,] 1 2
[2,] 1 2
[3,] 1 2
print(
      [,1]  [,2]
[1,] TRUE FALSE
[2,] FALSE TRUE
[3,] FALSE FALSE

对矩阵的行和列调用函数

*apply 系列函数:apply, tapply, lapply

apply() 函数在各行或各列上调用指定函数

使用 apply() 函数

一般形式

apply(m, dimcode, f, fargs)

其中

  • m:矩阵

  • dimcode:维度编号,1 表示对每行应用函数,2 表示对每列应用函数

  • f:应用在行或列上的函数,第一个参数为一行或一列

  • fargs:可选参数,用于传递给 f 函数

 matrix(1:6,
nrow=3
)print(z)
     [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6

对每一列应用函数 mean

注:可以直接使用 colmean() 函数

print(
[1] 2 5

对每一行应用自定义函数

 function(x) x/c(2, 8)
y apply(z, 1, f)print(y)
     [,1]  [,2] [,3]
[1,] 0.5 1.000 1.50
[2,] 0.5 0.625 0.75

调用的函数返回 k 个元素的向量,则 apply 函数返回的结果就有 k 行。执行几次函数,则返回结果就有多少列。

如果函数返回一个标量,则 apply 函数返回向量。

t() 函数可以实现转置

print(
     [,1]  [,2]
[1,] 0.5 0.500
[2,] 1.0 0.625
[3,] 1.5 0.750

附加参数

0 和 1 构成的矩阵,对于每行,前 d 个元素中 1 较多则取 1,否则取 0

 function(rw, d) {
maj sum(rw[1:d]) / dreturn(if(maj > 0.5) 1 else 0)
}
 matrix(c(1, 1, 1, 0,0, 1, 0, 1,1, 1, 0, 1,1, 1, 1, 1,0, 0, 1, 0
),
nrow=4
)print(x)
     [,1] [,2] [,3] [,4] [,5]
[1,] 1 0 1 1 0
[2,] 1 1 1 1 0
[3,] 1 0 0 1 1
[4,] 0 1 1 1 0
print(
[1] 1 1 0 1
print(
[1] 0 1 0 0

扩展案例:寻找异常值

识别偏离最远的观测值

 function(x) {
findol function(xrow) {
mdn median(xrow)
devs abs(xrow-mdn)return(which.max(devs))
}return(apply(x, 1, findol))
}

增加或删除矩阵的行或列

矩阵无法改变大小,增加或删除实际上是创建新的对象

改变矩阵大小

重新赋值可以改变向量的大小

 c(12, 5, 13, 16, 8)print(x)
[1] 12  5 13 16  8
 c(x, 20)print(x)
[1] 12  5 13 16  8 20
 c(x[1:3], 20, x[4:6])print(x)
[1] 12  5 13 20 16  8 20
 x[-2:-4]print(x)
[1] 12 16  8 20

使用 rbind() 和 cbind() 为矩阵增加行或列

 c(1, 1, 1, 1)
z matrix(c(1, 2, 3 ,4,1, 1, 0, 0,1, 0, 1, 0
),
nrow=4
)print(z)
     [,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 0
[3,] 3 0 1
[4,] 4 0 0
print(
     one
[1,] 1 1 1 1
[2,] 1 2 1 0
[3,] 1 3 0 1
[4,] 1 4 0 0

可能会用到循环补齐 (recycling)

print(
     [,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 1 0
[3,] 1 3 0 1
[4,] 1 4 0 0

cbind() 和 rbind() 可以用来快速生成小矩阵

 cbind(c(1, 2),c(3, 4)
)print(q)
     [,1] [,2]
[1,] 1 3
[2,] 2 4

注意:创建向量和矩阵均耗费时间,请谨慎使用 cbind 和 rbind

循环添加列或行的方式不可取,应该首先生成一个大矩阵,然后再往里填值。

可以通过重新赋值删除矩阵的行或列

 matrix(1:6, nrow=3)print(m)
     [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
 m[c(1,3),]print(m)
     [,1] [,2]
[1,] 1 4
[2,] 3 6

扩展案例:找到图中距离最近的一对端点

矩阵元素 (i, j) 表示城市 i 和 城市 j 间的距离

imin 函数寻找当前行的最小值,因为矩阵是对称的,所以仅搜索从行号 i + 1 开始到结尾。

返回最小值的索引(即矩阵中的列数)和最小值

 function(x) {
lx length(x)
i x[lx]
j which.min(x[(i+1):(lx-1)])
k i + jreturn(c(k, x[k]))
}

mind 函数为每行分别求最小值 (wmins),再找出全局最小值,返回最小值,行号和列号

 function(d) {
n nrow(d)
dd cbind(d, 1:n)
wmins apply(dd[-n,], 1, imin)
i which.min(wmins[2,])
j wmins[1, i]return(c(d[i, j], i, j))
}

示例

 matrix(c(0, 12, 13, 8, 20,12, 0, 15, 28, 88,13, 15, 0, 6, 9,8, 28, 6, 0, 33,20, 88, 9, 33, 0
),
nrow=5,
byrow=T,
)print(q)
     [,1] [,2] [,3] [,4] [,5]
[1,] 0 12 13 8 20
[2,] 12 0 15 28 88
[3,] 13 15 0 6 9
[4,] 8 28 6 0 33
[5,] 20 88 9 33 0
print(
[1] 6 3 4

分析

增加记录当前行号的一列

 nrow(q)
qq cbind(q, 1:n)print(qq)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 12 13 8 20 1
[2,] 12 0 15 28 88 2
[3,] 13 15 0 6 9 3
[4,] 8 28 6 0 33 4
[5,] 20 88 9 33 0 5

为前 n-1 行,即前 4 行,寻找最小值,并记录最小值(第 2 行)及其对应的列数(第 1 行)。

返回结果中每列代表矩阵中的行,即第 1 列表示矩阵的第 1 行,最后一列表示矩阵的第 4 行

 apply(qq[-n,], 1, imin)print(wmins)
     [,1] [,2] [,3] [,4]
[1,] 4 3 4 5
[2,] 8 15 6 33

求最小值,返回列号,即矩阵中的行号 i

 which.min(wmins[2,])print(i)
[1] 3

从 wmins 中取出矩阵的列号 j

 wmins[1, i]print(j)
[1] 4
print(
[1] 6 3 4

如果最小值唯一,则有更简单的方法

注:有两次矩阵计算,效率不是最优的

 function(d) {
smallest min(d)
ij which(d == smallest, arr.ind=TRUE)return(c(smallest, ij))
}

示例,需要将 0 值和重复值都替换为比较大的数

= q
      [,1]  [,2]  [,3]  [,4]  [,5]
[1,] 99999 12 13 8 20
[2,] 99999 99999 15 28 88
[3,] 99999 99999 99999 6 9
[4,] 99999 99999 99999 99999 33
[5,] 99999 99999 99999 99999 99999
print(
[1] 6 3 4

向量与矩阵的差异

矩阵实际上就是向量,只不过多了两个属性:行数和列数

 matrix(1:8, nrow=4)print(z)
     [,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8

求长度

print(
[1] 8

但 z 不仅仅是向量

print(
[1] "matrix" "array"
print(
$dim
[1] 4 2

用 dim() 函数访问 dim 属性

print(
[1] 4 3

nrow() 和 ncol() 是 dim() 的简单封装

print(
[1] 4
print(

[1] 2
nrow
function (x)
dim(x)[1L]

避免意外降维

从矩阵中提取一行,会减少一维,即返回向量

 matrix(1:8, nrow=4)print(z)
     [,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8
 z[2,]print(r)
[1] 2 6
print(
$dim
[1] 3 2
print(
NULL
str(z)
 int [1:3, 1:2] 1 2 3 4 5 6
str(r)
 int [1:2] 2 6

使用 drop 参数禁止自动减少维度

 z[2,, drop=FALSE]print(r)
     [,1] [,2]
[1,] 2 5
print(
[1] 1 2

可见中括号 [ 实际上也是一个函数

print(z[3, 
[1] 6
print(
[1] 6

as.matrix() 函数可以将向量转为矩阵

 1:3print(u)
[1] 1 2 3
 as.matrix(u)print(v)
     [,1]
[1,] 1
[2,] 2
[3,] 3
print(
NULL
print(
$dim
[1] 3 1

矩阵的行和列的命名问题

矩阵的行和列也可以命名

 matrix(1:4, nrow=2)print(z)
     [,1] [,2]
[1,] 1 3
[2,] 2 4
print(
NULL
colnames(z) 
     a b
[1,] 1 3
[2,] 2 4
print(
[1] "a" "b"
print(z[, 
[1] 1 2

高维数组

三个学生,两次考试,每次考试记录两个成绩

第一次考试的数据,行表示学生,列表示单次考试的两个成绩

 matrix(c(46, 30,21, 25,50, 50
),
nrow=3,
byrow=T
)print(firsttest)
     [,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50

第二次考试的成绩

 matrix(c(46, 43,41, 35,50, 50
),
nrow=3,
byrow=T
)print(secondtest)
     [,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50

合并为一个数组,两次考试形成第三维,叫做层 (layer)

 array(
data=c(firsttest, secondtest),
dim=c(3, 2, 2)
)
print(
$dim
[1] 3 2 2

学生 3 (dim 1) 在第一次考试 (dim 3) 的第二个部分 (dim 2) 中的得分

print(tests[3, 
[1] 50
print(tests)
, , 1

[,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50

, , 2

[,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50

参考

《学习 R 语言:快速开始》

《学习 R 语言:向量》

本文的 Jupyter Notebook 代码请访问如下项目:

https://github.com/perillaroc/the-art-of-r-programming


8e9bf6563456672acf25ab65c64d68dd.png

题图来自 Pixabay。
  • 1
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值