All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
est-pca.cc File Reference
Include dependency graph for est-pca.cc:

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 25 of file est-pca.cc.

References SpMatrix< Real >::AddMat2(), VectorBase< Real >::AddMatVec(), VectorBase< Real >::AddRowSumMat(), VectorBase< Real >::AddVec(), SpMatrix< Real >::AddVec2(), MatrixBase< Real >::CopyColFromVec(), count, VectorBase< Real >::Dim(), SequentialTableReader< Holder >::Done(), SpMatrix< Real >::Eig(), ParseOptions::GetArg(), rnnlm::i, KALDI_ERR, KALDI_LOG, KALDI_WARN, kaldi::kCopyData, kaldi::kNoTrans, kaldi::kTrans, SequentialTableReader< Holder >::Next(), ParseOptions::NumArgs(), MatrixBase< Real >::NumCols(), ParseOptions::PrintUsage(), ParseOptions::Read(), ParseOptions::Register(), SpMatrix< Real >::Resize(), Vector< Real >::Resize(), Matrix< Real >::Resize(), MatrixBase< Real >::Row(), PackedMatrix< Real >::Scale(), VectorBase< Real >::Scale(), kaldi::SortSvd(), SequentialTableReader< Holder >::Value(), and kaldi::WriteKaldiObject().

25  {
26  using namespace kaldi;
27  typedef kaldi::int32 int32;
28  try {
29  const char *usage =
30  "Estimate PCA transform; dimension reduction is optional (if not specified\n"
31  "we don't reduce the dimension; if you specify --normalize-variance=true,\n"
32  "we normalize the (centered) covariance of the features, and if you specify\n"
33  "--normalize-mean=true the mean is also normalized. So a variety of transform\n"
34  "types are supported. Because this type of transform does not need too much\n"
35  "data to estimate robustly, we don't support separate accumulator files;\n"
36  "this program reads in the features directly. For large datasets you may\n"
37  "want to subset the features (see example below)\n"
38  "By default the program reads in matrices (e.g. features), but with\n"
39  "--read-vectors=true, can read in vectors (e.g. iVectors).\n"
40  "\n"
41  "Usage: est-pca [options] (<feature-rspecifier>|<vector-rspecifier>) <pca-matrix-out>\n"
42  "e.g.:\n"
43  "utils/shuffle_list.pl data/train/feats.scp | head -n 5000 | sort | \\\n"
44  " est-pca --dim=50 scp:- some/dir/0.mat\n";
45 
46  bool binary = true;
47  bool read_vectors = false;
48  bool normalize_variance = false;
49  bool normalize_mean = false;
50  int32 dim = -1;
51  std::string full_matrix_wxfilename;
52  ParseOptions po(usage);
53  po.Register("binary", &binary, "Write accumulators in binary mode.");
54  po.Register("dim", &dim, "Feature dimension requested (if <= 0, uses full "
55  "feature dimension");
56  po.Register("read-vectors", &read_vectors, "If true, read in single vectors "
57  "instead of feature matrices");
58  po.Register("normalize-variance", &normalize_variance, "If true, make a "
59  "transform that normalizes variance to one.");
60  po.Register("normalize-mean", &normalize_mean, "If true, output an affine "
61  "transform that subtracts the data mean.");
62  po.Register("write-full-matrix", &full_matrix_wxfilename,
63  "Write full version of the matrix to this location (including "
64  "rejected rows)");
65  po.Read(argc, argv);
66 
67  if (po.NumArgs() != 2) {
68  po.PrintUsage();
69  exit(1);
70  }
71 
72  std::string rspecifier = po.GetArg(1),
73  pca_mat_wxfilename = po.GetArg(2);
74 
75  int32 num_done = 0, num_err = 0;
76  int64 count = 0;
77  Vector<double> sum;
78  SpMatrix<double> sumsq;
79 
80  if (!read_vectors) {
81  SequentialBaseFloatMatrixReader feat_reader(rspecifier);
82 
83  for (; !feat_reader.Done(); feat_reader.Next()) {
84  Matrix<double> mat(feat_reader.Value());
85  if (mat.NumRows() == 0) {
86  KALDI_WARN << "Empty feature matrix";
87  num_err++;
88  continue;
89  }
90  if (sum.Dim() == 0) {
91  sum.Resize(mat.NumCols());
92  sumsq.Resize(mat.NumCols());
93  }
94  if (sum.Dim() != mat.NumCols()) {
95  KALDI_WARN << "Feature dimension mismatch " << sum.Dim() << " vs. "
96  << mat.NumCols();
97  num_err++;
98  continue;
99  }
100  sum.AddRowSumMat(1.0, mat);
101  sumsq.AddMat2(1.0, mat, kTrans, 1.0);
102  count += mat.NumRows();
103  num_done++;
104  }
105  KALDI_LOG << "Accumulated stats from " << num_done << " feature files, "
106  << num_err << " with errors; " << count << " frames.";
107  } else {
108  // read in vectors, not matrices
109  SequentialBaseFloatVectorReader vec_reader(rspecifier);
110 
111  for (; !vec_reader.Done(); vec_reader.Next()) {
112  Vector<double> vec(vec_reader.Value());
113  if (vec.Dim() == 0) {
114  KALDI_WARN << "Empty input vector";
115  num_err++;
116  continue;
117  }
118  if (sum.Dim() == 0) {
119  sum.Resize(vec.Dim());
120  sumsq.Resize(vec.Dim());
121  }
122  if (sum.Dim() != vec.Dim()) {
123  KALDI_WARN << "Feature dimension mismatch " << sum.Dim() << " vs. "
124  << vec.Dim();
125  num_err++;
126  continue;
127  }
128  sum.AddVec(1.0, vec);
129  sumsq.AddVec2(1.0, vec);
130  count += 1.0;
131  num_done++;
132  }
133  KALDI_LOG << "Accumulated stats from " << num_done << " vectors, "
134  << num_err << " with errors.";
135  }
136  if (num_done == 0)
137  KALDI_ERR << "No data accumulated.";
138  sum.Scale(1.0 / count);
139  sumsq.Scale(1.0 / count);
140 
141  sumsq.AddVec2(-1.0, sum); // now sumsq is centered covariance.
142 
143  int32 full_dim = sum.Dim();
144  if (dim <= 0) dim = full_dim;
145  if (dim > full_dim)
146  KALDI_ERR << "Final dimension " << dim << " is greater than feature "
147  << "dimension " << full_dim;
148 
149  Matrix<double> P(full_dim, full_dim);
150  Vector<double> s(full_dim);
151 
152  sumsq.Eig(&s, &P);
153  SortSvd(&s, &P);
154 
155  KALDI_LOG << "Eigenvalues in PCA are " << s;
156  KALDI_LOG << "Sum of PCA eigenvalues is " << s.Sum() << ", sum of kept "
157  << "eigenvalues is " << s.Range(0, dim).Sum();
158 
159 
160  Matrix<double> transform(P, kTrans); // Transpose of P. This is what
161  // appears in the transform.
162 
163  if (normalize_variance) {
164  for (int32 i = 0; i < full_dim; i++) {
165  double this_var = s(i), min_var = 1.0e-15;
166  if (this_var < min_var) {
167  KALDI_WARN << "--normalize-variance option: very tiny variance " << s(i)
168  << "encountered, treating as " << min_var;
169  this_var = min_var;
170  }
171  double scale = 1.0 / sqrt(this_var); // scale on features that will make
172  // the variance unit.
173  transform.Row(i).Scale(scale);
174  }
175  }
176 
177  Vector<double> offset(full_dim);
178 
179  if (normalize_mean) {
180  offset.AddMatVec(-1.0, transform, kNoTrans, sum, 0.0);
181  transform.Resize(full_dim, full_dim + 1, kCopyData); // Add column to transform.
182  transform.CopyColFromVec(offset, full_dim);
183  }
184 
185  Matrix<BaseFloat> transform_float(transform);
186 
187  if (full_matrix_wxfilename != "") {
188  WriteKaldiObject(transform_float, full_matrix_wxfilename, binary);
189  }
190 
191  transform_float.Resize(dim, transform_float.NumCols(), kCopyData);
192  WriteKaldiObject(transform_float, pca_mat_wxfilename, binary);
193 
194  return 0;
195  } catch(const std::exception &e) {
196  std::cerr << e.what();
197  return -1;
198  }
199 }
Relabels neural network egs with the read pdf-id alignments.
Definition: chain.dox:20
void AddRowSumMat(Real alpha, const MatrixBase< Real > &M, Real beta=1.0)
Does *this = alpha * (sum of rows of M) + beta * *this.
void AddMat2(const Real alpha, const MatrixBase< Real > &M, MatrixTransposeType transM, const Real beta)
rank-N update: if (transM == kNoTrans) (*this) = beta*(*this) + alpha * M * M^T, or (if transM == kTr...
Definition: sp-matrix.cc:1134
void Scale(Real c)
void Resize(MatrixIndexT length, MatrixResizeType resize_type=kSetZero)
Set vector to a specified size (can be zero).
const size_t count
The class ParseOptions is for parsing command-line options; see Parsing command-line options for more...
Definition: parse-options.h:36
void Resize(MatrixIndexT nRows, MatrixResizeType resize_type=kSetZero)
Definition: sp-matrix.h:81
A templated class for reading objects sequentially from an archive or script file; see The Table conc...
Definition: kaldi-table.h:287
#define KALDI_ERR
Definition: kaldi-error.h:127
#define KALDI_WARN
Definition: kaldi-error.h:130
void Scale(Real alpha)
Multiplies all elements by this constant.
void AddVec2(const Real alpha, const VectorBase< OtherReal > &v)
rank-one update, this <– this + alpha v v'
Definition: sp-matrix.cc:970
void WriteKaldiObject(const C &c, const std::string &filename, bool binary)
Definition: kaldi-io.h:257
void Eig(VectorBase< Real > *s, MatrixBase< Real > *P=NULL) const
Solves the symmetric eigenvalue problem: at end we should have (*this) = P * diag(s) * P^T...
Definition: qr.cc:433
#define KALDI_LOG
Definition: kaldi-error.h:133
void AddVec(const Real alpha, const VectorBase< OtherReal > &v)
Add vector : *this = *this + alpha * rv (with casting between floats and doubles) ...
MatrixIndexT Dim() const
Returns the dimension of the vector.
Definition: kaldi-vector.h:59
void SortSvd(VectorBase< Real > *s, MatrixBase< Real > *U, MatrixBase< Real > *Vt, bool sort_on_absolute_value)
Function to ensure that SVD is sorted.