• R语言
  • 问个白痴问题——矩阵的N次方怎么算?

从网上搜了半天,好像只有 %*%可以做矩阵乘法。

如果我想算一个矩阵的4,5,6次方,我应该怎么操作?试了试直接^4好像不行(记得matlab里面好像是这样的),只是把矩阵里的每个数分别算了个四次方。

问题白痴,还请别见怪。或许R不应该用来计算这种东西,但是我觉得应该能算吧?懒的再装matlab了。

最笨的方法,呵呵:

<br />
X%*%X  ##二次方<br />
X%*%X*%X ##三次方<br />
X%*%X%*%X%*%X  ##四次方<br />
</p>

比较赞、值得收藏的方法:http://cos.name/cn/topic/5619

根据 http://bm2.genes.nig.ac.jp/RGM2/R_current/library/Biodem/man/mtx.exp.html

<br />
> library(Biodem)<br />
> test<-matrix(c(1:16), 4,4)<br />
> pow.test<-mtx.exp(test,10)<br />
> pow.test<br />
             [,1]         [,2]         [,3]         [,4]<br />
[1,] 2.568153e+14 5.933637e+14 9.299121e+14 1.266461e+15<br />
[2,] 2.908200e+14 6.719305e+14 1.053041e+15 1.434152e+15<br />
[3,] 3.248247e+14 7.504974e+14 1.176170e+15 1.601843e+15<br />
[4,] 3.588295e+14 8.290642e+14 1.299299e+15 1.769534e+15</p>
<p>
</p>

……我想了半天,搜了立方、N次方、power等等,就没想到“乘方”,汗。

yp版主的解决办法真赞,这个函数挺好用的,嘿嘿。对我这种懒人来说,拷贝一个程序不如下个包来得快。

顺便说一下在看到这个帖子前我的解决办法……别见笑。

<br />
b <- a %*% a<br />
b <- b%*%b<br />


然后我就不停的按上箭头和回车重复……汗。其实主要是想看一个转移矩阵的长期的稳定值……后来发现干脆用手算解方程,再汗。
</p>

回复 第4楼 的 Cloudly:

其实可以用递归的方法构造矩阵乘方的函数

<br />
mpower  <- function(mat,n){<br />
	if(n == 1)<br />
		return(mat)<br />
	else<br />
		return(mat %*% mpower(mat,n-1))<br />
}<br />
</p>

但是我觉得如果 n很大的话 一来速度会变慢 二来误差也会积累下去

刚才下载了Biodem的包 看到mtx.exp的代码 原来也这么简单 循环比递归快

<br />
mtx.exp <- function (X, n)<br />
{<br />
    if (n != round(n)) {<br />
        n <- round(n)<br />
        warning("rounding exponent `n' to", n)<br />
    }<br />
    phi <- diag(nrow = nrow(X))<br />
    pot <- X<br />
    while (n > 0) {<br />
        if (n%%2)<br />
            phi <- phi %*% pot<br />
        n <- n%/%2<br />
        pot <- pot %*% pot<br />
    }<br />
    return(phi)<br />
}</p>
<p>
</p>

PS: 转移矩阵是是对称矩阵 它的乘方是有规律可寻的

该矩阵可以分解成 A=LDL' D是个对角阵

A^n = L(D^n)L' 而对角阵的n次方也还是对角阵 对应元素做乘方就行

回复 第5楼 的 ypchen:……偏偏俺遇到的不是个对称阵……可能没说明,是马尔可夫转移矩阵,各种状态之间转移的……所以不一定会是对角阵吧?(状态A->B的概率一定等于B->A的概率?)

其实马尔可夫转移矩阵的话,手动解方程就能求长期稳定状态了,汗。我开始只是不想解方程,就想算算N次方,后来还是乖乖的解方程去了(其实貌似也可以用R解)……

至于递归和循环,指数和线性的关系么?