Gaussian Multivariate Distribution -Part 1

Introduction

In this article, we will look at the multivariate Gaussian distribution.

• The mutivariate Gaussian distribution is parameterized by mean vector µ and covariance matrix Σ
  f x ( x 1 , … , x k ) = 1 √ ( 2 π ) k | Σ | exp ( − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) ) ,
  w h e r e
  μ is a Nx1 mean vector
  Σ is NxN matrix covariance matrix
  | Σ | is the determinant of Σ
• The mutivariate Gaussian is also used as probability density function of the vector <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow class="MJX-TeXAtom-ORD"> <mi>X
• A naive way of implementation is to directly perform the computation as in the equation to get the result.
• A Gaussian is defined
• Since ill be using a lot of opencv conversion routines are defined for Opencv cv::Mat data structure to Eigen::MatrixXf data structure.
• When loading the covarince matrix ,the determinant, inverse, etc. are also computed which would be later required for probability calculation:
  C#

  void setMean(Mat &v)
  {
      Map<Eigen::Matrix<float,Dynamic,Dynamic,Eigen::RowMajor> > 
      mappedMat((float *)v.data(),1,v.size()) ;
      _mu=mappedMat;
  }
  
  void setSigma(cv::Mat &v)
  {
      Map<Eigen::Matrix<float,Dynamic,Dynamic,Eigen::RowMajor> > 
      mappedMat((float *)v.data,v.rows,v.cols) ;
      _Sigma=mappedMat;
      _dim=v.rows;
      _det=_Sigma.determinant();
      _scale=1.f/(pow(2*PI*_dim*_det,0.5));
      _invsigma=_Sigma.inverse();
  }
  
  MatrixXf setData(cv::Mat &v)
  {
      Map<Eigen::Matrix<float,Dynamic,Dynamic,Eigen::RowMajor> > 
      mappedMat((float *)v.data,1,v.cols) ;
      MatrixXf  ref=mappedMat;
      return ref;
  }

• The probability computation is simply writing out the formula:
  C#

     void validate(float &res)
     {
  if(res<0)
             res=0;
         if(res>1)
             res=1;
         if(isnan(res))
             res=0;
      }
      
     float Prob(Mat &x)
     {
         MatrixXf tmp=setData(x);
         float res=0;
         MatrixXf tmp1=(tmp-_mu);
         MatrixXf tmp2=tmp1*_invsigma*tmp1.transpose();
         res=scale*tmp2(0,0);        
  validate(res);   
         return res;
     }

• The covariance matrix is positive semidefinate and symmetric
• The Cholesky decomposition of a symmetric, positive definite matrix <math xmlns="http://www.w3.org/1998/Math/MathML"> <mi mathvariant="bold">Σ decomposes A into a product of a lower triangular matrix L and its transpose.
• Thus the product <math display="block" xmlns="http://www.w3.org/1998/Math/MathML"> <mtable columnalign="right center left" columnspacing="0 thickmathspace" displaystyle="true" rowspacing="3pt"> <mtr> <mtd> <mtd> <mi mathvariant="bold">Σ=LL^T
  Σ⁻¹=(LL^T)⁻¹=L^−TL⁻¹
  y^TΣ⁻¹y
  y^TΣ⁻¹y
  y^TL^−TL⁻¹y
  (L⁻¹y)^T(L⁻¹y)<mtr><mtd><mrow><mo>
• Thus, we can only compute <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow> <mo>(L⁻¹y)<mo> and take its transpose
• Thus, we can perform the Cholskey factorization of the covariance matrix to find <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow class="MJX-TeXAtom-ORD"> <mi mathvariant="bold">L
• Then, we take the inverse of <math xmlns="http://www.w3.org/1998/Math/MathML"> <mrow class="MJX-TeXAtom-ORD"> <mi mathvariant="bold">L:
  C#

  float Prob(Mat &x)
  {
      MatrixXf tmp=setData(x);
  
      float res=0;
      MatrixXf tmp1=(tmp-_mu).transpose();
      tmp1=_LI*tmp1;
      MatrixXf tmp2=tmp1.transpose()*tmp1;
      res=tmp2(0,0);
      validate(res);
      return res;
  }
  

Source Code

The code can be found at git repository https://github.com/pi19404/OpenVisionin files ImgML/gaussian.hpp and ImgML.gaussian.cpp files.

The PDF Version of the document can be found below:

• Gaussian Multivariate Distribution -Part 1 by pi194043