multinormal.h

Go to the documentation of this file.

/*

Copyright (c) 2003, Cornell University
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

   - Redistributions of source code must retain the above copyright notice,
       this list of conditions and the following disclaimer.
   - Redistributions in binary form must reproduce the above copyright
       notice, this list of conditions and the following disclaimer in the
       documentation and/or other materials provided with the distribution.
   - Neither the name of Cornell University nor the names of its
       contributors may be used to endorse or promote products derived from
       this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.

*/

// -*- C++ -*-

#if !defined _CLUS_MULTINORMAL_H_
#define _CLUS_MULTINORMAL_H_

#include "general.h"
#include "distribution.h"
#include "vec.h"
#include "cholesky.h"
#include <math.h>
#include "linearregressor.h"
#include "misc.h"
#include "trislv.h"

using namespace TNT;

namespace CLUS
{

/** Implements a multidimentional normal distribution. */
class MultiDimNormal : public Distribution
{
protected:
    /// the expected value
    Vector<double> mu;
    
    /// the equation of the hyperplane
    Vector<double> C;
    
    /// the covariance matrix's cholesky, upper diagonal
    /// lower diagonal
    Fortran_Matrix<double> sigmaChol;

    /// the determinant of the cholesky factor of Sigma
    double alpha_over_detChol_inDimp1;    
    double alpha_over_detChol_inDim;

    /// number of datapoints
    int N;

    /// sum probabilities
    double S;

    /// sum probab*x
    Vector<double> S_x;

    /// sum prob*xx^T, upper symm diag
    Fortran_Matrix<double> S_xxT;

    /// the radius of the distribution. Used to generate distributions close by
    double radius; 

    /** Computes alpha_over_detChol_inDimp1 and alpha_over_detChol_inDim given sigmaChol */
    void ComputeAlphas(void)
    {
        double detChol=1.0;
        for (int i=1; i<=inDim; i++)
            detChol*=sigmaChol(i,i);
        alpha_over_detChol_inDim=weight*pow(2*M_PI,-inDim/2.0)/detChol;
        detChol*=sigmaChol(inDim+1,inDim+1);
        alpha_over_detChol_inDimp1=weight*pow(2*M_PI,-(inDim+1)/2.0)/detChol;
        // print debug info
        if (alpha_over_detChol_inDimp1<0.0)
        {
            cerr << "Problem with alphas" << endl;
            cerr << "mu: " << mu << " sigmaChol: " << sigmaChol << " ao=" << alpha_over_detChol_inDim;
            cerr << " aop1=" << alpha_over_detChol_inDimp1 << endl;
            cerr << "S_xxT=" << S_xxT;
        }
    }

    /** Computes the equation of the hyperplane using Sigma=S_xxT.
    Must be called when S_xxT contains Sigma not the partial sums only.
    
    @return true if anything is fine, false if something is wrong and the computation is bogus.*/
    bool ComputeCCoef(Fortran_Matrix<double>& Sigma, bool isCholesky)
    {
        if (weight==0.0)
            return false;

#if 0

        Fortran_Matrix<double> S(inDim+1,inDim+1);
        if (isCholesky)
        {
            for (int i=1; i<=inDim+1; i++)
            {
                for (int j=i; j<=inDim+1; j++)
                {
                    // compute G(i,:)*G'(j,:)
                    double tmp=0.0;
                    for (int k=1; k<=i; k++)
                        tmp+=Sigma(i,k)*Sigma(j,k);
                    S(i,j)=tmp;
                }
            }
        }
        else
        {
            S=Sigma;
        }

        Vector<double> u(inDim+1);

        // compute u
        double lambda,lambda2;
        double slimness;
        if (MinEigenvalueUpperSymmetric(S,lambda,lambda2,slimness))
        {
            if (!EigenVector(S,lambda,u,1))
            {
                cerr << "OOPS could not find the eigenvector for " << lambda << endl;
                return false;
            }
        }
        else
        {
            cerr << "OOPS-eigenvalue problem not solved" << endl;
            return false;
        }
        if (lambda < -SMALL_POZ_VALUE)
        {
            cerr << "Matrix not positive definite since the smallest eigv=" << lambda << endl;
            return false;
        }

        // slimness is the average eigvalue nonminimal. The square of this is the size of the cluster
        radius = sqrt(slimness);

        // compute C from u
        double uy=u[inDim];

        if ( fabs(uy) < 1.0/(inDim*1000) )
        {
            cerr << "Smallest eigenvector oriented badly. Something is wrong, killing distrib" << endl;
            return false;
        }

        C[0]=mu[inDim];
        for(int j=0;j<inDim;j++)
        {
            C[j+1]=-u[j]/uy;
            C[0]+=mu[j]*u[j]/uy;
        }
        return true;

#else

        ExitIf(!isCholesky,"Not implemented");

        // use least square
        // form b
        Vector<double> b(inDim+1);
        // compute in radius the trace of original S
        radius=0.0;
        for (int i=1; i<=inDim; i++)
        {
            b(i)=Sigma(inDim+1,i);
            radius+=pow2(Sigma(i,i));
        }
        b(inDim+1)=0.0;
        radius+=pow2(Sigma(inDim+1,inDim+1));
        radius=sqrt(radius/inDim); // suppose one of the eigenvalues is small.

        // solve for c=G\b
        Vector<double> c=Transposed_Upper_triangular_solve< Fortran_Matrix<double>,Vector<double> >(Sigma,b);

        // form C
        C[0]=mu[inDim];
        for(int i=1; i<=inDim; i++)
        {
            C[0]-=mu[i-1]*c[i-1];
            C[i]=c[i-1];
        }
        return true;

#endif

    }

public:
    MultiDimNormal(int InDim=0):Distribution(InDim),mu(InDim+1),C(InDim+1),
            sigmaChol(InDim+1,InDim+1),S_x(InDim+1),S_xxT(InDim+1,InDim+1)
    {
        alpha_over_detChol_inDimp1=1.0;
        alpha_over_detChol_inDim=1.0;
        N=0;
        S=0;
        S_x=0;
        S_xxT=0;
        weight=1.0;
    }

    MultiDimNormal(int InDim,int):Distribution(InDim),mu(InDim+1),C(InDim+1),
            sigmaChol(InDim+1,InDim+1),S_x(0),S_xxT(0,0)
    {
        alpha_over_detChol_inDimp1=1.0;
        alpha_over_detChol_inDim=1.0;
        N=0;
        S=0;
        S_x=0;
        S_xxT=0;
        weight=1.0;
    }

    MultiDimNormal(int InDim, double Alpha, Vector<double>& Mu, Fortran_Matrix<double>& SigmaChol):
            Distribution(InDim), mu(Mu), C(0), sigmaChol(SigmaChol), S_x(0),S_xxT(0,0)
    {
        N=0;
        S=0;
        S_x=0;
        S_xxT=0;
        weight=Alpha;
        ComputeAlphas();
    }

    Regressor* CreateRegressor(void)
    {
        if ( ComputeCCoef(sigmaChol,true) )
            return new LinearRegressor(C);
        else
            return new LinearRegressor();
    }

    /** Reverses the normalization transformation in the original space
    mu=pD.sigmaChol*mu+pD.mu;
    sigmaChol=pD.sigmaChol*sigmaChol;
    */
    void DenormalizeParameters(MultiDimNormal* pD)
    {
        // compute mu=pD.sigmaChol*mu+pD.mu
        // do it in reverse order to avoid keeping a copy of mu
        for(int i=inDim+1; i>=1; i--)
        {
            double tmp=0.0;
            for (int j=1; j<=i; j++)
                tmp+=pD->sigmaChol(i,j)*mu[j-1];
            mu[i-1]=tmp+pD->mu[i-1];
        }

        // compute sigmaChol=pD.sigmaChol*sigmaChol
        // do it in reverse order on i to avoid keeping a copy of sigmaChol
        for(int i=inDim+1; i>=1; i--)
            for (int j=1; j<=i; j++)
            {
                double tmp=0.0;
                for (int k=j; k<=i; k++)
                    tmp+=pD->sigmaChol(i,k)*sigmaChol(k,j);
                sigmaChol(i,j)=tmp;
            }
    }

    /** Transform data such that it looks like a 0 centered sphere
    of radius 1
    */
    void NormalizeData(const double* DataCache, double* X)
    {
        // X=(DataCache-mu)*sigmaChol^-1
        for (int i=0; i<=inDim; i++)
            X[i]=DataCache[i]-mu[i];

        // solve X=GY for Y into X
        for(int i=1; i<=inDim+1; i++)
        {
            double tmp=0.0;
            for (int j=1; j<i; j++)
                tmp+=sigmaChol(i,j)*X[j-1];

            X[i-1]=(X[i-1]-tmp)/sigmaChol(i,i);
        }
    }

    double LearnProbability(const double* DataCache)
    {
        dataCache=DataCache;
        double X[MAX_IN_DIM];
        // double* Xm1=X-1; // pointer to the element before the start of X such that Xm1[i]=X[i+1]

        if (weight==0.0)
        {
            probabilityLearn=0.0;
            return 0.0;
        }

        for (int i=0; i<=inDim; i++)
            X[i]=dataCache[i]-mu[i];
        double norm2=0.0;
        // solve X=GY for Y into X and compute X^T*X
        for(int i=1; i<=inDim+1; i++)
        {
            double tmp=0.0;
            for (int j=1; j<i; j++)
                tmp+=sigmaChol(i,j)*X[j-1];

            double sCii=sigmaChol(i,i);
            if (sCii == 0.0)
                if (X[i-1]!=tmp)
                    return (probabilityLearn=0.0); // since norm2=inf => e^{-norm2/2)=0.0
                else // 0.0/0.0 situation, point right on the hyperplane, probability=HUGE
                    return (probabilityLearn=MAXREAL);

            X[i-1]=(X[i-1]-tmp)/sCii;
            norm2+=pow2(X[i-1]);
        }

        probabilityLearn=alpha_over_detChol_inDimp1*exp(-norm2/2.0);

        // assert(probabilityLearn>=0.0 && finite(probabilityLearn) );

        return probabilityLearn;
    }

    double InferProbability(const double* DataCache)
    {
        dataCache=DataCache;
        double X[MAX_IN_DIM];
        // double* Xm1=X-1;

        if (weight==0.0)
            return 0.0;

        for (int i=0; i<inDim; i++)
            X[i]=dataCache[i]-mu[i];
        double norm2=0.0;
        // solve X=GY for Y into X and compute X^T*X, ignore last line in G
        for(int i=1; i<=inDim; i++)
        {
            double tmp=0.0;
            for (int j=1; j<i; j++)
                tmp+=sigmaChol(i,j)*X[j-1];
            X[i-1]=(X[i-1]-tmp)/sigmaChol(i,i);
            norm2+=pow2(X[i-1]);
        }

        probabilityInfer=alpha_over_detChol_inDim*exp(-norm2/2.0);
        return probabilityInfer;
    }

    void UpdateStatistics(const double* DataCache, double prob)
    {
        N++;
        S+=prob;
        for (int i=0; i<=inDim; i++)
        {
            double aux=prob*DataCache[i];
            S_x[i]+=aux;
            // fill just the upperpart of S_xxT
            for (int j=i; j<=inDim; j++)
                S_xxT(i+1,j+1)+=aux*DataCache[j];
        }
    }

    double NormalizeLearnProbability(double Coef, int nrClus=1)
    {
        double pLearn=Distribution::NormalizeLearnProbability(Coef,nrClus);

        assert(pLearn>=0.0 && pLearn <=1.0 && finite(pLearn) );

        UpdateStatistics(dataCache, pLearn);
        return 0.0;
    }

    double UpdateParameters(void)
    {
        double oldMu, distP=0.0;

        // if the cluster is dead do nothing
        if (weight==0.0)
            return 0.0;

        if (S<=(inDim+1))
        {
            cerr << "Cluster got too small. We have to throw it out since is not anymore reliable." << endl;
            weight=0.0;
            distP=1.0;
            goto cleanup;
        }

        // don't change the weight to have clusters of about the same size
        // weight=S/N; // S is the number of points on which this clusterdepends


        for (int i=1; i<=inDim+1; i++)
        {
            oldMu=mu(i);
            mu(i)=S_x(i)/S;
            for (int j=i; j <=inDim+1; j++)
            {
                S_xxT(i,j)=(S_xxT(i,j)-mu(i)*S_x(j))/S;

                // fix in S_xxT the almost 0.0 entries. This helpes a lot
                if (fabs(S_xxT(i,j))<SMALL_POZ_VALUE)
                {
                    S_xxT(i,j)=0.0;
                }
            }
            distP+=pow2(oldMu-mu(i));
        }


        // compute Cholesky factorization of S_xxT into sigmaChol
        int cholRet;
        if ( (cholRet=Cholesky_upper_factorization(S_xxT,sigmaChol)) )
        {
            // If we got in here we have trouble
            // If the cluster does not depend on at least 2*(inDim+1) points we just throw it out
            if (S<2*(inDim+1))
            {
                cerr << "Troublesome cluster too small to bother. We throw it out" << endl;
                weight=0.0;
                distP=1.0;
                goto cleanup;
            }
            else
            {
                cerr << "Matrix S_xxT is not positive definite: " << S_xxT ;
                cerr << "sigmaChol: " << sigmaChol;
                cerr << "weight=" << weight << endl;
                // fix the last value of sigmaChol
                if (cholRet==2)
                {
                    cout << "Fixing sigmaChol" << endl;
                    sigmaChol(inDim+1,inDim+1)=0.0; // only the last value is bad, we can fix it
                    distP=-1.0; // this will make the returned value nan and terminate anything
                }
                else
                {
                    // probably the cluster is a "squeezy"(has one less dimention in the input space).
                    cerr << "The cluster is a squeezy, setting weight to 0." << endl;
                    weight=0.0;
                    distP=1.0;
                    goto cleanup;
                }
            }
        }

        ComputeAlphas();

        /* // need this to have all the time a regressor
        if (!ComputeCCoef(S_xxT))
        weight=0.0;

        */
cleanup:
        // erase N, S, S_x, S_xxT
        N=0;
        S=0.0;
        S_x=0; // the operator= takes care of this
        S_xxT=0;

        return sqrt(distP)/(inDim+1); //how much the center has moved
    }

    void AdjustYwithYoverD(double& y)
    {
        double dy=C[0];
        for(int j=0; j<inDim; j++)
            dy+=C[j+1]*dataCache[j];
        dy*=probabilityInfer;
        y+=dy;
    }

    static string TypeName(void)
    {
        return string("MultiDimNormal");
    }

    void RandomDistribution(int NrClusters)
    {
        const double smallness=0.1;
        const double sizeDist=1/pow(NrClusters,1.0/(inDim+1));
        int d=inDim+1;

        // reset weithts to 1.
        weight = 1.0/NrClusters;
        // choose mu random
        for(int i=0; i<d ;i++)
        {
            mu[i]=RANDOM01FLOAT*2-1.0;
        }

        do
        {
            // choose n randomly and build V, set of vectors orthogonal on n
            Vector<double> n(d);
            for (int i=0; i<d; i++)
                n[i]=RANDOM01FLOAT*2-1.0;
            // normalize n
            double norm_n=dot_prod(n,n);
            for (int i=0; i<d; i++)
                n[i]/=norm_n;

            // set Sigma as I-0.1n*n'
            // in this way Sigma has n as the eigenvector coresp to eightvalue 0.1
            // and d-1 other eigenvectors orth. on n with eigvals 1.0
            for (int i=1; i<=d; i++)
            {
                sigmaChol(i,i)=sizeDist*(1-smallness*pow2(n(i)));
                for (int j=i+1; j<=d; j++)
                    sigmaChol(i,j)=-sizeDist*(1-smallness)*n(i)*n(j);
            }

        }
        while (Cholesky_upper_factorization(sigmaChol,sigmaChol)!=0);

        ComputeAlphas();
    }

    /** For hierarchies or trees we need to produce offsprings close to the parent*/
    void RandomDistribution(int NrClusters, MultiDimNormal& parent)
    {
        // reset weitht to 1.0. It might have been 0.
        weight = 1.0;

        // Create a distribution 1/NrClusters as small as the parent with a center
        // around the center of the parent
        radius=parent.radius/pow(NrClusters,1.0/inDim);
        mu[inDim]=parent.C[0];
        for(int i=0;i<inDim;i++)
        {
            mu[i]=parent.mu[i]+(2*parent.radius-radius)*RANDOM01FLOAT;
            mu[inDim]+=mu[i]*parent.C[i+1];
        }

        // comment this line if the parent.C is used
        // mu[inDim]=parent.mu[inDim];

        // put sigma=I => sigmaChol=I*radius
        sigmaChol=0.0;
        for( int i=1; i<=inDim+1; i++)
            sigmaChol(i,i)=radius;

        ComputeAlphas();
    }

    void SaveToStream(ostream& out)
    {
        out << " [ ";
        // write mu
        for (int i=0; i<inDim+1; i++)
            out << mu[i] << " ";
        out << " ] ( ";
        // write sigmaChol
        for (int i=1; i<=inDim+1; i++)
            for (int j=1; j<=i; j++)
                out << sigmaChol(i,j) << " ";
        out << ") { " << weight << " } " << endl;
        /* for (int i=0; i<inDim+1; i++)
        out << C[i] << " ";
        out << endl;*/
    }

    void LoadFromStream(istream& in)
    {
        string s;

        in >> s;
        if (s!="[")
            throw ErrMsg("[ missing in MultiDimNormal");
        for (int i=0; i<inDim+1; i++)
            in >> mu[i];
        in >> s;
        if (s!="]")
            throw ErrMsg("] missing in MultiDimNormal");
        in >> s;
        if (s!="(")
            throw ErrMsg("( missing in MultiDimNormal");
        for (int i=1; i<=inDim+1; i++)
            for (int j=1; j<=i; j++)
                in >>  sigmaChol(i,j);
        in >> s;
        if (s!=")")
            throw ErrMsg(") missing MultiDimNormal");
        in >> s;
        if (s!="{")
            throw ErrMsg("{ missing in MultiDimNormal");
        in >> weight;
        in >> s;
        if (s!="}")
            throw ErrMsg("} missing in MultiDimNormal");

        ComputeCCoef(sigmaChol,true);
        ComputeAlphas();

    }
};

}
#endif // _CLUS_MULTINORMAL_H_

---