cholreduce.r

这是核学习的一个基础软件包

setClass("inc.chol",representation("matrix",pivots="vector",diag.residues="vector",maxresiduals="vector"),prototype=structure(.Data=matrix(),pivots=vector(),diag.residues=vector(), maxresiduals=vector()))


chol.reduce <- function(x, kernel="rbfdot",kpar=list(sigma=0.1), tol= 0.001, max.iter = dim(x)[1], verbose = 0)
{
## Usage: 
## 
## [ Tk, pivots, diag.residues, maxresiduals ] ... 
## <- chol.reduce( data, kernel, tol, max.iter, verbose )
## 
## Description:
## 
## Find the incomplete Cholesky decomposition of the kernel matrix
## 
## Parameters:
## 
## data      : 
## kernel    : kernlab object
## tol       : algo stops when remaining pivots < tol 
## max.iter  : maximum number of colums in Tk
##
## Return:
## 
## Tk     : K \approx Tk * Tk'
## pivots : Indices on which we pivoted
## diag.residues : Residuals left on the diagonal
## maxresiduals  : Residuals we picked for pivoting
## 
## Authors : S.V.N. Vishwanathan / Alex Smola
## R Version : Alexandros Karatzoglou

## For aggressive memory allocation
  BLOCKSIZE <- 50

if(!is.matrix(x))
	stop("x must be a matrix")

if(!is(kernel,"kernel"))
    {
      if(is(kernel,"function")) kernel <- deparse(substitute(kernel))
        kernel <- do.call(kernel, kpar)
    }

     if(!is(kernel,"kernel")) stop("kernel must inherit from class `kernel'")
## Begin initialization

  Tk <- T <- pivots <- maxresiduals <- padded.veck <- matrix(0,0,0)
  counter <- 0

## End initialization

  m <- dim(x)[1] 

## Compute the diagonal of kernel matrix

  diag.residues <- matrix(0, m, 1)

  for (i  in  1:m) 
    diag.residues[i] <- kernel(x[i,],x[i,])
  
## Choose first pivot

  residue <- max(diag.residues)
  index <- which.max(diag.residues == residue)

  dot.squarex <- t(rowSums(x^2))

  while( residue > tol && counter < max.iter )
    {
      ## Aggressively allocate memory
      if(counter %% BLOCKSIZE  == 0)
        {
          Tktmp <- matrix(0, m, dim(Tk)[2] + BLOCKSIZE)
          Tktmp[1:m > 0, 1:(dim(Tk)[2] + BLOCKSIZE) <= dim(Tk)[2]] <- Tk
          Tk <- Tktmp

          Ttmp <- matrix(0, dim(T)[1]+BLOCKSIZE, BLOCKSIZE+counter)
          ind <- 1:(dim(T)[1]+BLOCKSIZE) <= dim(T)[1]
          ind2 <- 1:(BLOCKSIZE + counter) <= counter
          Ttmp[ind , ind2] <- T
          Ttmp[ind == FALSE, ind2 == FALSE] <- diag(1, BLOCKSIZE)
          T <- Ttmp
          
          padded.veck.tmp <- matrix(0,dim(padded.veck)[1]+BLOCKSIZE)
          padded.veck.tmp[1:(dim(padded.veck)[1]+BLOCKSIZE) <= dim(padded.veck)[1]] <- padded.veck
          padded.veck <-  padded.veck.tmp
    
          pivots.tmp <- matrix(0, dim(pivots)[1]+BLOCKSIZE)
          pivots.tmp[1:(dim(pivots)[1] + BLOCKSIZE)<= dim(pivots)[1]] <- pivots
          pivots <- pivots.tmp
     
          maxresiduals.tmp <- matrix(0,dim(maxresiduals)[1]+BLOCKSIZE)
          maxresiduals.tmp[1:(dim(maxresiduals)[1]+BLOCKSIZE) <= dim(maxresiduals)[1]] <- maxresiduals
          maxresiduals <- maxresiduals.tmp
        } 
    
      veck <- kernelMatrix(kernel, x, x[index, ,drop=FALSE])
  
      if (counter == 0)
        {
          ## No need to compute t here
          tau <- sqrt(veck[index])
          
          ## Update T
          T[1, 1] <- tau
          ## Compute the update for Tk
          update <- veck/tau
        }
      else
        {
          padded.veck[1:counter] <- veck[pivots[1:counter]]
          
          ## First compute t
          t <- t(crossprod(padded.veck,backsolve(T,diag(rep(1,dim(T)[1])))))
          
          ## Now compute tau
          tau <- as.vector(sqrt(veck[index] - crossprod(t)))

          ## Update T

          T[1:counter, counter+1] <- t[1:counter]
          T[counter + 1, counter + 1] <- tau

          ## Compute the update for Tk
          update <- (1/tau) * (veck - Tk %*% t)
        }  

      ## Update Tk
      Tk[,counter + 1] <- update
      
      ## Update diagonal residuals
      diag.residues <- diag.residues - update^2			
      
      ## Update pivots
      pivots[counter + 1]  <- index
      
      ## Monitor residuals
      maxresiduals[counter + 1] <- residue 
      
      ## Choose next candidate
      residue <- max( diag.residues )
      index <- which.max(diag.residues)

      ## Update counter
      counter <- counter + 1
      
      ## Report progress to the user
      if(counter%%50 == 0 && (verbose == TRUE))
        cat("counter = ",counter," ", "residue = ", residue, "\n")
    } 

  ## Throw away extra columns which we might have added
  Tk <- Tk[, 1:counter] 
  
  pivots <- pivots[1:counter]
  
  maxresiduals <- maxresiduals[1:counter]

  return(new("inc.chol",.Data=Tk,pivots=pivots,diag.residues = diag.residues, maxresiduals = maxresiduals))
  
}