R语言里的矩阵处理学习笔记
关于矩阵,通常都会使用matlab来做处理。其实使用R也可以对矩阵做出一些简单的处理。而R语言中提供的matrix,matlab包也提供了不少关于矩阵处理的东西(可以通过??matlab来查看具体函数)。
一、矩阵的输入
通常我们使用函数matrix来创建矩阵,函数的介绍如下:
matrix(data = NA,nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
如果想将数据按行输入矩阵,参数byrow须改为T。
由于矩阵也是特殊的数组,我们也可以用生成数组的函数array()。具体格式如下:
array(data = NA, dim= length(data), dimnames = NULL)
这里的dim是一个二维数组,生成的就是矩阵了。
当然dim()函数,attr()(这个是一个格式转换的函数)也是可以生成矩阵的。
还有如diag()可以产生对角矩阵这样特殊矩阵的函数。
例如生成下面的这个矩阵(为了便于下面的叙述,这个矩阵记为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)
二、矩阵的基本操作
这里主要的有:矩阵的加法与乘法,矩阵求秩,矩阵的转置,方阵求行列式,矩阵求逆,解线性方程组
1、矩阵的加法与乘法
加法直接使用加号,实现对应元素相加。但是要注意两个矩阵必须可加
矩阵的乘法:和matlab类似,R也给出了两种乘法:
“*”:对应位置元素做乘法,如x*x得到结果:
> 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)得到结果:
> x%*%t(x)
[,1] [,2]
[1,] 680 744
[2,] 744 816
这里乘法也必须是有意义的才行。
通常我们也使用crossprod()函数来做乘法,crossprod(x,x)效果与x%*%t(x)相同
2、矩阵求秩
这里可以利用qr分解来解决求秩的问题
qr(x)$rank
可以得到x的秩
3、矩阵的转置
常用的命令为t().
R中还可以使用命令aperm()来实现矩阵的广义转置,函数用法如下:
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、矩阵求逆与解线性方程组
使用函数solve()
对于线性方程组b<-A%*%y
的解使用函数solve(A,b)即可
从而我们知道取b为单位阵时即可求逆,通常简化为solve(A)
但是值得注意的是,用直接求逆的办法解线性方程组y<-solve(A)%*%b不仅稳定性低,效率也不咋地。
三、矩阵的分解
这里主要介绍矩阵的lu分解,Cholesky分解,特征值与特征向量,QR分解,奇异值分解
1、LU分解
将系数矩阵A转变成等价两个矩阵L和U的乘积 ,其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时,矩阵A可以分解为A=LU,且分解唯一。其中L是单位下三角矩阵,U是上三角矩阵。
library(Matrix) #这里引入函数包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(这个矩阵的意思是保持被变换的矩阵第一行不变,二三行对调)
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)当限定的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分解求矩阵的逆,这时用函数chol2inv(),这种用法更有效。
函数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可以快速得到结果
4、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 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,可以用函数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 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分解
利用函数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包中的函数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
验证性质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
对于矩阵,我们还可以使用cbind(),rbind()构造分块矩阵。
数据分析咨询请扫描二维码
若不方便扫码,搜微信号:CDAshujufenxi
在当今这个数据驱动的时代,几乎每一个业务决策都离不开对数据的深入分析。而其中,指标波动归因分析更是至关重要的一环。无论是 ...
2025-02-18当数据开始说谎:那些年我们交过的学费 你有没有经历过这样的场景?熬了三个通宵做的数据分析报告,在会议上被老板一句"这数据靠 ...
2025-02-17数据分析作为一门跨学科领域,融合了统计学、编程、业务理解和可视化技术。无论是初学者还是有一定经验的从业者,系统化的学习路 ...
2025-02-17挖掘用户价值本质是让企业从‘赚今天的钱’升级为‘赚未来的钱’,同时让用户从‘被推销’变为‘被满足’。询问deepseek关于挖 ...
2025-02-17近来deepseek爆火,看看deepseek能否帮我们快速实现数据看板实时更新。 可以看出这对不知道怎么动手的小白来说是相当友好的, ...
2025-02-14一秒精通 Deepseek,不用找教程,不用买资料,更不用报一堆垃圾课程,所有这么去做的,都是舍近求远,因为你忽略了 deepseek 的 ...
2025-02-12自学 Python 的关键在于高效规划 + 实践驱动。以下是一份适合零基础快速入门的自学路径,结合资源推荐和实用技巧: 一、快速入 ...
2025-02-12“我们的利润率上升了,但销售额却没变,这是为什么?” “某个业务的市场份额在下滑,到底是什么原因?” “公司整体业绩 ...
2025-02-08活动介绍 为了助力大家在数据分析领域不断精进技能,我们特别举办本期打卡活动。在这里,你可以充分利用碎片化时间在线学习,让 ...
2025-02-071、闺女,醒醒,媒人把相亲的带来了。 我。。。。。。。 2、前年春节相亲相了40个, 去年春节相亲50个, 祖宗,今年你想相多少个 ...
2025-02-06在数据科学的广阔领域中,统计分析与数据挖掘占据了重要位置。尽管它们常常被视为有关联的领域,但两者在理论基础、目标、方法及 ...
2025-02-05在数据分析的世界里,“对比”是一种简单且有效的方法。这就像两个女孩子穿同一款式的衣服,效果不一样。 很多人都听过“货比三 ...
2025-02-05当我们只有非常少量的已标记数据,同时有大量未标记数据点时,可以使用半监督学习算法来处理。在sklearn中,基于图算法的半监督 ...
2025-02-05考虑一种棘手的情况:训练数据中大部分样本没有标签。此时,我们可以考虑使用半监督学习方法来处理。半监督学习能够利用这些额 ...
2025-02-04一、数学函数 1、取整 =INT(数字) 2、求余数 =MOD(除数,被除数) 3、四舍五入 =ROUND(数字,保留小数位数) 4、取绝对值 =AB ...
2025-02-03作者:CDA持证人 余治国 一般各平台出薪资报告,都会哀嚎遍野。举个例子,去年某招聘平台发布《中国女性职场现状调查报告》, ...
2025-02-02真正的数据分析大神是什么样的呢?有人认为他们能轻松驾驭各种分析工具,能够从海量数据中找到潜在关联,或者一眼识别报告中的数 ...
2025-02-01现今社会,“转行”似乎成无数职场人无法回避的话题。但行业就像座围城:外行人看光鲜,内行人看心酸。数据分析这个行业,近几年 ...
2025-01-31本人基本情况: 学校及专业:厦门大学经济学院应用统计 实习经历:快手数据分析、字节数据分析、百度数据分析 Offer情况:北京 ...
2025-01-3001专家简介 徐杨老师,CDA数据科学研究院教研副总监,主要负责CDA认证项目以及机器学习/人工智能类课程的研发与授课,负责过中 ...
2025-01-29