问个白痴问题——矩阵的N次方怎么算？

Cloudly

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

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

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

aristotle

最笨的方法，呵呵：

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

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

ypchen

根据 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>

Cloudly

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

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

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

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

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

ypchen

回复第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 />

但是我觉得如果 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>

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

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

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

Cloudly

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

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

至于递归和循环，指数和线性的关系么？