Arquivos
2014-11-21 18:52:56 +05:30

2009 linhas
83 KiB
Python

# Copyright (c) 2006, National ICT Australia
# All rights reserved.
#
# The contents of this file are subject to the Mozilla Public License Version
# 1.1 (the 'License'); you may not use this file except in compliance with
# the License. You may obtain a copy of the License at
# http://www.mozilla.org/MPL/
#
# Software distributed under the License is distributed on an 'AS IS' basis,
# WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
# for the specific language governing rights and limitations under the
# License.
#
# Authors: Le Song (lesong@it.usyd.edu.au) and Alex Smola
# (alex.smola@nicta.com.au)
# Created: (20/10/2006)
# Last Updated: (dd/mm/yyyy)
#
##\package elefant.kernels.vector
# This module contains kernel classes for vectorial data.
#
# The CVectorKernel class provides common interface for all kernel classes
# for vectorial data. This kernel should never be instantiated as well.
# All other kernel classes are derived from the CVectorKernel class. The
# hierarchy for the the classes in this package is as follows:
# --CKernel (abstract)
# ------CVectorKernel (abstract)
# ----------CDeltaKernel
# ----------CLinearKernel
# ----------CDotProductKernel (abstract)
# --------------CPolynomialKernel
# --------------CTanhKernel
# --------------CExpKernel
# ----------CRBFKernel (abstract)
# --------------CGaussKernel
# --------------CLaplaceKernel
# --------------CInvDisKernel
# --------------CInvSqDisKernel
# --------------CBesselKernel
#
__version__ = '$Revision: $'
# $Source$
import numpy
import numpy.random as random
from generic import CKernel
from my_exceptions import CElefantConstraintException
## Generic kernel class for vectorial data
#
# This kernel provide common interface for all kernels operating on
# vectorial data. This interface includes the following key kernel
# manipulations (functions):
# --Dot(x1, x2): $K(x1, x2)$
# --Expand(x1, x2, alpha): $sum_r K(x1_i,x2_r) \times alpha2_r$
# --Tensor(x1, y1, x2, y2): $K(x1_i,x2_j) \times (y1_i \times y1_j)$
# --TensorExpand(x1, y1, x2, y2, alpha2):
# $sum_r K(x1_i,x2_r) \times (y1_i \times y1_r) \times alpha2_r$
# --Remember(x): Remember data x
# --Forget(x): Remove remembered data x
# To design a specific kernel, simply overload these methods. The generic
# kernel itself should never be instantiated, although methods
# Remember and Forget are implemented in this class. The Remember method
# stores the inner product of an input vector.
#
class CVectorKernel(CKernel):
def __init__( self, blocksize = 128 ):
CKernel.__init__( self, blocksize )
self._name = 'Vector kernel'
## @var _cacheKernel
# Cache that store the base part for the kernel matrix.
# This cache facilates the incremental and decremental
# computational of the kernel matrix.
#
self._cacheKernel = {}
## @var _typicalParam
# Typical parameter for the kernel. Many kernels do not
# have parameters, such linear kernel. In this case, set
# zero as the typical parameter. This variable will be
# usefull when optimizing the kernel matrix with respect
# to the kernel matrix.
#
self._typicalParam = numpy.array([0])
## Compute the kernel between two data points x1 and x2.
# It returns a scale value of dot product between x1 and x2.
# @param x1 [read] The first data point.
# @param x2 [read] The second data point.
#
def K(self, x1, x2):
raise NotImplementedError, \
'CVectorKernel.K in abstract class is not implemented'
## Compute the kernel between the data points in x1 and those in x2.
# It returns a matrix with entry $(ij)$ equal to $K(x1_i, x1_j)$.
# If index1/index2 is
# specified, only those data points in x1/x2 with indices corresponding
# to index1/index2 are used to compute the kernel matrix. Furthermore,
# if output is specified, the provided buffer is used explicitly to
# store the kernel matrix.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Dot(self, x1, x2, index1=None, index2=None, output=None):
raise NotImplementedError, \
'CVectorKernel.Dot in abstract class is not implemented'
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix by alpha2.
# It returns a matrix with entry $(ij)$ equal to
# $sum_r K(x1_i,x2_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param alpha2 [read] The set of coefficients.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Expand(self, x1, x2, alpha2, index1=None, index2=None):
raise NotImplementedError, \
'CVectorKernel.Expand in abstract class is not implemented'
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2. It returns a matrix
# with entry $(ij)$ equal to $K(x1_i,x2_j) \times (y1_i \times y1_j)$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Tensor(self, x1, y1, x2, y2, index1=None, index2=None):
raise NotImplementedError, \
'CVectorKernel.Tensor in abstract class is not implemented'
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2, and final multiply
# the resulting matrix by alpha2. It returns a matrix with entry $(ij)$
# equal to $sum_r K(x1_i,x2_r) \times (y1_i \times y1_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def TensorExpand(self, x1, y1, x2, y2, alpha2, index1=None, index2=None, \
output=None):
raise NotImplementedError, \
'CVectorKernel.K in abstract class is not implemented'
## Remember x by computing the inner product of the data points contained
# in x, storing them in the cache and indexing them by the id of
# x. If x have already been remembered,
# the old stored information is simply overwritten.
# @param x [read] The data to be remembered.
#
def Remember(self, x):
# default behavior
assert x is not None, 'x is None'
assert len(x.shape) == 2, 'x is not a matrix'
self._cacheData[id(x)] = (x**2).sum(axis=1)
## Remove the remembered data x. If x is not given, then all remembered
# data in the cache is removed. If x has never been remembered, then
# this function does nothing and return False instead.
# @param x [read] The data to be removed.
#
def Forget(self, x=None):
# default behavior
if x is not None:
assert len(x.shape) == 2, 'Argument 1 is has wrong shape'
if self._cacheData.has_key(id(x)) is False:
return False
else:
del self._cacheData[id(x)]
else:
self._cacheData.clear()
return True
## Method that operates on the base part x of a kernel.
# The derived classes overload this method to generate new
# kernels.
# @param x Base part of the kernel.
#
def Kappa(self, x):
# default behavior
return x
## Gradient of the kernel with respect to the kernel
# parameter evaluated in the base part x of a kernel.
# The derived classes overload this method to generate new
# gradients of the kernels.
# @param x Base part of the kernel.
#
def KappaGrad(self, x):
# default behavior
return numpy.zeros(x.shape)
## Function that set the parameter of the kernel.
# If the derived classes have parameters, overload this
# method to set the parameters.
# @param param Parameters to be set.
#
def SetParam(self, param):
# default behavior
pass
## Clear the base part of the kernel computed for data x.
# If x is not given, then all remembered data in the cache is removed.
# If x has never been remembered, then this function does nothing
# and return False instead.
# @param x [read] The data whose base part is to be removed from the cache.
#
def ClearCacheKernel(self, x=None):
# default behavior
if x is not None:
assert len(x.shape) == 2, "Argument 1 has wrong shape"
if self._cacheKernel.has_key(id(x)) is False:
return False
else:
del self._cacheKernel[id(x)]
else:
self._cacheKernel.clear()
return True
## Create the cache for the base part of the kernel computed for
# data x, and index them by the id of x. If x have already been
# remembered, the old stored information is simply overwritten.
# Overload this method to store different base part for different
# kernels.
# @param x [read] The data whose base part is to be cached.
#
def CreateCacheKernel(self, x):
raise NotImplementedError, \
'CVectorKernel.K in abstract class is not implemented'
## Dot product of x with itself with the cached base part of the kernel.
# Overload this method to use the base part differently for different
# kernel. If param is given, the kernel matrix is computed using
# the given parameter and the current base part. Otherwise, the old
# parameters are used.
# @param x The data set.
# @param param The new parameters.
# @param output The output buffer.
#
def DotCacheKernel(self, x, param=None, output=None):
raise NotImplementedError, \
'CVectorKernel.K in abstract class is not implemented'
## Decrement the base part of the kernel for x1 stored in the cache
# by x2. Overload this method to define the decrement of the base
# part differently for different kernels. Note that this method
# updates the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
#
def DecCacheKernel(self, x1, x2):
raise NotImplementedError, \
'CVectorKernel.K in abstract class is not implemented'
## Decrement the base part of the kernel for x1 stored in the cache
# by x2, and return the resulting kernel matrix. If param is given,
# the kernel matrix is computed using the given parameter and the
# current base part. Otherwise, the old parameters are used. Overload
# this method to have different behavior for different kernel. Note
# that this method does NOT change the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
# @param param The new parameters.
#
def DecDotCacheKernel(self, x1, x2, param=None, output=None):
raise NotImplementedError, \
'CVectorKernel.K in abstract class is not implemented'
## Gradient of the kernel matrix with respect to the kernel parameter.
# If param is given, the kernel matrix is computed using the given
# parameter and the current base part. Otherwise, the old parameters
# are used. Overload this method to have different behavior.
# @param x The data set for the kernel matrix.
# @param param The kernel parameters.
#
def GradDotCacheKernel(self, x, param=None, output=None):
# default behavior
assert len(x.shape)==2, "Argument 1 has wrong shape"
assert self._cacheKernel.has_key(id(x)) == True, \
"Argument 1 has not been cached"
if param is not None:
self.SetParam(param)
n = x.shape[0]
output = numpy.zeros((n,n), numpy.float64)
return output
#------------------------------------------------------------------------------
## Linear kernel class
#
# The methods are implemente efficiently by computing the resulting
# matrices block by block.
#
class CLinearKernel(CVectorKernel):
def __init__(self, blocksize=128):
CVectorKernel.__init__(self, blocksize)
self._name = 'Linear kernel'
## Compute the kernel between two data points x1 and x2.
# It returns a scale value of dot product between x1 and x2.
# @param x1 [read] The first data point.
# @param x2 [read] The second data point.
#
def K(self, x1, x2):
assert len(x1.squeeze().shape) == 1, 'x1 is not a vector'
assert len(x2.squeeze().shape) == 1, 'x2 is not a vector'
return (x1.squeeze()*x2.squeeze()).sum()
## Compute the kernel between the data points in x1 and those in x2.
# It returns a matrix with entry $(ij)$ equal to $K(x1_i, x1_j)$.
# If index1/index2 is
# specified, only those data points in x1/x2 with indices corresponding
# to index1/index2 are used to compute the kernel matrix. Furthermore,
# if output is specified, the provided buffer is used explicitly to
# store the kernel matrix.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Dot(self, x1, x2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and Argument 2 have different dimensions'
if index1 is not None:
x1 = x1[index1,]
if index2 is not None:
x2 = x2[index2,]
# number of data points in x1.
n1 = x1.shape[0]
# number of data points in x2.
n2 = x2.shape[0]
# number of blocks.
nb = n1 / self._blocksize
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.dot(x1, numpy.transpose(x2))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,])))
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix by alpha2.
# It returns a matrix with entry $(ij)$ equal to
# $sum_r K(x1_i,x2_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param alpha2 [read] The set of coefficients.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Expand(self, x1, x2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert len(alpha2.shape) == 2, 'Argument 3 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and 2 has different dimesions'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 2 and 3 has different number of data points'
if index1 is not None:
x1 = x1[index1,]
if index2 is not None:
x2 = x2[index2,]
alpha2 = alpha2[index2,]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.dot(numpy.dot(x1, numpy.transpose(x2)), alpha2)
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.dot(numpy.transpose(numpy.dot(x2,numpy.transpose(x1[lower_limit:upper_limit,]))), alpha2)
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.dot(numpy.transpose(numpy.dot(x2,numpy.transpose(x1[lower_limit:n1,]))), alpha2)
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2. It returns a matrix
# with entry $(ij)$ equal to $K(x1_i,x2_j) \times (y1_i \times y1_j)$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Tensor(self, x1, y1, x2, y2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(y1.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == y1.shape[0], \
'Argument 1 and 2 has different dimensions'
assert len(x2.shape) == 2, 'Argument 3 has wrong shape'
assert len(y2.shape) == 2, 'Argument 4 has wrong shape'
assert x2.shape[0] == y2.shape[0], \
'Argument 2 and 3 has different dimensions'
if index1 is not None:
x1 = x1[index1,]
y1 = y1[index1,]
if index2 is not None:
x2 = x2[index2,]
y2 = y2[index2,]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = x2.shape[0]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.transpose(y1[:,0]*numpy.transpose(y2[:,0]*numpy.dot(x1, numpy.transpose(x2))))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = y2[:,0]*numpy.transpose(y1[lower_limit:upper_limit,0]*numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = y2[:,0]*numpy.transpose(y1[lower_limit:n1,0]*numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,])))
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2, and final multiply
# the resulting matrix by alpha2. It returns a matrix with entry $(ij)$
# equal to $sum_r K(x1_i,x2_r) \times (y1_i \times y1_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def TensorExpand(self, x1, y1, x2, y2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(y1.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == y1.shape[0], \
'Argument 1 and 2 have different dimensions'
assert len(x2.shape) == 2, 'Argument 3 has wrong shape'
assert len(y2.shape) == 2, 'Argument 4 has wrong shape'
assert x2.shape[0] == y2.shape[0], \
'Argument 3 and 4 have different dimensions'
assert len(alpha2.shape) == 2, 'Argument 5 has wrong shape'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 3 and 5 have different number of data points'
if index1 is not None:
x1 = x1[index1,]
y1 = y1[index1,]
if index2 is not None:
x2 = x2[index2,]
y2 = y2[index2,]
alpha2 = alpha2[index2,]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.transpose(y1[:,0] * numpy.transpose(numpy.dot(y2[:,0]*numpy.dot(x1,numpy.transpose(x2)), alpha2)))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.dot(y2[:,0]*numpy.transpose(y1[lower_limit:upper_limit,0] * numpy.dot(x2,numpy.transpose(x1[lower_limit:upper_limit,]))), alpha2)
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.dot(y2[:,0]*numpy.transpose(y1[lower_limit:n1,0] * numpy.dot(x2,numpy.transpose(x1[lower_limit:n1,]))), alpha2)
return output
## Create the cache for the base part of the kernel computed for
# data x, and index them by the id of x. If x have already been
# remembered, the old stored information is simply overwritten.
# @param x [read] The data whose base part is to be cached.
#
def CreateCacheKernel(self, x):
assert len(x.shape) == 2, 'Argument 1 has wrong shape'
n = x.shape[0]
nb = n / self._blocksize
# create the cache space
if self._cacheKernel.has_key(id(x)):
self.ClearCacheKernel(x)
tmpCacheKernel = numpy.zeros((n,n), numpy.float64)
self._cacheKernel[id(x)] = tmpCacheKernel
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
tmpCacheKernel[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x, numpy.transpose(x[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n:
tmpCacheKernel[lower_limit:n,] = numpy.transpose(numpy.dot(x, numpy.transpose(x[lower_limit:n,])))
return True
## Dot product of x with itself with the cached base part of the kernel.
# If param is given, the kernel matrix is computed using
# the given parameter and the current base part. Otherwise, the old
# parameters are used.
# @param x The data set.
# @param param The new parameters.
# @param output The output buffer.
#
def DotCacheKernel(self, x, param=None, output=None):
assert len(x.shape)==2, 'Argument 1 has wrong shape'
assert self._cacheKernel.has_key(id(x)) == True, \
'Argument 1 has not been cached'
n = x.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x)]
# set parameters.
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = tmpCacheKernel[lower_limit:upper_limit,]
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = tmpCacheKernel[lower_limit:n,]
return output
## Decrement the base part of the kernel for x1 stored in the cache
# by x2. Note that this method updates the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
#
def DecCacheKernel(self, x1, x2):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and 2 have different number of data points'
assert self._cacheKernel.has_key(id(x1)) == True, \
'Argument 1 has not been cached'
n = x1.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x1)]
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
tmpCacheKernel[lower_limit:upper_limit,] = tmpCacheKernel[lower_limit:upper_limit,] \
- numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n:
tmpCacheKernel[lower_limit:n,] = tmpCacheKernel[lower_limit:n,] \
- numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:n,])))
return True
## Decrement the base part of the kernel for x1 stored in the cache
# by x2, and return the resulting kernel matrix. If param is given,
# the kernel matrix is computed using the given parameter and the
# current base part. Otherwise, the old parameters are used. Note
# that this method does NOT change the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
# @param param The new parameters.
#
def DecDotCacheKernel(self, x1, x2, param=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and 2 have different number of data points'
assert self._cacheKernel.has_key(id(x1)) == True, \
'Argument 1 has not been cached'
n = x1.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x1)]
# set parameters.
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = tmpCacheKernel[lower_limit:upper_limit,] \
- numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = tmpCacheKernel[lower_limit:n,] \
- numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:n,])))
return output
#------------------------------------------------------------------------------
## Dot Product kernel class
#
# All kernels that are function of the dot product between the two
# data points, ie. $K(x1,x2)=Kappa(x1^\top x2).
# The methods are implemente efficiently by computing the resulting
# matrices block by block. All kernel classes derived from CDotProductKernel
# differ only in their choice of Kappa function.
#
class CDotProductKernel(CVectorKernel):
def __init__(self, blocksize=128):
CVectorKernel.__init__(self, blocksize)
self._name = 'Dot Product kernel'
## Compute the kernel between two data points x1 and x2.
# It returns a scale value of dot product between x1 and x2.
# @param x1 [read] The first data point.
# @param x2 [read] The second data point.
#
def K(self, x1, x2):
assert len(x1.squeeze().shape) == 1, 'x1 is not a vector'
assert len(x2.squeeze().shape) == 1, 'x2 is not a vector'
return (x1.squeeze()*x2.squeeze()).sum()
## Compute the kernel between the data points in x1 and those in x2.
# It returns a matrix with entry $(ij)$ equal to $K(x1_i, x1_j)$.
# If index1/index2 is
# specified, only those data points in x1/x2 with indices corresponding
# to index1/index2 are used to compute the kernel matrix. Furthermore,
# if output is specified, the provided buffer is used explicitly to
# store the kernel matrix.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Dot(self, x1, x2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and Argument 2 have different dimensions'
if index1 is not None:
x1 = x1[index1,]
if index2 is not None:
x2 = x2[index2,]
# number of data points in x1.
n1 = x1.shape[0]
# number of data points in x2.
n2 = x2.shape[0]
# number of blocks.
nb = n1 / self._blocksize
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = self.Kappa(numpy.dot(x1, numpy.transpose(x2)))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = self.Kappa(numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,]))))
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = self.Kappa(numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,]))))
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix by alpha2.
# It returns a matrix with entry $(ij)$ equal to
# $sum_r K(x1_i,x2_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param alpha2 [read] The set of coefficients.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Expand(self, x1, x2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert len(alpha2.shape) == 2, 'Argument 3 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and 2 has different dimesions'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 2 and 3 has different number of data points'
if index1 is not None:
x1 = x1[index1,]
if index2 is not None:
x2 = x2[index2,]
alpha2 = alpha2[index2,]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is not None:
output = output
else:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.dot(self.Kappa(numpy.dot(x1, numpy.transpose(x2))), alpha2)
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.dot(self.Kappa(numpy.transpose(numpy.dot(x2,numpy.transpose(x1[lower_limit:upper_limit,])))), alpha2)
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.dot(self.Kappa(numpy.transpose(numpy.dot(x2,numpy.transpose(x1[lower_limit:n1,])))), alpha2)
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2. It returns a matrix
# with entry $(ij)$ equal to $K(x1_i,x2_j) \times (y1_i \times y1_j)$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Tensor(self, x1, y1, x2, y2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(y1.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == y1.shape[0], \
'Argument 1 and 2 has different dimensions'
assert len(x2.shape) == 2, 'Argument 3 has wrong shape'
assert len(y2.shape) == 2, 'Argument 4 has wrong shape'
assert x2.shape[0] == y2.shape[0], \
'Argument 2 and 3 has different dimensions'
if index1 is not None:
x1 = x1[index1,]
y1 = y1[index1,]
if index2 is not None:
x2 = x2[index2,]
y2 = y2[index2,]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = x2.shape[0]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.transpose(y1[:,0]*numpy.transpose(y2[:,0]*self.Kappa(numpy.dot(x1, numpy.transpose(x2)))))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = y2[:,0]*numpy.transpose(y1[lower_limit:upper_limit,0]*self.Kappa(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,]))))
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = y2[:,0]*numpy.transpose(y1[lower_limit:n1,0]*self.Kappa(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,]))))
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2, and final multiply
# the resulting matrix by alpha2. It returns a matrix with entry $(ij)$
# equal to $sum_r K(x1_i,x2_r) \times (y1_i \times y1_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def TensorExpand(self, x1, y1, x2, y2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(y1.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == y1.shape[0], \
'Argument 1 and 2 have different dimensions'
assert len(x2.shape) == 2, 'Argument 3 has wrong shape'
assert len(y2.shape) == 2, 'Argument 4 has wrong shape'
assert x2.shape[0] == y2.shape[0], \
'Argument 3 and 4 have different dimensions'
assert len(alpha2.shape) == 2, 'Argument 5 has wrong shape'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 3 and 5 have different number of data points'
if index1 is not None:
x1 = x1[index1,]
y1 = y1[index1,]
if index2 is not None:
x2 = x2[index2,]
y2 = y2[index2,]
alpha2 = alpha2[index2,]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.transpose(y1[:,0] * numpy.transpose(numpy.dot(y2[:,0]*self.Kappa(numpy.dot(x1,numpy.transpose(x2))), alpha2)))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.dot(y2[:,0]*numpy.transpose(y1[lower_limit:upper_limit,0] * self.Kappa(numpy.dot(x2,numpy.transpose(x1[lower_limit:upper_limit,])))), alpha2)
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.dot(y2[:,0]*numpy.transpose(y1[lower_limit:n1,0] * self.Kappa(numpy.dot(x2,numpy.transpose(x1[lower_limit:n1,])))), alpha2)
return output
## Create the cache for the base part of the kernel computed for
# data x, and index them by the id of x. If x have already been
# remembered, the old stored information is simply overwritten.
# @param x [read] The data whose base part is to be cached.
#
def CreateCacheKernel(self, x):
assert len(x.shape) == 2, 'Argument 1 has wrong shape'
n = x.shape[0]
nb = n / self._blocksize
# create the cache space
if self._cacheKernel.has_key(id(x)):
self.ClearCacheKernel(x)
tmpCacheKernel = numpy.zeros((n,n), numpy.float64)
self._cacheKernel[id(x)] = tmpCacheKernel
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
tmpCacheKernel[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x, numpy.transpose(x[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n:
tmpCacheKernel[lower_limit:n,] = numpy.transpose(numpy.dot(x, numpy.transpose(x[lower_limit:n,])))
return True
## Dot product of x with itself with the cached base part of the kernel.
# If param is given, the kernel matrix is computed using
# the given parameter and the current base part. Otherwise, the old
# parameters are used.
# @param x The data set.
# @param param The new parameters.
# @param output The output buffer.
#
def DotCacheKernel(self, x, param=None, output=None):
assert len(x.shape)==2, 'Argument 1 has wrong shape'
assert self._cacheKernel.has_key(id(x)) == True, \
'Argument 1 has not been cached'
n = x.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x)]
# set parameters.
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = self.Kappa(tmpCacheKernel[lower_limit:upper_limit,])
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = self.Kappa(tmpCacheKernel[lower_limit:n,])
return output
## Decrement the base part of the kernel for x1 stored in the cache
# by x2. Note that this method updates the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
#
def DecCacheKernel(self, x1, x2):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and 2 have different number of data points'
assert self._cacheKernel.has_key(id(x1)) == True, \
'Argument 1 has not been cached'
n = x1.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x1)]
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
tmpCacheKernel[lower_limit:upper_limit,] = tmpCacheKernel[lower_limit:upper_limit,] \
- numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:upper_limit,])))
lower_limit = upper_limit
if lower_limit <= n:
tmpCacheKernel[lower_limit:n,] = tmpCacheKernel[lower_limit:n,] \
- numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:n,])))
return True
## Decrement the base part of the kernel for x1 stored in the cache
# by x2, and return the resulting kernel matrix. If param is given,
# the kernel matrix is computed using the given parameter and the
# current base part. Otherwise, the old parameters are used. Note
# that this method does NOT change the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
# @param param The new parameters.
#
def DecDotCacheKernel(self, x1, x2, param=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and 2 have different number of data points'
assert self._cacheKernel.has_key(id(x1)) == True, \
'Argument 1 has not been cached'
n = x1.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x1)]
# set parameters.
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:upper_limit,])))
output[lower_limit:upper_limit,] = self.Kappa(tmpCacheKernel[lower_limit:upper_limit,] \
- output[lower_limit:upper_limit,])
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:n,])))
output[lower_limit:n,] = self.Kappa(tmpCacheKernel[lower_limit:n,] \
- output[lower_limit:n,])
return output
## Gradient of the kernel matrix with respect to the kernel parameter.
# If param is given, the kernel matrix is computed using the given
# parameter and the current base part. Otherwise, the old parameters
# are used.
# @param x The data set for the kernel matrix.
# @param param The kernel parameters.
#
def GradDotCacheKernel(self, x, param=None, output=None):
assert len(x.shape) == 2, "Argument 1 has wrong shape"
assert self._cacheKernel.has_key(id(x)) == True, \
"Argument 1 has not been cached"
n = x.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x)]
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = self.KappaGrad(tmpCacheKernel[lower_limit:upper_limit,])
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = self.KappaGrad(tmpCacheKernel[lower_limit:n,])
return output
## Polynomial kernel
#
# Kernel of the form: $K(x1,x2)=(scale (x1^\top x2) + offset)^degree$.
# This is implemented by the Kappa function.
#
class CPolynomialKernel(CDotProductKernel):
def __init__(self, degree=2, offset=1.0, scale=1.0, blocksize=128):
CDotProductKernel.__init__(self, blocksize)
self._name = 'Polynomial kernel'
self._scale = scale
self._offset = offset
self._degree = degree
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return (self._scale * x + self._offset) ** self._degree
## Tanh kernel
#
# Kernel of the form : $K(x1,x2)=\tanh(scale (x1^\top x2) + offset)$.
# This is implemented by the Kappa function.
#
class CTanhKernel(CDotProductKernel):
def __init__(self, offset=1.0, scale=1.0, blocksize=128):
CDotProductKernel.__init__(self, blocksize)
self._name = 'Tahn kernel'
self._scale = scale
self._offset = offset
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return numpy.tanh(self._scale * x + self._offset)
## Exponential kernel
#
# Kernel of the form: $K(x1,x2)=\exp(scale (x1^\top x2) + offset)
# This is implemented by the Kappa function.
#
class CExpKernel(CDotProductKernel):
def __init__(self, offset=0, scale=1.0, blocksize=128):
CDotProductKernel.__init__(self, blocksize)
self._name = 'Expential kernel'
self._scale = scale
self._offset = offset
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return numpy.exp(self._scale * x + self._offset)
#------------------------------------------------------------------------------
## RBF kernels.
#
# Distance-based kernel, where Gauss, Laplace and Bessel kernels are derived.
# This class of kernels has the form: $K(x1,x2)=f(\|x1-x2\|)$
# The methods are implemente efficiently by computing the resulting
# matrices block by block. All kernel classes derived from CRBFKernel
# differ only in their choice of Kappa function.
#
class CRBFKernel(CVectorKernel):
def __init__(self, blocksize=128):
CVectorKernel.__init__(self, blocksize=blocksize)
self._name = 'RBF kernel'
## Method that operates on the base part of a kernel.
# The derived classes overload this method to generate new
# kernels.
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return x
## Compute the kernel between two data points x1 and x2.
# It returns a scale value of dot product between x1 and x2.
# @param x1 [read] The first data point.
# @param x2 [read] The second data point.
#
def K(self, x1, x2):
##assert len(x1.squeeze().shape) == 1, 'x1 is not a vector'
##assert len(x2.squeeze().shape) == 1, 'x2 is not a vector'
v = (x1 - x2).squeeze()
return self.Kappa(numpy.sum(v*v))
## Compute the square distance using inner products. It computes
# the distance by following formular:
# $(x1-x2)^2 = x1^2 - 2 x1 x2 + x2^2$
# @param dest [write] The destination where the results are stored.
# @param x1_x2 [read] The inner product matrix between x1 and x2.
# @param dot_x1 [read] The inner product of the data in x1.
# @param dot_x2 [read] The inner product of the data in x2.
#
def KappaSqDis(self, dest, x1_x2, dot_x1, dot_x2):
dest[:,] = self.Kappa(- 2.0 * x1_x2 + numpy.add.outer(dot_x1, dot_x2))
## Compute the kernel between the data points in x1 and those in x2.
# It returns a matrix with entry $(ij)$ equal to $K(x1_i, x1_j)$.
# If index1/index2 is
# specified, only those data points in x1/x2 with indices corresponding
# to index1/index2 are used to compute the kernel matrix. Furthermore,
# if output is specified, the provided buffer is used explicitly to
# store the kernel matrix.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Dot(self, x1, x2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and Argument 2 have different dimensions'
# retrieve remembered data from the cache.
flg1 = False
if self._cacheData.has_key(id(x1)):
dot_x1 = self._cacheData[id(x1)]
flg1 = True
flg2 = False
if self._cacheData.has_key(id(x2)):
dot_x2 = self._cacheData[id(x2)]
flg2 = True
if index1 is not None:
x1 = x1[index1,]
if flg1 is True:
dot_x1 = dot_x1[index1,]
else:
dot_x1 = numpy.sum(x1**2, 1)
else:
if flg1 is False:
dot_x1 = numpy.sum(x1**2, 1)
if index2 is not None:
x2 = x2[index2,]
if flg2 is True:
dot_x2 = dot_x2[index2,]
else:
dot_x2 = numpy.sum(x2**2, 1)
else:
if flg2 is False:
dot_x2 = numpy.sum(x2**2, 1)
# number of data points in x1.
n1 = x1.shape[0]
# number of data points in x2.
n2 = x2.shape[0]
# number of blocks.
nb = n1 / self._blocksize
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.dot(x1, numpy.transpose(x2))
self.KappaSqDis(output, output, dot_x1, dot_x2)
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,])))
self.KappaSqDis(output[lower_limit:upper_limit,], output[lower_limit:upper_limit,], dot_x1[lower_limit:upper_limit], dot_x2)
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,])))
self.KappaSqDis(output[lower_limit:n1,], output[lower_limit:n1,], dot_x1[lower_limit:n1], dot_x2)
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix by alpha2.
# It returns a matrix with entry $(ij)$ equal to
# $sum_r K(x1_i,x2_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param alpha2 [read] The set of coefficients.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Expand(self, x1, x2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert len(alpha2.shape) == 2, 'Argument 3 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and 2 has different dimesions'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 2 and 3 has different number of data points'
# retrieve remembered data from cache.
flg1 = False
if self._cacheData.has_key(id(x1)):
dot_x1 = self._cacheData[id(x1)]
flg1 = True
flg2 = False
if self._cacheData.has_key(id(x2)):
dot_x2 = self._cacheData[id(x2)]
flg2 = True
if index1 is not None:
x1 = x1[index1,]
if flg1 is True:
dot_x1 = dot_x1[index1,]
else:
dot_x1 = numpy.sum(x1**2, 1)
else:
if flg1 is False:
dot_x1 = numpy.sum(x1**2, 1)
if index2 is not None:
x2 = x2[index2,]
if flg2 is True:
dot_x2 = dot_x2[index2,]
else:
dot_x2 = numpy.sum(x2**2, 1)
alpha2 = alpha2[index2,]
else:
if flg2 is False:
dot_x2 = numpy.sum(x2**2, 1)
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases.
if index2 is not None:
if len(index2) <= self._blocksize:
x1_x2 = numpy.dot(x1, numpy.transpose(x2))
self.KappaSqDis(x1_x2, x1_x2, dot_x1, dot_x2)
output = numpy.dot(x1_x2,alpha2)
return output
# blocking.
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
x1_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,])))
self.KappaSqDis(x1_x2, x1_x2, dot_x1[lower_limit:upper_limit], dot_x2)
output[lower_limit:upper_limit,] = numpy.dot(x1_x2, alpha2)
lower_limit = upper_limit
if lower_limit <= n1:
x1_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,])))
self.KappaSqDis(x1_x2, x1_x2, dot_x1[lower_limit:n1], dot_x2)
output[lower_limit:n1,] = numpy.dot(x1_x2, alpha2)
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2. It returns a matrix
# with entry $(ij)$ equal to $K(x1_i,x2_j) \times (y1_i \times y1_j)$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Tensor(self, x1, y1, x2, y2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(y1.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == y1.shape[0], \
'Argument 1 and 2 has different dimensions'
assert len(x2.shape) == 2, 'Argument 3 has wrong shape'
assert len(y2.shape) == 2, 'Argument 4 has wrong shape'
assert x2.shape[0] == y2.shape[0], \
'Argument 2 and 3 has different dimensions'
# retrieve remembered data from the cache
flg1 = False
if self._cacheData.has_key(id(x1)):
dot_x1 = self._cacheData[id(x1)]
flg1 = True
flg2 = False
if self._cacheData.has_key(id(x2)):
dot_x2 = self._cacheData[id(x2)]
flg2 = True
if index1 is not None:
x1 = x1[index1,]
if flg1 is True:
dot_x1 = dot_x1[index1,]
else:
dot_x1 = numpy.sum(x1**2, 1)
y1 = y1[index1,]
else:
if flg1 is False:
dot_x1 = numpy.sum(x1**2, 1)
if index2 is not None:
x2 = x2[index2,]
if flg2 is True:
dot_x2 = dot_x2[index2,]
else:
dot_x2 = numpy.sum(x2**2, 1)
y2 = y2[index2,]
else:
if flg2 is False:
dot_x2 = numpy.sum(x2**2, 1)
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = x2.shape[0]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
x1_x2 = numpy.dot(x1, numpy.transpose(x2))
self.KappaSqDis(x1_x2, x1_x2, dot_x1, dot_x2)
output = numpy.transpose(y1[:,0] * numpy.transpose(y2[:,0] * x1_x2))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
x1_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,])))
self.KappaSqDis(x1_x2, x1_x2, dot_x1[lower_limit:upper_limit], dot_x2)
output[lower_limit:upper_limit,] = numpy.transpose(y1[lower_limit:upper_limit,0] * numpy.transpose(y2[:,0] * x1_x2))
lower_limit = upper_limit
if lower_limit <= n1:
x1_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,])))
self.KappaSqDis(x1_x2, x1_x2, dot_x1[lower_limit:n1], dot_x2)
output[lower_limit:n1,] = numpy.transpose(y1[lower_limit:n1,0] * numpy.transpose(y2[:,0] * x1_x2))
return output
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix elementwiesely by the
# the outer-product matrix between y1 and y2, and final multiply
# the resulting matrix by alpha2. It returns a matrix with entry $(ij)$
# equal to $sum_r K(x1_i,x2_r) \times (y1_i \times y1_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param y1 [read] The first set of labels.
# @param x2 [read] The second set of data points.
# @param y2 [read] The second set of labels.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def TensorExpand(self, x1, y1, x2, y2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(y1.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == y1.shape[0], \
'Argument 1 and 2 have different dimensions'
assert len(x2.shape) == 2, 'Argument 3 has wrong shape'
assert len(y2.shape) == 2, 'Argument 4 has wrong shape'
assert x2.shape[0] == y2.shape[0], \
'Argument 3 and 4 have different dimensions'
assert len(alpha2.shape) == 2, 'Argument 5 has wrong shape'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 3 and 5 have different number of data points'
# retrieve remembered data from the cache
flg1 = False
if self._cacheData.has_key(id(x1)):
dot_x1 = self._cacheData[id(x1)]
flg1 = True
flg2 = False
if self._cacheData.has_key(id(x2)):
dot_x2 = self._cacheData[id(x2)]
flg2 = True
if index1 is not None:
x1 = x1[index1,]
if flg1 is True:
dot_x1 = dot_x1[index1,]
else:
dot_x1 = numpy.sum(x1**2, 1)
y1 = y1[index1,]
else:
if flg1 is False:
dot_x1 = numpy.sum(x1**2, 1)
if index2 is not None:
x2 = x2[index2,]
if flg2 is True:
dot_x2 = dot_x2[index2,]
else:
dot_x2 = numpy.sum(x2**2, 1)
y2 = y2[index2,]
alpha2 = alpha2[index2,]
else:
if flg2 is False:
dot_x2 = numpy.sum(x2**2, 1)
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
x1_x2 = numpy.dot(x1, numpy.transpose(x2))
self.KappaSqDis(x1_x2, x1_x2, dot_x1, dot_x2)
output = numpy.transpose(y1[:,0] * numpy.transpose(numpy.dot(y2[:,0]*x1_x2, alpha2)))
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
x1_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:upper_limit,])))
self.KappaSqDis(x1_x2, x1_x2, dot_x1[lower_limit:upper_limit], dot_x2)
output[lower_limit:upper_limit,] = numpy.transpose(y1[lower_limit:upper_limit,0] * numpy.transpose(numpy.dot(y2[:,0]*x1_x2, alpha2)))
lower_limit = upper_limit
if lower_limit <= n1:
x1_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x1[lower_limit:n1,])))
self.KappaSqDis(x1_x2, x1_x2, dot_x1[lower_limit:n1], dot_x2)
output[lower_limit:n1,] = numpy.transpose(y1[lower_limit:n1,0] * numpy.transpose(numpy.dot(y2[:,0]*x1_x2, alpha2)))
return output
## Create the cache for the base part of the kernel computed for
# data x, and index them by the id of x. If x have already been
# remembered, the old stored information is simply overwritten.
# @param x [read] The data whose base part is to be cached.
#
def CreateCacheKernel(self, x):
assert len(x.shape) == 2, 'Argument 1 has wrong shape'
n = x.shape[0]
nb = n / self._blocksize
# create the cache space
if self._cacheKernel.has_key(id(x)):
self.ClearCacheKernel(x)
tmpCacheKernel = numpy.zeros((n,n), numpy.float64)
self._cacheKernel[id(x)] = tmpCacheKernel
if self._cacheData.has_key(id(x)):
dot_x = self._cacheData[id(x)]
else:
dot_x = numpy.sum(x**2,1)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
tmpCacheKernel[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x, numpy.transpose(x[lower_limit:upper_limit,])))
tmpCacheKernel[lower_limit:upper_limit,] = numpy.add.outer(dot_x[lower_limit:upper_limit], dot_x) \
- 2*tmpCacheKernel[lower_limit:upper_limit,]
lower_limit = upper_limit
if lower_limit <= n:
tmpCacheKernel[lower_limit:n,] = numpy.transpose(numpy.dot(x, numpy.transpose(x[lower_limit:n,])))
tmpCacheKernel[lower_limit:n,] = numpy.add.outer(dot_x[lower_limit:n], dot_x) \
- 2*tmpCacheKernel[lower_limit:n,]
return True
## Dot product of x with itself with the cached base part of the kernel.
# If param is given, the kernel matrix is computed using
# the given parameter and the current base part. Otherwise, the old
# parameters are used.
# @param x The data set.
# @param param The new parameters.
# @param output The output buffer.
#
def DotCacheKernel(self, x, param=None, output=None):
assert len(x.shape)==2, 'Argument 1 has wrong shape'
assert self._cacheKernel.has_key(id(x)) == True, \
'Argument 1 has not been cached'
n = x.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x)]
# set parameters.
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = self.Kappa(tmpCacheKernel[lower_limit:upper_limit,])
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = self.Kappa(tmpCacheKernel[lower_limit:n,])
return output
## Decrement the base part of the kernel for x1 stored in the cache
# by x2. Note that this method updates the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
#
def DecCacheKernel(self, x1, x2):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and 2 have different number of data points'
assert self._cacheKernel.has_key(id(x1)) == True, \
'Argument 1 has not been cached'
n = x1.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x1)]
if self._cacheData.has_key(id(x2)):
dot_x2 = self._cacheData[id(x2)]
else:
dot_x2 = numpy.sum(x2**2,1)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
x2_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:upper_limit,])))
tmpCacheKernel[lower_limit:upper_limit,] = tmpCacheKernel[lower_limit:upper_limit,] \
+ 2*x2_x2 \
- numpy.add.outer(dot_x2[lower_limit:upper_limit], dot_x2)
lower_limit = upper_limit
if lower_limit <= n:
x2_x2 = numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:n,])))
tmpCacheKernel[lower_limit:n,] = tmpCacheKernel[lower_limit:n,] \
+ 2*x2_x2 \
- numpy.add.outer(dot_x2[lower_limit:n], dot_x2)
return True
## Decrement the base part of the kernel for x1 stored in the cache
# by x2, and return the resulting kernel matrix. If param is given,
# the kernel matrix is computed using the given parameter and the
# current base part. Otherwise, the old parameters are used. Note
# that this method does NOT change the cache for the kernel part.
# @param x1 The data set whose base part has been cached.
# @param x2 The data set who is to be decremented from x1.
# @param param The new parameters.
#
def DecDotCacheKernel(self, x1, x2, param=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and 2 have different number of data points'
assert self._cacheKernel.has_key(id(x1)) == True, \
'Argument 1 has not been cached'
n = x1.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x1)]
# set parameters.
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
if self._cacheData.has_key(id(x2)):
dot_x2 = self._cacheData[id(x2)]
else:
dot_x2 = numpy.sum(x2*x2,1)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:upper_limit,])))
output[lower_limit:upper_limit,] = self.Kappa(tmpCacheKernel[lower_limit:upper_limit,] \
+ 2*output[lower_limit:upper_limit,] \
- numpy.add.outer(dot_x2[lower_limit:upper_limit], dot_x2))
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = numpy.transpose(numpy.dot(x2, numpy.transpose(x2[lower_limit:n,])))
output[lower_limit:n,] = self.Kappa(tmpCacheKernel[lower_limit:n,] \
+ 2*output[lower_limit:n,] \
- numpy.add.outer(dot_x2[lower_limit:n], dot_x2))
return output
## Gradient of the kernel matrix with respect to the kernel parameter.
# If param is given, the kernel matrix is computed using the given
# parameter and the current base part. Otherwise, the old parameters
# are used.
# @param x The data set for the kernel matrix.
# @param param The kernel parameters.
#
def GradDotCacheKernel(self, x, param=None, output=None):
assert len(x.shape) == 2, "Argument 1 has wrong shape"
assert self._cacheKernel.has_key(id(x)) == True, \
"Argument 1 has not been cached"
n = x.shape[0]
nb = n / self._blocksize
tmpCacheKernel = self._cacheKernel[id(x)]
if param is not None:
self.SetParam(param)
if output is None:
output = numpy.zeros((n,n), numpy.float64)
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = self.KappaGrad(tmpCacheKernel[lower_limit:upper_limit,])
lower_limit = upper_limit
if lower_limit <= n:
output[lower_limit:n,] = self.KappaGrad(tmpCacheKernel[lower_limit:n,])
return output
## Gauss kernel
#
# Kernel of the form: $K(x1,x2)=\exp(-scale*\|x1-x2\|^2)$.
# This is implemented by the Kappa function.
#
class CGaussKernel(CRBFKernel):
def __init__(self, scale=1.0, blocksize=128):
CRBFKernel.__init__(self, blocksize)
self._name = 'Gauss kernel'
if scale < 0.0:
raise CElefantConstraintException(scale, "param must not be negative")
self._scale = scale
self._typicalParam = 10**(numpy.arange(10)-7.0)
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return numpy.exp(-self._scale * x)
## Gradient of the kernel with respect to the kernel
# parameter evaluated at the base part x of the kernel.
# @param x Base part of the kernel.
#
def KappaGrad(self, x):
return -x * numpy.exp(-self._scale * x)
## Function that set the parameter of the kernel.
# @param param Parameters to be set.
#
def SetParam(self, param):
if param < 0.0:
raise CElefantConstraintException(param, "param must not be negative")
self._scale = param
## Laplace Kernel
#
# Kernel of the form: $K(x1,x2)=\exp(-scale*\|x1-x2\|)$.
# This is implemented by the Kappa function.
#
class CLaplaceKernel(CRBFKernel):
def __init__(self, scale=1.0, blocksize=128):
CRBFKernel.__init__(self, blocksize)
self._name = 'Laplace kernel'
if scale < 0.0:
raise CElefantConstraintException(scale, "param must not be negative")
self._scale = scale
self._typicalParam = 10**(numpy.arange(10)-7.0)
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return numpy.exp(-self._scale * numpy.sqrt(x.clip(min = 0.0, max = numpy.inf)))
## Gradient of the kernel with respect to the kernel
# parameter evaluated at the base part x of the kernel.
# @param x Base part of the kernel.
#
def KappaGrad(self, x):
tmp = numpy.sqrt(x.clip(min = 0.0, max = numpy.inf))
return -tmp * numpy.exp(-self._scale * tmp)
## Function that set the parameter of the kernel.
# @param param Parameters to be set.
#
def SetParam(self, param):
# Enforce nonnegativity of parameter
# assert param >= 0.0, 'param must not be negative'
if param < 0.0:
raise CElefantConstraintException(param, "param must not be negative")
self._scale = param
## Inverse square distance kernel
#
# Kernel of the form: $K(x1,x2)=\frac{1}{\|x1-x2\|^2}
# This is implemented by the Kappa function. Note that
# if the distance between x1 and x2 is very small, the
# kernel is computed as zero.
#
class CInvSqDisKernel(CRBFKernel):
def __init__(self, blocksize=128, epsilon = 1.0e-6):
CRBFKernel.__init__(self, blocksize)
self._name = 'Inverse Square Distance kernel'
self._epsilon = epsilon
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return 1.0 / (x + self._epsilon)
## Inverse distance kernel
#
# Kernel of the form: $K(x1,x2)=\frac{1}{\|x1-x2\|}
# This is implemented by the Kappa function. Note that
# if the distance between x1 and x2 is very small, the
# kernel is computed as zero.
#
class CInvDisKernel(CRBFKernel):
def __init__(self, blocksize=128, epsilon = 1.0e-6):
CRBFKernel.__init__(self, blocksize)
self._name = 'Inverse Distance kernel'
self._epsilon = epsilon
## Method that operates on the base part of the kernel to
# generate the new kernel
# @param x Base part of the kernel matrix.
#
def Kappa(self, x):
return 1.0/numpy.sqrt(x + self._epsilon)
#------------------------------------------------------------------------------
## Delta kernel class
#
# The methods are implemente efficiently by computing the resulting
# matrices block by block. Currently this class only works for the
# case where the inputs are integer.
#
class CDeltaKernel(CVectorKernel):
def __init__(self, blocksize=128):
CVectorKernel.__init__(self, blocksize)
self._name = 'Delta kernel'
## Compute the kernel between two data points x1 and x2.
# It returns a value 1 if x1 and x2 are equal, otherwise 0.
# @param x1 [read] The first data point.
# @param x2 [read] The second data point.
#
def K(self, x1, x2):
assert len(x1.squeeze().shape) == 1, 'x1 is not a vector'
assert len(x2.squeeze().shape) == 1, 'x2 is not a vector'
return numpy.equal.outer(x1.squeeze(), x2.squeeze()).astype(numpy.int)
## Compute the kernel between the data points in x1 and those in x2.
# It returns a matrix with entry $(ij)$ equal to $K(x1_i, x1_j)$.
# If index1/index2 is
# specified, only those data points in x1/x2 with indices corresponding
# to index1/index2 are used to compute the kernel matrix. Furthermore,
# if output is specified, the provided buffer is used explicitly to
# store the kernel matrix.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Dot(self, x1, x2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert x1.shape[0] == x2.shape[0], \
'Argument 1 and Argument 2 have different dimensions'
if index1 is not None:
x1 = x1[index1,]
if index2 is not None:
x2 = x2[index2,]
# number of data points in x1.
n1 = x1.shape[0]
# number of data points in x2.
n2 = x2.shape[0]
# number of blocks.
nb = n1 / self._blocksize
if output is None:
output = numpy.zeros((n1,n2), numpy.int)
# handle special cases:
if index2 is not None:
if len(index2) <= self._blocksize:
output = numpy.equal.outer(x1.squeeze(), x2.squeeze()).astype(numpy.int)
return output
# blocking
lower_limit = 0
upper_limit = 0
for i in range(nb):
upper_limit = upper_limit + self._blocksize
output[lower_limit:upper_limit,] = numpy.equal.outer(x1[lower_limit:upper_limit,].squeeze(), x2.squeeze()).astype(numpy.int)
lower_limit = upper_limit
if lower_limit <= n1:
output[lower_limit:n1,] = numpy.equal.outer(x1[lower_limit:n1,].squeeze(), x2.squeeze()).astype(numpy.int)
return output
#------------------------------------------------------------------------------
## Joint kernel class for k(x,x')l(y,y')
#
# The methods are implemente efficiently by computing the resulting
# matrices block by block. Currently this class only works for the
# case where the inputs are integer.
#
class CJointKernel(CVectorKernel):
def __init__(self, y, xkernel=CLinearKernel(), ykernel=CDeltaKernel(), blocksize=128):
CVectorKernel.__init__(self, blocksize)
self._name = 'Joint kernel'
self.xkernel = xkernel
self.ykernel = ykernel
# remember the kernel matrix on y.
assert y is not None, 'y has be to provided for the initialization'
self.kyy = self.ykernel.Dot(y,y)
## Compute the kernel between the data points in x1 and those in x2,
# then multiply the resulting kernel matrix by alpha2.
# It returns a matrix with entry $(ij)$ equal to
# $sum_r K(x1_i,x2_r) \times alpha2_r$.
# Other parameters are defined similarly as those in Dot.
# @param x1 [read] The first set of data points.
# @param x2 [read] The second set of data points.
# @param alpha2 [read] The set of coefficients.
# @param index1 [read] The indices into the first set of data points.
# @param index2 [read] The indices into the second set of data points.
# @param output [write] The buffer where the output matrix is written into.
#
def Expand(self, x1, x2, alpha2, index1=None, index2=None, output=None):
assert len(x1.shape) == 2, 'Argument 1 has wrong shape'
assert len(x2.shape) == 2, 'Argument 2 has wrong shape'
assert len(alpha2.shape) == 2, 'Argument 3 has wrong shape'
assert x1.shape[1] == x2.shape[1], \
'Argument 1 and 2 has different dimesions'
assert x2.shape[0] == alpha2.shape[0], \
'Argument 2 and 3 has different number of data points'
L = self.kyy
if index1 is not None:
x1 = x1[index1,]
L = L[index1,]
if index2 is not None:
x2 = x2[index2,]
alpha2 = alpha2[index2,]
L = L[:,index2]
n1 = x1.shape[0]
nb = n1 / self._blocksize
n2 = alpha2.shape[1]
if output is None:
output = numpy.zeros((n1,n2), numpy.float64)
output = numpy.transpose(numpy.dot(L, self.xkernel.Expand(x1, x2, alpha2)))
return output
if __name__== '__main__':
n = 5000
x = numpy.random.rand(n,1)
y = (numpy.random.rand(n,1)*10).astype(numpy.int)
alpha = numpy.random.rand(n,4)
import time
t1 = time.clock()
kernel = CJointKernel(y)
K = kernel.Expand(x, x, alpha)
t2 = time.clock()
print K.shape, t2-t1