add implementation for kaldi io (by ymz)

1 files changed, 166 insertions, 0 deletions

diff --git a/kaldi_io/src/kaldi/matrix/kaldi-gpsr.h b/kaldi_io/src/kaldi/matrix/kaldi-gpsr.h
new file mode 100644
index 0000000..c294bdd
--- /dev/null
+++ b/kaldi_io/src/kaldi/matrix/kaldi-gpsr.h
@@ -0,0 +1,166 @@
+// matrix/kaldi-gpsr.h
+
+// Copyright 2012  Arnab Ghoshal
+
+// See ../../COPYING for clarification regarding multiple authors
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//  http://www.apache.org/licenses/LICENSE-2.0
+//
+// THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
+// WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
+// MERCHANTABLITY OR NON-INFRINGEMENT.
+// See the Apache 2 License for the specific language governing permissions and
+// limitations under the License.
+
+#ifndef KALDI_MATRIX_KALDI_GPSR_H_
+#define KALDI_MATRIX_KALDI_GPSR_H_
+
+#include <string>
+#include <vector>
+
+#include "base/kaldi-common.h"
+#include "matrix/matrix-lib.h"
+#include "itf/options-itf.h"
+
+namespace kaldi {
+
+/// This is an implementation of the GPSR algorithm. See, Figueiredo, Nowak and
+/// Wright, "Gradient Projection for Sparse Reconstruction: Application to
+/// Compressed Sensing and Other Inverse Problems," IEEE Journal of Selected
+/// Topics in Signal Processing, vol. 1, no. 4, pp. 586-597, 2007.
+/// http://dx.doi.org/10.1109/JSTSP.2007.910281
+
+/// The GPSR algorithm, described in Figueiredo, et al., 2007, solves:
+/// \f[ \min_x 0.5 * ||y - Ax||_2^2 + \tau ||x||_1, \f]
+/// where \f$ x \in R^n, y \in R^k \f$, and \f$ A \in R^{n \times k} \f$.
+/// In this implementation, we solve:
+/// \f[ \min_x 0.5 * x^T H x - g^T x + \tau ||x||_1, \f]
+/// which is the more natural form in which such problems arise in our case.
+/// Here, \f$ H = A^T A \in R^{n \times n} \f$ and \f$ g = A^T y \in R^n \f$.
+
+
+/** \struct GpsrConfig
+ *  Configuration variables needed in the GPSR algorithm.
+ */
+struct GpsrConfig {
+  bool use_gpsr_bb;  ///< Use the Barzilai-Borwein gradient projection method
+
+  /// The following options are common to both the basic & Barzilai-Borwein
+  /// versions of GPSR
+  double stop_thresh;  ///< Stopping threshold
+  int32 max_iters;  ///< Maximum number of iterations
+  double gpsr_tau;  ///< Regularization scale
+  double alpha_min;  ///< Minimum step size in the feasible direction
+  double alpha_max;  ///< Maximum step size in the feasible direction
+  double max_sparsity;  ///< Maximum percentage of dimensions set to 0
+  double tau_reduction;  ///< Multiply tau by this if max_sparsity reached
+
+  /// The following options are for the backtracking line search in basic GPSR.
+  /// Step size reduction factor in backtracking line search. 0 < beta < 1
+  double gpsr_beta;
+  /// Improvement factor in backtracking line search, i.e. the new objective
+  /// function must be less than the old one by mu times the gradient in the
+  /// direction of the change in x. 0 < mu < 1
+  double gpsr_mu;
+  int32 max_iters_backtrak;  ///< Max iterations for backtracking line search
+
+  bool debias;  ///< Do debiasing, i.e. unconstrained optimization at the end
+  double stop_thresh_debias;  ///< Stopping threshold for debiasing stage
+  int32 max_iters_debias;  ///< Maximum number of iterations for debiasing stage
+
+  GpsrConfig() {
+    use_gpsr_bb = true;
+
+    stop_thresh = 0.005;
+    max_iters = 100;
+    gpsr_tau = 10;
+    alpha_min = 1.0e-10;
+    alpha_max = 1.0e+20;
+    max_sparsity = 0.9;
+    tau_reduction = 0.8;
+
+    gpsr_beta = 0.5;
+    gpsr_mu = 0.1;
+    max_iters_backtrak = 50;
+
+    debias = false;
+    stop_thresh_debias = 0.001;
+    max_iters_debias = 50;
+  }
+
+  void Register(OptionsItf *po);
+};
+
+inline void GpsrConfig::Register(OptionsItf *po) {
+  std::string module = "GpsrConfig: ";
+  po->Register("use-gpsr-bb", &use_gpsr_bb, module+
+               "Use the Barzilai-Borwein gradient projection method.");
+
+  po->Register("stop-thresh", &stop_thresh, module+
+               "Stopping threshold for GPSR.");
+  po->Register("max-iters", &max_iters, module+
+               "Maximum number of iterations of GPSR.");
+  po->Register("gpsr-tau", &gpsr_tau, module+
+               "Regularization scale for GPSR.");
+  po->Register("alpha-min", &alpha_min, module+
+               "Minimum step size in feasible direction.");
+  po->Register("alpha-max", &alpha_max, module+
+               "Maximum step size in feasible direction.");
+  po->Register("max-sparsity", &max_sparsity, module+
+               "Maximum percentage of dimensions set to 0.");
+  po->Register("tau-reduction", &tau_reduction, module+
+               "Multiply tau by this if maximum sparsity is reached.");
+
+  po->Register("gpsr-beta", &gpsr_beta, module+
+               "Step size reduction factor in backtracking line search (0<beta<1).");
+  po->Register("gpsr-mu", &gpsr_mu, module+
+               "Improvement factor in backtracking line search (0<mu<1).");
+  po->Register("max-iters-backtrack", &max_iters_backtrak, module+
+               "Maximum number of iterations of backtracking line search.");
+
+  po->Register("debias", &debias, module+
+               "Do final debiasing step.");
+  po->Register("stop-thresh-debias", &stop_thresh_debias, module+
+               "Stopping threshold for debiaisng step.");
+  po->Register("max-iters-debias", &max_iters_debias, module+
+               "Maximum number of iterations of debiasing.");
+}
+
+/// Solves a quadratic program in \f$ x \f$, with L_1 regularization:
+/// \f[ \min_x 0.5 * x^T H x - g^T x + \tau ||x||_1. \f]
+/// This is similar to SolveQuadraticProblem() in sp-matrix.h with an added
+/// L_1 term.
+template<typename Real>
+Real Gpsr(const GpsrConfig &opts, const SpMatrix<Real> &H,
+          const Vector<Real> &g, Vector<Real> *x,
+          const char *debug_str = "[unknown]") {
+  if (opts.use_gpsr_bb)
+    return GpsrBB(opts, H, g, x, debug_str);
+  else
+    return GpsrBasic(opts, H, g, x, debug_str);
+}
+
+/// This is the basic GPSR algorithm, where the step size is determined by a
+/// backtracking line search. The line search is called "Armijo rule along the
+/// projection arc" in Bertsekas, Nonlinear Programming, 2nd ed. page 230.
+template<typename Real>
+Real GpsrBasic(const GpsrConfig &opts, const SpMatrix<Real> &H,
+               const Vector<Real> &g, Vector<Real> *x,
+               const char *debug_str = "[unknown]");
+
+/// This is the paper calls the Barzilai-Borwein variant. This is a constrained
+/// Netwon's method where the Hessian is approximated by scaled identity matrix
+template<typename Real>
+Real GpsrBB(const GpsrConfig &opts, const SpMatrix<Real> &H,
+            const Vector<Real> &g, Vector<Real> *x,
+            const char *debug_str = "[unknown]");
+
+
+}  // namespace kaldi
+
+#endif  // KALDI_MATRIX_KALDI_GPSR_H_