matrix.h

这是一款用C++编写的实现gmm算法的程序

        _[j*column+i] -= scalar;
    return *this;
  }

  inline Matrix& operator *= (float scalar)
  {
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] *= scalar;
    return *this;
  }

  inline Matrix& operator /= (float scalar)
  {
    scalar = 1.0f/scalar;
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] *= scalar;
    return *this;
  }
  
  inline Matrix operator + (const Matrix &matrix) const
  {
    Matrix result(row,column,false);  
    return Add(matrix,result);
  }
  
  inline Matrix& Add(const Matrix &matrix, Matrix &result) const
  {   
    result.Resize(row,column,false);
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++){
      for (unsigned int i = 0; i < ki; i++)
        result._[j*column+i] = _[j*column+i] + matrix._[j*column+i];
      for (unsigned int i = ki; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    }
    for (unsigned int j = kj; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    return result;
  }
  
  inline Matrix operator - (const Matrix &matrix) const
  {
    Matrix result(row,column,false);  
    return Sub(matrix,result);
  }
  
  inline Matrix& Sub(const Matrix &matrix, Matrix &result) const
  {   
    result.Resize(row,column,false);
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++){
      for (unsigned int i = 0; i < ki; i++)
        result._[j*column+i] = _[j*column+i] - matrix._[j*column+i];
      for (unsigned int i = ki; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    }
    for (unsigned int j = kj; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    return result;
  }
  
  inline Matrix operator ^ (const Matrix &matrix) const
  {
    Matrix result(row,column,false);  
    return PMult(matrix,result);
  }
  
  inline Matrix& PMult(const Matrix &matrix, Matrix &result) const
  {   
    result.Resize(row,column,false);
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++){
      for (unsigned int i = 0; i < ki; i++)
        result._[j*column+i] = _[j*column+i] * matrix._[j*column+i];
      for (unsigned int i = ki; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    }
    for (unsigned int j = kj; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    return result;
  }
  
  inline Matrix operator / (const Matrix &matrix) const
  {
    Matrix result(row,column,false);  
    return PDiv(matrix,result);
  }
  
  inline Matrix& PDiv(const Matrix &matrix, Matrix &result) const
  {   
    result.Resize(row,column,false);
    const unsigned int kj = (row<=matrix.row?row:matrix.row);
    const unsigned int ki = (column<=matrix.column?column:matrix.column);
    for (unsigned int j = 0; j < kj; j++){
      for (unsigned int i = 0; i < ki; i++)
        result._[j*column+i] = _[j*column+i] / matrix._[j*column+i];
      for (unsigned int i = ki; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    }
    for (unsigned int j = kj; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i];  
    return result;
  }

  inline Matrix operator + (float scalar) const
  {
    Matrix result(row,column,false);  
    return Add(scalar,result);    
  }

  inline Matrix& Add(float scalar, Matrix& result) const
  {
    result.Resize(row,column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i] + scalar;    
    return result;
  }

  inline Matrix operator - (float scalar) const
  {
    Matrix result(row,column,false);  
    return Sub(scalar,result);    
  }
  
  inline Matrix& Sub(float scalar, Matrix& result) const
  {
    result.Resize(row,column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i] - scalar;    
    return result;
  }

  inline Matrix operator * (float scalar) const
  {
    Matrix result(row,column,false);  
    return Mult(scalar,result);    
  }

  inline Matrix& Mult(float scalar, Matrix& result) const
  {
    result.Resize(row,column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i] * scalar;    
    return result;
  }


  inline Matrix operator / (float scalar) const
  {
    Matrix result(row,column,false);  
    return Div(scalar,result);    
  }

  inline Matrix& Div(float scalar, Matrix& result) const
  {
    result.Resize(row,column,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*column+i] = _[j*column+i] / scalar;    
    return result;    
  }

  inline bool operator == (const Matrix& matrix) const
  {
    if((row!=matrix.row)||(column!=matrix.column)) return false;
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        if(_[j*column+i] != matrix._[j*column+i]) return false;
    return true;
  }

  inline bool operator != (const Matrix& matrix) const
  {
    return !(*this ==  matrix);
  }

  inline Vector operator * (const Vector &vector) const
  {
    Vector result(row,false);  
    return Mult(vector,result);    
  }  

  inline Vector Mult(const Vector &vector) const
  {
    Vector result(row,false);  
    return Mult(vector,result);    
  }

  inline Vector& Mult(const Vector &vector, Vector &result) const
  {
    result.Resize(row,false);
    const unsigned int ki = (column<=vector.row?column:vector.row);
    for (unsigned int j = 0; j < row; j++){
      result._[j] = 0.0f;
      for (unsigned int i = 0; i < ki; i++)
        result._[j] += _[j*column+i] * vector._[i];
    }
    return result;
  }


  inline Matrix operator * (const Matrix &matrix) const  
  {
    Matrix result(row,matrix.column,false);  
    return Mult(matrix,result);
  }  

  inline Matrix& Mult(const Matrix &matrix, Matrix &result) const
  {
    result.Resize(row,matrix.column,false);
    const unsigned int rrow = result.row;
    const unsigned int rcol = result.column;
    const unsigned int kk = (column<=matrix.row?column:matrix.row);
    for (unsigned int j = 0; j < rrow; j++){
      for (unsigned int i = 0; i < rcol; i++){
        result._[j*rcol+i] = 0.0f;
        for(unsigned int k = 0; k< kk; k++)    
          result._[j*rcol+i] += _[j*column+k] * matrix._[k*rcol+i];
      }
    }    
    return result;
  }



  inline Matrix& Identity()
  {
    const unsigned int k = (row>column?column:row);
    Zero();
    for (unsigned int i = 0; i < k; i++)
      _[i*column+i] = 1.0f;
    return *this;    
  }

  inline Matrix& Diag(const Vector &vector)
  {
    const unsigned int k = (row>column?column:row);
    const unsigned int k2 = (k>vector.row?vector.row:k);
    Zero();
    for (unsigned int i = 0; i < k2; i++)
      _[i*column+i] = vector._[i];
    return *this;    
  }

  inline Matrix& Random(){
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        _[j*column+i] =((float)rand())/((float)(RAND_MAX+1.0));    
    return *this;    
  }

	inline Matrix Transpose() const
	{
    Matrix result(row,column,false);
    return Transpose(result);    
	}

  inline Matrix& Transpose(Matrix &result) const
  {    
    result.Resize(column,row,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[i*row+j] = _[j*column+i];
    return result;    
  }

  inline Matrix VCat(const Matrix& matrix)
  {
    Matrix result;
    return VCat(matrix,result);    
  }
  
  inline Matrix& VCat(const Matrix& matrix, Matrix & result)
  {
    unsigned int k1 = (column>matrix.column?column:matrix.column);
    result.Resize(row+matrix.row,k1,false);
    for (unsigned int j = 0; j < row; j++){
      for (unsigned int i = 0; i < column; i++)
        result._[j*k1+i] = _[j*column+i];
      for (unsigned int i = column; i < k1; i++)
        result._[j*k1+i] = 0.0f;
    }
    for (unsigned int j = 0; j < matrix.row; j++){
      for (unsigned int i = 0; i < matrix.column; i++)
        result._[(row+j)*k1+i] = matrix._[j*matrix.column+i];
      for (unsigned int i = matrix.column; i < k1; i++)
        result._[(row+j)*k1+i] = 0.0f;
    }
    return result;
  }

  inline Matrix HCat(const Matrix& matrix)
  {
    Matrix result;
    return HCat(matrix,result);    
  }
  
  inline Matrix& HCat(const Matrix& matrix, Matrix & result)
  {
    unsigned int k1 = (row>matrix.row?row:matrix.row);
    unsigned int k2 = column+matrix.column;
    result.Resize(k1,k2,false);
    for (unsigned int j = 0; j < row; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*k2+i] = _[j*column+i];
    for (unsigned int j = row; j < k1; j++)
      for (unsigned int i = 0; i < column; i++)
        result._[j*k2+i] = 0.0f;

    for (unsigned int j = 0; j < matrix.row; j++)
      for (unsigned int i = 0; i < matrix.column; i++)
        result._[j*k2+i+column] = matrix._[j*matrix.column+i];
    for (unsigned int j = matrix.row; j < k1; j++)
      for (unsigned int i = 0; i < matrix.column; i++)
        result._[j*k2+i+column] = 0.0f;
    
    return result;
  }

  inline static int IsInverseOk(){
    return bInverseOk;
  }

  inline Matrix Inverse(float *determinant=NULL) const
  {
    Matrix result;
    return Inverse(result,determinant);
  }  

  inline Matrix& Inverse(Matrix &result, float *determinant=NULL) const
  {
    bInverseOk = true;       
    if(row==column){ // Square matrix
      if(determinant!=NULL) *determinant = 1.0f;
      result.Resize(row,column,false);
      const unsigned int n = row;
      Matrix MM(*this);
      result.Identity();
      for(unsigned int i=0;i<n;i++){
        float pivot = MM._[i*column+i]; 
        if(fabs(pivot)<=EPSILON){
          for(unsigned int j=i+1;j<n;j++){
            if((pivot = MM._[j*column+i])!=0.0f){
              MM.SwapRow(i,j);
              result.SwapRow(i,j);
              break;  
            }
          }            
          if(fabs(pivot)<=EPSILON){
            bInverseOk = false;
            if(determinant!=NULL) *determinant = 0.0f;
            return result;
          }                      
        }
        if(determinant!=NULL) *determinant *= pivot;
        pivot = 1.0f/pivot;
        for(unsigned int j=0;j<n;j++){
          MM._[i*column+j]   *= pivot;
          result._[i*column+j] *= pivot;
        }
        for(unsigned int k=0;k<n;k++){
          if(k!=i){
            const float mki = MM._[k*column+i];
            for(unsigned int j=0;j<n;j++){
               MM._[k*column+j]   -= MM._[i*column+j]   *mki;
               result._[k*column+j] -= result._[i*column+j] *mki;              
            }            
          }
        }
      }        
    }else{ // Moore-Penrose pseudo inverse
      if(determinant!=NULL) *determinant = 0.0f;
      if(row>column){ // (JtJ)^(-1)Jt
        Matrix MT,SQ,SQInv;
        Transpose(MT);
        MT.Mult(*this,SQ);
        SQ.Inverse(SQInv);
        SQInv.Mult(MT,result); 
      }else{ // Jt(JJt)^(-1)
        Matrix MT,SQ,SQInv;
        Transpose(MT);
        Mult(MT,SQ);
        SQ.Inverse(SQInv);
        MT.Mult(SQInv,result);         
      }
    }
    return result;    
  }

  inline Matrix& SwapRow(unsigned int j1, unsigned int j2){
    if((j1<row)&&(j2<row)){
      float tmp;
      for (unsigned int i = 0; i < column; i++){
        tmp            = _[j1*column+i];
        _[j1*column+i] = _[j2*column+i];
        _[j2*column+i] = tmp;        
      }        
    }
    return *this; 
  }