NormalizeComponentsWhitening.h
Go to the documentation of this file.
1 //===========================================================================
2 /*!
3  *
4  *
5  * \brief Data normalization to zero mean, unit variance and zero covariance
6  *
7  *
8  *
9  *
10  * \author T. Glasmachers
11  * \date 2010
12  *
13  *
14  * \par Copyright 1995-2015 Shark Development Team
15  *
16  * <BR><HR>
17  * This file is part of Shark.
18  * <http://image.diku.dk/shark/>
19  *
20  * Shark is free software: you can redistribute it and/or modify
21  * it under the terms of the GNU Lesser General Public License as published
22  * by the Free Software Foundation, either version 3 of the License, or
23  * (at your option) any later version.
24  *
25  * Shark is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
28  * GNU Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public License
31  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
32  *
33  */
34 //===========================================================================
35 
36 
37 #ifndef SHARK_ALGORITHMS_TRAINERS_NORMALIZECOMPONENTSWHITENING_H
38 #define SHARK_ALGORITHMS_TRAINERS_NORMALIZECOMPONENTSWHITENING_H
39 
40 
43 #include <shark/Data/Statistics.h>
45 
46 namespace shark {
47 
48 
49 /// \brief Train a linear model to whiten the data
50 class NormalizeComponentsWhitening : public AbstractUnsupervisedTrainer<LinearModel<RealVector> >
51 {
53 public:
54 
56  NormalizeComponentsWhitening(double targetVariance = 1.0){
57  SHARK_CHECK(targetVariance > 0.0, "[NormalizeComponentsWhitening::NormalizeComponentsWhitening] target variance must be positive");
58  m_targetVariance = targetVariance;
59  }
60 
61  /// \brief From INameable: return the class name.
62  std::string name() const
63  { return "NormalizeComponentsWhitening"; }
64 
65  void train(ModelType& model, UnlabeledData<RealVector> const& input){
66  std::size_t dc = dataDimension(input);
67  SHARK_CHECK(input.numberOfElements() >= dc + 1, "[NormalizeComponentsWhitening::train] input needs to contain more points than there are input dimensions");
68  SHARK_CHECK(m_targetVariance > 0.0, "[NormalizeComponentsWhitening::train] target variance must be positive");
69 
70  // dense model with bias having input and output dimension equal to data dimension
71  model.setStructure(dc, dc, true);
72 
73  RealVector mean;
74  RealMatrix covariance;
75  meanvar(input, mean, covariance);
76 
77  RealMatrix whiteningMatrix = createWhiteningMatrix(covariance);
78  whiteningMatrix *= std::sqrt(m_targetVariance);
79 
80  RealVector offset = -prod(trans(whiteningMatrix),mean);
81 
82  model.setStructure(RealMatrix(trans(whiteningMatrix)), offset);
83  }
84 
85 private:
86  RealMatrix createWhiteningMatrix(
87  RealMatrix& covariance
88  ){
89  using namespace blas;
90 
91  SIZE_CHECK(covariance.size1() == covariance.size2());
92  std::size_t m = covariance.size1();
93  //we use the inversed cholesky decomposition for whitening
94  //since we have to assume that covariance does not have full rank, we use
95  //the generalized decomposition
96  RealMatrix whiteningMatrix(m,m,0.0);
97 
98  //do a pivoting cholesky decomposition
99  //this destroys the covariance matrix as it is not neeeded anymore afterwards.
100  PermutationMatrix permutation(m);
101  std::size_t rank = pivotingCholeskyDecompositionInPlace(covariance,permutation);
102  //only take the nonzero columns as C
103  matrix_range<RealMatrix> C = columns(covariance,0,rank);
104 
105  //full rank, means that we can use the typical cholesky inverse with pivoting
106  //so U is P C^-1 P^T
107  if(rank == m){
108  noalias(whiteningMatrix) = blas::identity_matrix<double>( m );
109  solveTriangularSystemInPlace<SolveXAB,upper>(trans(C),whiteningMatrix);
110  swap_full_inverted(permutation,whiteningMatrix);
111  return whiteningMatrix;
112  }
113  //complex case.
114  //A' = P C(C^TC)^-1(C^TC)^-1 C^T P^T
115  //=> P^T U P = C(C^TC)^-1
116  //<=> P^T U P (C^TC) = C
117  RealMatrix CTC(rank,rank);
118  symm_prod(trans(C),CTC);
119 
120  matrix_range<RealMatrix> submat = columns(whiteningMatrix,0,rank);
121  solveSymmPosDefSystem<SolveXAB>(CTC,submat,C);
122  swap_full_inverted(permutation,whiteningMatrix);
123 
124  return whiteningMatrix;
125  }
126 };
127 
128 
129 }
130 #endif