Nothing
setClass(
Class="CURobj",
representation=representation(
C='matrix',
U='matrix',
R='matrix',
C.leverage.score='numeric',
R.leverage.score='numeric',
C.index='integer',
R.index='integer',
Error='numeric'
)
)
CUR <- function(A, c=dim(A)[2], r=dim(A)[1], k=NULL, sv=NULL,
method="random", alpha=1, weighted=FALSE, beta=4,
matrix.return=TRUE, error.return=FALSE)
{
#function to compute the leverage scores
levscores <- function(v,k) {
if (k==1) {
v[,1]^2
} else {
if (weighted) {qq<-sv$d[1:k]^beta; scores=apply(apply(v[,1:k]^2,1,function(x) x*qq),2,sum); scores/sum(scores)} else apply(v[,1:k]^2,1,sum)/k
}
}
method=match.arg(method,c("random","exact.num.random","top.scores","ortho.top.scores","highest.ranks"))
m=dim(A)[1];n=dim(A)[2]
ScoresC=NULL
ScoresR=NULL
Error=NULL
if(is.null(k)) {
if(is.null(sv)) sv=svd(A)
cs=cumsum(sv$d)
k=length(cs[cs<cs[length(cs)]*0.8])
}
if(is.null(sv)) sv=svd(A,k,k)
if(is.null(rownames(A))){
rownames(A)=1:dim(A)[1]
}
if(is.null(colnames(A))){
colnames(A)=1:dim(A)[2]
}
if (method=="random") {
#select columns randomly
if (c==n) IndC = 1:c else {ScoresC=levscores(sv$v,k); repeat { IndC = which(c*ScoresC>=runif(n)); if (length(IndC)>0) break }} #we repeat until at least 1 column is selected (it is a concern if c is small)
if (r==m) IndR = 1:r else {ScoresR=levscores(sv$u,k); repeat { IndR = which(r*ScoresR>=runif(m)); if (length(IndR)>0) break }}
} else if (method=="exact.num.random") {
#select randomly, but the threshold is adjusted, to pick exatly c columns
if (c==n) IndC = 1:c else {ScoresC=levscores(sv$v,k); IndC = order(c*ScoresC-runif(n),decreasing=TRUE)[1:c]}
if (r==m) IndR = 1:r else {ScoresR=levscores(sv$u,k); IndR = order(r*ScoresR-runif(m),decreasing=TRUE)[1:r]}
} else if (method=="top.scores") {
#select the columns with highest leverage score
if (c==n) IndC = 1:c else {ScoresC=levscores(sv$v,k); IndC = order(ScoresC,decreasing=TRUE)[1:c]}
if (r==m) IndR = 1:r else {ScoresR=levscores(sv$u,k); IndR = order(ScoresR,decreasing=TRUE)[1:r]}
} else if (method=="ortho.top.scores") {
if (c==n) {
IndC = 1:c
} else {
pi=levscores(sv$v,k)
ScoresC=pi
IndC=integer();
vnormsqr=apply(A,2,function(x) sum(x^2)) #Norm square of the columns of A: <c_i|c_i>.
vort=A #Orthogonal part of the vectors. Initially the whole vector is orthogonal.
vortnormsqr=vnormsqr #Norm square of the orthogonal part.
for (i in 1:c)
{
IndC[i]=which.max(pi+alpha*sqrt(vortnormsqr/vnormsqr)) #The next column is selected
sn = vortnormsqr[IndC[i]] # sn = <v_S|v_S>
if (sn>0) {
delta = (vort[,IndC[i]] %*% vort); #Scalar product of the selected vector with all of the vectors: <v_i|v_S>
vort = vort - vort[,IndC[i]] %*% delta /sn #The part, parallel with the new vector is substracted: v_i -> v_i-<v_i|v_S>*v_S/<v_S|v_S>
vortnormsqr = pmax(vortnormsqr-delta^2/sn,0) #The norm square is adjusted accordingly: <v_i|v_i> -> <v_i|v_i>-<v_i|v_S>^2/<v_S|v_S>. Rounding errors may set vortnormsqr negative, pmax mitigates this.
} #(sn>0)
pi[IndC[i]] = 0 #To avoid being selected again.
}
} # c<n
#and the rows
if (r==m) {
IndR = 1:r
} else {
pi=levscores(sv$u,k)
ScoresR=pi
IndR=integer();
vnormsqr=apply(A,1,function(x) sum(x^2))
vort=A #orthogonal part of the vectors
vortnormsqr=vnormsqr
for (i in 1:r)
{
IndR[i]=which.max(pi+alpha*sqrt(vortnormsqr/vnormsqr))
sn = vortnormsqr[IndR[i]]
if (sn>0) {
delta = (vort %*% vort[IndR[i],]);
vort = vort - delta %*% vort[IndR[i],] /sn
vortnormsqr = pmax(vortnormsqr-delta^2/sn,0)
} #(sn>0)
pi[IndR[i]] = 0 #To avoid being selected again.
}
} #r<m
} else if (method=="highest.ranks") {
if (k==1) stop("For k=1 use method=top.scores instead.")
if (c==n) {
IndC = 1:c
} else {
crk=apply(apply(sv$v[,1:k]^2,1,cumsum),1,rank) #cumulative rank of the leverage scores
w=apply(crk,1,which.max) #position of maxima
mx=apply(crk,1,max) #value of maxima
ScoresC=w/(n+1)-mx
IndC=order(-mx,w)[1:c]
crk=NULL #to release memory
}
#and the rows
if (r==m) {
IndR = 1:r
} else {
crk=apply(apply(sv$u[,1:k]^2,1,cumsum),1,rank) #cumulative rank of the leverage scores
w=apply(crk,1,which.max) #position of maxima
mx=apply(crk,1,max) #value of maxima
ScoresR=w/(m+1)-mx
IndR=order(-mx,w)[1:r]
crk=NULL #to release memory
}
} else stop("Unknown method: ",method)
#computation of U
if(matrix.return) {
if (length(IndC)==1) {
C = matrix(A[,IndC],ncol=1) #submatrix must be recoerced from vector
} else {
C = A[,IndC]; # Choose c' columns from A.
}
if (length(IndR)==1) {
R = matrix(A[IndR,],nrow=1) #submatrix must be recoerced from vector
} else {
R = A[IndR,] ; # Choose r' rows from A.
}
# Compute U.
if (c==n) {
U = MASS::ginv(R)
} else {
if (r==m) {
U = MASS::ginv(C)
} else {
U = MASS::ginv(C) %*% A %*% MASS::ginv(R) ;
}
}
if (error.return) {
if (r<c) {
Error=Matrix::norm(A - C %*% U %*% R,"F")
} else {
Error=Matrix::norm(A - C %*% (U %*% R),"F")
}
}
} else {
U=NULL
C=NULL
R=NULL
}
res = new(Class='CURobj')
#result=list(C=C, U=U, R=R, IndC=IndC, IndR=IndR, ScoresC=ScoresC, ScoresR=ScoresR, Error=Error)
if (!is.null(C)) res@C = C
if (!is.null(ScoresC)) res@C.leverage.score = ScoresC
if (!is.null(IndC)) res@C.index =IndC
if (!is.null(R)) res@R = R
if (!is.null(ScoresR)) res@R.leverage.score = ScoresR
if (!is.null(IndR)) res@R.index =IndR
if (!is.null(U)) res@U = U
if (!is.null(Error)) res@Error = Error
return(res)
}
setGeneric(
name= "leverage",
def=function(object, C=TRUE){
standardGeneric("leverage")
}
)
setMethod(
f="leverage",
signature= "CURobj",
definition=function(object, C=TRUE){
if (C==FALSE){
res = object@R.leverage.score
} else if (C==TRUE){
res = object@C.leverage.score
}
return(res)
}
)
setGeneric(
name= "plotLeverage",
def=function(x, C=TRUE, mplr=1000, top.n=100, top.col='red', top.pch=16, ul.col='red', ul.lty=2, ...){
standardGeneric("plotLeverage")
}
)
setMethod(
f="plotLeverage",
signature= "CURobj",
definition=function(x, C=TRUE, mplr=1000, top.n=100, top.col='red', top.pch=16, ul.col='red', ul.lty=2, ...){
if (mplr!=1 && mplr!=0){
lab=paste('Leverage score (*', mplr, ')', sep='')
l.s = leverage(x, C=C)*mplr
} else {
lab = 'Leverage score'
l.s = leverage(x, C=C)
}
plot(l.s, type='h', ylab=lab, ...)
idx = rev(sort.list(l.s))[1:top.n]
points(idx, l.s[idx], pch=top.pch, col=top.col)
abline(h=mean(l.s), col=ul.col, lty=ul.lty)
}
)
setGeneric(
name= "topLeverage",
def=function(object, C=TRUE, top.n=100, sort=TRUE){
standardGeneric("topLeverage")
}
)
setMethod(
f="topLeverage",
signature= "CURobj",
definition=function(object, C=TRUE, top.n=100, sort=TRUE){
l.s = leverage(object, C=C)
idx = rev(sort.list(l.s))[1:top.n]
if (sort==TRUE){
idx = sort(idx)
}
return(idx)
}
)
setGeneric(
name= "getR",
def=function(object){
standardGeneric("getR")
}
)
setMethod(
f="getR",
signature= "CURobj",
definition=function(object){
return(object@R)
}
)
setGeneric(
name= "getU",
def=function(object){
standardGeneric("getU")
}
)
setMethod(
f="getU",
signature= "CURobj",
definition=function(object){
return(object@U)
}
)
setGeneric(
name= "getC",
def=function(object){
standardGeneric("getC")
}
)
setMethod(
f="getC",
signature= "CURobj",
definition=function(object){
return(object@C)
}
)
setGeneric(
name= "getError",
def=function(object){
standardGeneric("getError")
}
)
setMethod(
f="getError",
signature= "CURobj",
definition=function(object){
return(object@Error)
}
)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.