// 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_