/**
* @file KCNullSpace.cpp
* @brief Classification and novelty detection with kernel null space methods (Kernel Null Foley Sammon Transform - KNFST)
* @author Paul Bodesheim
* @date 26/11/2012
*/
#include <iostream>
#include <sstream>
#include "core/vector/Algorithms.h"
#include "KCNullSpace.h"
#include <limits>
#undef DEBUG
using namespace NICE;
using namespace std;
using namespace OBJREC;
KCNullSpace::KCNullSpace( const Config *conf, Kernel *kernelFunction, const string & section )
: KernelClassifier ( conf, kernelFunction )
{
this->maxClassNo = 0;
this->verbose = conf->gB( section, "verbose", false );
}
KCNullSpace::KCNullSpace( const KCNullSpace &vcova ): KernelClassifier(vcova)
{
verbose = vcova.verbose;
dimNullSpace = vcova.dimNullSpace;
oneClassSetting = vcova.oneClassSetting;
trainingSetStatistic.clear();
for ( std::map<int,int>::const_iterator it = vcova.trainingSetStatistic.begin(); it != vcova.trainingSetStatistic.end(); it++ )
{
trainingSetStatistic.insert(pair<int,int>( (*it).first,(*it).second ));
}
nullProjectionDirections.resize( vcova.nullProjectionDirections.rows(),vcova.nullProjectionDirections.cols() );
nullProjectionDirections.set( 0.0 );
for(int i = 0; i < (int)vcova.nullProjectionDirections.rows(); i++)
{
for(int j = 0; j < (int)vcova.nullProjectionDirections.cols(); j++)
{
nullProjectionDirections(i,j) = vcova.nullProjectionDirections(i,j);
}
}
targetPoints.clear();
for( std::map<int,NICE::Vector>::const_iterator it = vcova.targetPoints.begin(); it != vcova.targetPoints.end(); it++ )
{
targetPoints.insert(pair<int,NICE::Vector>( (*it).first,(*it).second ));
}
eigenBasis.resize( vcova.eigenBasis.rows(),vcova.eigenBasis.cols() );
eigenBasis.set( 0.0 );
for(int i = 0; i < (int)vcova.eigenBasis.rows(); i++)
{
for(int j = 0; j < (int)vcova.eigenBasis.cols(); j++)
{
eigenBasis(i,j) = vcova.eigenBasis(i,j);
}
}
}
KCNullSpace::~KCNullSpace()
{
}
void KCNullSpace::teach ( KernelData *kernelData, const NICE::Vector & y )
{
NICE::Vector labels(y);
int maxLabel( (int)labels.Max());
/** check if we are in a one-class setting */
int minLabel( (int)labels.Min());
if (maxLabel == minLabel)
{
oneClassSetting = true;
maxClassNo = 1;
if (verbose)
std::cerr << "KCNullSpace::teach: one-class setting" << std::endl;
/** one-class setting: add a row and a column of zeros to the kernel matrix representing dot products with origin in kernel feature space*/
kernelData->increase_size_by_One();
kernelData->getKernelMatrix() (kernelData->getKernelMatrix().rows()-1, kernelData->getKernelMatrix().cols()-1) = 0.0;
labels.append(minLabel+1);
computeTrainingSetStatistic(labels);
}
else
{
oneClassSetting = false;
if (verbose)
std::cerr << "KCNullSpace::teach: multi-class setting" << std::endl;
computeTrainingSetStatistic(labels);
maxClassNo = trainingSetStatistic.size()-1; // number of classes, start counting at 0
}
computeNullProjectionDirections(kernelData,labels);
computeTargetPoints(kernelData,labels);
if (oneClassSetting)
{
/** remove the value that corresponds to the dot product with the origin in the kernel feature space, since this is always equal to zero */
nullProjectionDirections.deleteRow(nullProjectionDirections.rows()-1);
}
}
std::map<int,int> * KCNullSpace::getTrainingSetStatistic()
{
return &trainingSetStatistic;
}
NICE::Matrix KCNullSpace::getNullProjectionDirections()
{
return nullProjectionDirections;
}
std::map<int,NICE::Vector> * KCNullSpace::getTargetPoints()
{
return &targetPoints;
}
int KCNullSpace::getNullSpaceDimension()
{
return dimNullSpace;
}
bool KCNullSpace::isOneClass()
{
return oneClassSetting;
}
void KCNullSpace::computeTrainingSetStatistic(const NICE::Vector & y)
{
trainingSetStatistic.clear();
std::map<int,int>::iterator it;
for ( uint i = 0 ; i < y.size(); i++ )
{
it = trainingSetStatistic.find ( (int)y[i] );
if ( it == trainingSetStatistic.end() )
{
trainingSetStatistic.insert(it, pair<int,int>( (int)y[i],1 ));
}
else
{
it->second += 1;
}
}
}
void KCNullSpace::computeBasisUsingKernelPCA(const KernelData *kernelData)
{
if (verbose)
std::cerr << "KCNullSpace::computeBasisUsingKernelPCA: compute kernel PCA basis..." << std::endl;
NICE::Matrix K (kernelData->getKernelMatrix());
/** let K represent dot products of zero mean data in kernel feature space */
centerKernelMatrix(K);
/** get eigenvectors and eigenvalues (descreasing order) of centered kernel matrix*/
NICE::Matrix eigenVectors(K.rows(), K.cols(), 0.0);
NICE::Vector eigenValues(K.rows(), 0.0);
K.addIdentity(1.0);
eigenvectorvalues(K, eigenVectors, eigenValues);
eigenValues -= 1.0;
// NICE::GMCovariance gm(&K);
// NICE::EigValuesTRLAN ev;
// ev.getEigenvalues(gm, eigenValues, eigenVectors, K.rows());
/** only use eigenvectors of non-zero eigenvalues */
int j(0);
for (size_t i=0; i<K.rows(); i++)
{
if ( eigenValues(i) < 1e-12 )
{
eigenVectors.deleteCol(i);
}
else
{
j++;
}
}
eigenValues.resize(j);
/** scale eigenvectors with eigenvalues */
double scale(0.0);
for (size_t c=0; c<eigenVectors.cols(); c++)
{
scale = 1.0/sqrt(eigenValues(c));
for (size_t r=0; r<eigenVectors.rows(); r++)
eigenVectors(r,c) *= scale;
}
eigenBasis = eigenVectors;
if (verbose)
std::cerr << "KCNullSpace::computeBasisUsingKernelPCA: computation done" << std::endl;
}
void KCNullSpace::centerKernelMatrix(NICE::Matrix & kernelMatrix)
{
NICE::Matrix onesK (kernelMatrix.rows(), kernelMatrix.cols(), 1.0/kernelMatrix.rows());
kernelMatrix = kernelMatrix - onesK*kernelMatrix - kernelMatrix*onesK + onesK*kernelMatrix*onesK;
}
void KCNullSpace::computeNullProjectionDirections ( const KernelData *kernelData, const NICE::Vector & y )
{
if (verbose)
std::cerr << "KCNullSpace::computeNullProjectionDirections: compute null projection directions..." << std::endl;
/** obtain Kernel PCA basis */
computeBasisUsingKernelPCA(kernelData);
/** set matrix IM=(I-M) with I being the unit matrix and M being a matrix with all entries equal to 1/n where n is the number of training samples */
NICE::Matrix IM (y.size(),y.size(),-1.0/y.size());
IM.addIdentity(1.0);
/** compute matrix IL=(I-L) with I being the unit matrix and L being a block matrix where in each block of class samples each entry is equal to 1/numClassSamples */
NICE::Matrix IL (y.size(),y.size(),0.0);
IL.addIdentity(1.0);
for (size_t c=0; c<IL.cols(); c++)
{
/** if sample with index r is in the same class as sample with index c, then insert the value 1/numClassSamples */
for (size_t r=0; r<IL.rows(); r++)
{
if ( y(r) == y(c) )
{
IL(r,c) -= 1.0/trainingSetStatistic[(int)y(r)];;
}
}
}
/** compute Matrix H = ((I-M)*basisvecs)^T * kernelMatrix * (I-L) with I being the unit matrix */
NICE::Matrix H (y.size(),y.size(),0.0);
IM = IM*eigenBasis;
H = IM.transpose() * kernelData->getKernelMatrix() * IL;
/** obtain matrix T = H * H^T */
NICE::Matrix T = H*H.transpose();
/** get eigenvectors and eigenvalues (descreasing order) of T */
NICE::Matrix eigenVectors(T.rows(), T.cols(), 0.0);
NICE::Vector eigenValues(T.rows(), 0.0);
T.addIdentity(1.0);
eigenvectorvalues(T, eigenVectors, eigenValues);
eigenValues -= 1.0;
// NICE::GMCovariance gm(&T);
// NICE::EigValuesTRLAN ev;
// ev.getEigenvalues(gm, eigenValues, eigenVectors, T.rows());
if (verbose)
std::cerr << "T: " << T.rows() << " x " << T.rows() << " nan: " << T.containsNaN()<< std::endl;
/** only use eigenvectors of zero eigenvalues (null space!!!) but at least one eigenvector according to the smallest eigenvalue (therefore start at index i=T.rows()-2)*/
for (int i=T.rows()-2; i>=0; i--)
{
if ( eigenValues(i) > 1e-12 )
{
eigenVectors.deleteCol(i);
}
}
/** compute null projection directions */
nullProjectionDirections = IM*eigenVectors;
dimNullSpace = nullProjectionDirections.cols();
if (verbose)
std::cerr << "KCNullSpace::computeNullProjectionDirections: computation done" << std::endl;
}
void KCNullSpace::computeTargetPoints ( const KernelData *kernelData, const NICE::Vector & y )
{
if (verbose)
std::cerr << "KCNullSpace::computeTargetPoints: compute target points..." << std::endl;
targetPoints.clear();
NICE::Vector targetPoint (dimNullSpace, 0.0);
int classLabel(0);
if (oneClassSetting)
{
std::map<int,int>::iterator it;
it = trainingSetStatistic.begin();
classLabel = (*it).first;
for (size_t i=0; i<y.size(); i++)
{
if ( classLabel == y(i) )
{
/** project each training sample in the null space */
targetPoint += nullProjectionDirections.transpose() * kernelData->getKernelMatrix().getColumn(i);
}
}
/** average the projections of all class samples due to numerical noise */
targetPoint /= (*it).second;
/** we only have one target point in an one-class setting */
targetPoints.insert(pair<int,NICE::Vector>( classLabel,targetPoint));
}
else
{
/** create averaging vectors for each class, necessary to compute target points at the end of this method */
std::map<int,NICE::Vector> averagingVectors;
averagingVectors.clear();
NICE::Vector averagingVector(y.size(),0.0);
/** insert one averaging vector for each class */
int numClassSamples(0);
for ( std::map<int,int>::iterator it = trainingSetStatistic.begin(); it != trainingSetStatistic.end(); it++ )
{
/** create current averaging vector */
classLabel = (*it).first;
numClassSamples = (*it).second;
for ( size_t j = 0; j < y.size(); j++ )
{
if ( classLabel == y[j] )
{
averagingVector(j) = 1.0/numClassSamples;
}
else
{
averagingVector(j) = 0.0;
}
}
/** insert averaging vector for current class*/
averagingVectors.insert(pair<int,NICE::Vector>( classLabel,averagingVector));
}
/** compute target points using previously created averaging vectors: average for each class the projections of the class samples in the null space */
for ( std::map<int,NICE::Vector>::iterator it = averagingVectors.begin(); it != averagingVectors.end(); it++ )
{
targetPoint = nullProjectionDirections.transpose() * kernelData->getKernelMatrix() * (*it).second;
targetPoints.insert(pair<int,NICE::Vector>( (*it).first,targetPoint));
}
}
if (verbose)
std::cerr << "KCNullSpace::computeTargetPoints: computation done" << std::endl;
}
ClassificationResult KCNullSpace::classifyKernel ( const NICE::Vector & kernelVector, double kernelSelf ) const
{
if ( targetPoints.size() <= 0 )
fthrow(Exception, "The classifier was not trained with training data (use teach(...))");
if (oneClassSetting)
{
return noveltyDetection ( kernelVector, kernelSelf );
}
NICE::Vector projection(dimNullSpace,0.0);
projection = nullProjectionDirections.transpose() * kernelVector;
FullVector scores ( maxClassNo+1 );
scores.set(- std::numeric_limits<double>::max());
int iter(0);
for ( std::map<int,NICE::Vector>::const_iterator it = targetPoints.begin(); it != targetPoints.end(); it++ )
{
scores[iter] = 1.0-(it->second - projection).normL2();
iter++;
}
ClassificationResult r ( scores.maxElement(), scores );
return r;
}
ClassificationResult KCNullSpace::noveltyDetection ( const NICE::Vector & kernelVector, double kernelSelf ) const
{
if ( targetPoints.size() <= 0 )
fthrow(Exception, "The classifier was not trained with training data (use teach(...))");
NICE::Vector projection(dimNullSpace,0.0);
projection = nullProjectionDirections.transpose() * kernelVector;
FullVector scores ( 2 );
scores.set(- std::numeric_limits<double>::max());
double tmp_score(0.0);
for ( std::map<int,NICE::Vector>::const_iterator it = targetPoints.begin(); it != targetPoints.end(); it++ )
{
tmp_score = 1.0-(it->second - projection).normL2();
if (tmp_score > scores[1])
{
scores[1] = tmp_score;
}
}
ClassificationResult r ( scores.maxElement(), scores );
return r;
}
KCNullSpace* KCNullSpace::clone(void) const
{
KCNullSpace *classifier = new KCNullSpace( *this );
return classifier;
}
void KCNullSpace::clear()
{
//nothing to clear
}
void KCNullSpace::restore(std::istream& ifs, int type)
{
ifs.precision (numeric_limits<double>::digits10 + 1);
ifs >> maxClassNo;
ifs >> oneClassSetting;
ifs >> dimNullSpace;
ifs >> eigenBasis;
int k(0);
int classLabel(0);
int numClassSamples(0);
ifs >> k;
trainingSetStatistic.clear();
for (int i=0; i<k; i++)
{
ifs >> classLabel;
ifs >> numClassSamples;
trainingSetStatistic.insert( pair<int,int>(classLabel,numClassSamples) );
}
ifs >> nullProjectionDirections;
NICE::Vector targetPoint;
ifs >> k;
for (int i=0; i<k; i++)
{
ifs >> classLabel;
ifs >> targetPoint;
targetPoints.insert( pair<int,NICE::Vector>(classLabel,targetPoint) );
}
KernelClassifier::restore(ifs,type);
}
void KCNullSpace::store(std::ostream& ofs, int type) const
{
ofs.precision (numeric_limits<double>::digits10 + 1);
ofs << maxClassNo << endl;
ofs << oneClassSetting << endl;
ofs << dimNullSpace << endl;
ofs << eigenBasis << endl;
ofs << trainingSetStatistic.size() << endl;
for (std::map<int,int>::const_iterator it = trainingSetStatistic.begin() ; it != trainingSetStatistic.end(); it++)
{
ofs << (*it).first << endl;
ofs << (*it).second << endl;
}
ofs << nullProjectionDirections << endl;
ofs << targetPoints.size() << endl;
for (std::map<int,NICE::Vector>::const_iterator it = targetPoints.begin() ; it != targetPoints.end(); it++)
{
ofs << (*it).first << endl;
ofs << (*it).second << endl;
}
KernelClassifier::store(ofs,type);
}