R程序调用中的下标出界

lrl

写了两个R代码，一个调用另一个，程序运行到调用代码时，会报“下标出界”的错误。

但是把第一个代码的计算结果保存下来，直接作为被调用代码的输入的时候，运行是正常的，

请问这是什么错误？[s:14]

nan.xiao

过于抽象的问题。。。贴代码吧。。。

lrl

回复第2楼的 nan.xiao：

代码有点长：

主程序：

</p>
<p>### test_lfda.r<br />
source("D:/lfda.r")</p>
<p>### 生成样本数据<br />
n1a=100<br />
n1b=100<br />
n2=100<br />
A=matrix(c(1,2),nrow=2)<br />
B=matrix(c(-6,0),nrow=2)<br />
C=matrix(c(6,0),nrow=2)<br />
D=matrix(c(0,0),nrow=2)</p>
<p>rand_1a=matrix(rnorm(2*n1a),2,n1a)<br />
rep_1a=matrix(rep(A,n1a),2,n1a,byrow=FALSE)<br />
rep_2a=matrix(rep(B,n1a),2,n1a,byrow=FALSE)<br />
X1a=rand_1a*rep_1a+rep_2a</p>
<p>rand_1b=matrix(rnorm(2*n1b),2,n1b)<br />
rep_1b=matrix(rep(A,n1b),2,n1b,byrow=FALSE)<br />
rep_2b=matrix(rep(C,n1b),2,n1b,byrow=FALSE)<br />
X1b=rand_1b*rep_1b+rep_2b</p>
<p>rand_2=matrix(rnorm(2*n2),2,n2)<br />
rep_2m=matrix(rep(A,n2),2,n2,byrow=FALSE)<br />
rep_2n=matrix(rep(D,n2),2,n2,byrow=FALSE)<br />
X2=rand_2*rep_2m+rep_2n</p>
<p>X=matrix(c(X1a,X1b,X2),2,n1a+n1b+n2,byrow=FALSE)<br />
Y=matrix(c(rep(1,n1a+n1b),rep(2,n2)))</p>
<p>### 调用lfda<br />
coords=lfda(X,Y,1,7)<br />

被调函数：

<br />
## lfda.r</p>
<p>lfda<-function(x,y,r,kNN){<br />
## x: data matrix, dimension*sample<br />
## y: class label matrix, sample*1</p>
<p>n=ncol(x)		## number of sample<br />
d=nrow(x)		## number of dimension</p>
<p>tSb=matrix(0,d,d)	##dimension * dimension<br />
tSw=matrix(0,d,d)	##dimension * dimension</p>
<p>for (c in unique(y)){<br />
xc=x[,y==c]			# samples in the same class, dimension*sample<br />
nc=dim(xc)[2]		# number of samples in the same class</p>
<p>##------------ Define classwise affinity matrix A ------------</p>
<p>A=matrix(0,nc,nc)</p>
<p>xc.distances=dist(xc)<br />
xc.distances=as.matrix(xc.distances)<br />
xc.sorted=apply(xc.distances,1,sort)<br />
kNNdist=xc.sorted[kNN+1,]<br />
sigma=kNNdist<br />
localscale=sigma%*%t(sigma)<br />
if(localscale!=0)<br />
{<br />
xc.distances2=xc.distances*xc.distances<br />
A=exp(-xc.distances2/localscale)		## A=exp((-||xi-xj||^2)/sigma^2)<br />
}</p>
<p>##------------ Within-class & between-class scatter matrix ---------</p>
<p>xc1=as.matrix(apply(xc,1,sum))</p>
<p>temp_A=as.matrix(apply(A,1,sum))<br />
An=kronecker(matrix(1,1,d),temp_A)<br />
G=xc%*%(An*t(xc))-xc%*%A%*%t(xc)<br />
tSb=tSb+G/n+xc%*%t(xc)*(1-nc/n)+xc1%*%t(xc1)/n<br />
tSw=tSw+G/nc<br />
}</p>
<p>x1=as.matrix(apply(x,1,sum))<br />
tSb=tSb-x1%*%t(x1)/n-tSw</p>
<p>tSb=(tSb+t(tSb))/2<br />
tSw=(tSw+t(tSw))/2</p>
<p>##------------ Egienevector ---------------------<br />
## AV=λBV  <=>  (B^(-1)*A)V=λV </p>
<p>tSw=solve(tSw)<br />
tS=tSw%*%tSb<br />
a=eigen(tS)$values<br />
b=eigen(tS)$vectors</p>
<p>b=b[,sort.list(a,decreasing=TRUE)]</p>
<p>coords=as.matrix(b[,1:r])</p>
<p>return(coords)</p>
<p>}<br />

lrl

回复第2楼的 nan.xiao：

直接运行会出错，然后我就把第一个程序生成的变量作为第二个程序的输入,在R console中用命令的方式传过去：

<br />
 x=X<br />
 y=Y<br />
 r=1<br />
 kNN=7<br />

再直接把lfda的代码站到R console里面一行行的运行，就没有报错了

lrl

<br />
traceback()<br />

了一下，给出了下列报错信息：

4: lfda(X, Y, 1, 7)

3: eval.with.vis(expr, envir, enclos)

2: eval.with.vis(ei, envir)

1: source("test_lfda.r")

nan.xiao

看了下，应该是 lfda() 那个大循环里的具体实现细节有问题。这个你是从 Masashi Sugiyama 的代码移植的吧，最好一步一步调试过来，保证两边结果完全一致。

我也移植过一个，但是 R 中缺个迭代法求解广义特征值的基础工具，所以后半部分没法写下去。不过可以参考一下我这个循环的细节。

https://github.com/road2stat/sdml/blob/master/R/lfda.R

<br />
x   = t(X)<br />
y   = as.vector(Y)<br />
r   = 1<br />
knn = 7</p>
<p># essential helper function</p>
<p>repmat <- function(A, N, M) {<br />
  kronecker(matrix(1, N, M), A)<br />
}</p>
<p># the for loop</p>
<p>x = t(as.matrix(x))<br />
y = t(as.matrix(y))<br />
d = nrow(x)<br />
n = ncol(x)<br />
if(is.null(r)) r = d</p>
<p>tSb = mat.or.vec(d, d)<br />
tSw = mat.or.vec(d, d)</p>
<p>for (i in unique(as.vector(t(y)))) {</p>
<p>  Xc = x[, y == i]<br />
  nc = ncol(Xc)</p>
<p>  # Define classwise affinity matrix<br />
  Xc2 = t(as.matrix(colSums(Xc^2)))<br />
  distance2 = repmat(Xc2, nc, 1) + repmat(t(Xc2), 1, nc) - 2 * t(Xc) %*% Xc</p>
<p>  sorted = apply(distance2, 2, sort)<br />
  index = apply(distance2, 2, order)</p>
<p>  kNNdist2 = t(as.matrix(sorted[knn + 1, ]))<br />
  sigma = sqrt(kNNdist2)<br />
  localscale = t(sigma) %*% sigma<br />
  flag = (localscale != 0)<br />
  A = mat.or.vec(nc, nc)<br />
  A[flag] = exp(-distance2[flag]/localscale[flag])</p>
<p>  Xc1 = as.matrix(rowSums(Xc))<br />
  G = Xc %*% (repmat(as.matrix(colSums(A)), 1, d) * t(Xc)) - Xc %*% A %*% t(Xc)<br />
  tSb = tSb + (G/n) + Xc %*% t(Xc) * (1 - nc/n) + Xc1 %*% (t(Xc1)/n)<br />
  tSw = tSw + G/nc<br />
}</p>
<p>tSb<br />
# [,1]       [,2]<br />
# [1,] 2459.427040  -1.445689<br />
# [2,]   -1.445689 486.245238<br />
tSw<br />
# [,1]        [,2]<br />
# [1,] 2.125304811 0.008843023<br />
# [2,] 0.008843023 2.407482248<br />

不过你这个 AV=λBV <=> (B^(-1)*A)V=λV 明显不靠谱，数据稍一稀疏方程组很可能就近奇异解不出来了，数值效率方面也有一些问题。所以现阶段估计还是得依赖 MATLAB 的 eigs。

lrl

回复第6楼的 nan.xiao：

看了版主的代码。学习了~~

我也觉得AV=λBV <=> (B^(-1)*A)V=λV太弱了，不过目前还不知道R里有什么解广义特征值的函数，R也不能调用Matlab的eigs函数[s:14]

想问下

r==d时

eigen(solve(tSw) %*% tSb)

是什么原理？

nan.xiao

R 中支持调用 MATLAB 函数做计算然后返回结果。参考 R.matlab 包。

这里的 solve, 只能处理 tSw 非奇异的情况。事实上应该加个判断，如果是奇异的，则用 QZ 法解广义特征值。

而 MATLAB 中的 eig 做得比较友好，自动实现了这个过程：

If B is nonsingular, the problem could be solved by reducing it to a standard eigenvalue problem B^(–1) A x = λ x.

Because B can be singular, an alternative algorithm, called the QZ method, is necessary.

所以尽早实现 ARPACK 的 R 接口是非常必要的。

lrl

回复第8楼的 nan.xiao：

借鉴版主的代码，我用R调用Matlab函数实现了一个LFDA

把下面的连接Matlab Server的代码一行一行在R console里面跑，很正常，但是放在代码里就有问题。

<br />
library(R.matlab)<br />
Matlab$startServer()<br />
matlab <- Matlab(host="localhost",port=9999)<br />
Isopen<-open(matlab)<br />

我用测试数据在Matlab的LFDA和自己实现的LFDA中同时跑，一行一行比对，最后得到的结果是一样的。但是直接调用函数跑结果就不太一样。

这是我的程序：

https://gist.github.com/3076933

继续完善中[s:13]

nan.xiao

用来测试的例子似乎过于正统了，不足以覆盖各种变态的数据。。。在这个数据下纯 R 和纯 MATLAB 得到一样的结果是很正常的。。。混编还是最终出路。

仔细检查一下返回的结果是不是细节上出了问题，比如结果默认排了序或者默认没排序这类。