sp-matrix.h

Go to the documentation of this file.

// matrix/sp-matrix.h

// Copyright 2009-2011   Ondrej Glembek;  Microsoft Corporation;  Lukas Burget;
//                       Saarland University;  Ariya Rastrow;  Yanmin Qian;
//                       Jan Silovsky

// 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_SP_MATRIX_H_
#define KALDI_MATRIX_SP_MATRIX_H_

#include <algorithm>
#include <vector>

#include "matrix/packed-matrix.h"

namespace kaldi {


template<typename Real> class SpMatrix;


template<typename Real>
class SpMatrix : public PackedMatrix<Real> {
  friend class CuSpMatrix<Real>;
 public:
  // so it can use our assignment operator.
  friend class std::vector<Matrix<Real> >;

  SpMatrix(): PackedMatrix<Real>() {}


  explicit SpMatrix(const CuSpMatrix<Real> &cu);

  explicit SpMatrix(MatrixIndexT r, MatrixResizeType resize_type = kSetZero)
      : PackedMatrix<Real>(r, resize_type) {}

  SpMatrix(const SpMatrix<Real> &orig)
      : PackedMatrix<Real>(orig) {}

  template<typename OtherReal>
  explicit SpMatrix(const SpMatrix<OtherReal> &orig)
      : PackedMatrix<Real>(orig) {}

#ifdef KALDI_PARANOID
  explicit SpMatrix(const MatrixBase<Real> & orig,
                    SpCopyType copy_type = kTakeMeanAndCheck)
      : PackedMatrix<Real>(orig.NumRows(), kUndefined) {
    CopyFromMat(orig, copy_type);
  }
#else
  explicit SpMatrix(const MatrixBase<Real> & orig,
                    SpCopyType copy_type = kTakeMean)
      : PackedMatrix<Real>(orig.NumRows(), kUndefined) {
    CopyFromMat(orig, copy_type);
  }
#endif

  void Swap(SpMatrix *other);

  inline void Resize(MatrixIndexT nRows, MatrixResizeType resize_type = kSetZero) {
    PackedMatrix<Real>::Resize(nRows, resize_type);
  }

  void CopyFromSp(const SpMatrix<Real> &other) {
    PackedMatrix<Real>::CopyFromPacked(other);
  }

  template<typename OtherReal>
  void CopyFromSp(const SpMatrix<OtherReal> &other) {
    PackedMatrix<Real>::CopyFromPacked(other);
  }

#ifdef KALDI_PARANOID
  void CopyFromMat(const MatrixBase<Real> &orig,
                   SpCopyType copy_type = kTakeMeanAndCheck);
#else  // different default arg if non-paranoid mode.
  void CopyFromMat(const MatrixBase<Real> &orig,
                   SpCopyType copy_type = kTakeMean);
#endif

  inline Real operator() (MatrixIndexT r, MatrixIndexT c) const {
    // if column is less than row, then swap these as matrix is stored
    // as upper-triangular...  only allowed for const matrix object.
    if (static_cast<UnsignedMatrixIndexT>(c) >
        static_cast<UnsignedMatrixIndexT>(r))
      std::swap(c, r);
    // c<=r now so don't have to check c.
    KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
                 static_cast<UnsignedMatrixIndexT>(this->num_rows_));
    return *(this->data_ + (r*(r+1)) / 2 + c);
    // Duplicating code from PackedMatrix.h
  }

  inline Real &operator() (MatrixIndexT r, MatrixIndexT c) {
    if (static_cast<UnsignedMatrixIndexT>(c) >
        static_cast<UnsignedMatrixIndexT>(r))
      std::swap(c, r);
    // c<=r now so don't have to check c.
    KALDI_ASSERT(static_cast<UnsignedMatrixIndexT>(r) <
                 static_cast<UnsignedMatrixIndexT>(this->num_rows_));
    return *(this->data_ + (r * (r + 1)) / 2 + c);
    // Duplicating code from PackedMatrix.h
  }

  SpMatrix<Real>& operator=(const SpMatrix<Real> &other) {
    PackedMatrix<Real>::operator=(other);
    return *this;
  }

  using PackedMatrix<Real>::Scale;

  void Invert(Real *logdet = NULL, Real *det_sign= NULL,
              bool inverse_needed = true);

  // Below routine does inversion in double precision,
  // even for single-precision object.
  void InvertDouble(Real *logdet = NULL, Real *det_sign = NULL,
                    bool inverse_needed = true);

  inline Real Cond() const {
    Matrix<Real> tmp(*this);
    return tmp.Cond();
  }

  void ApplyPow(Real exponent);

  void SymPosSemiDefEig(VectorBase<Real> *s, MatrixBase<Real> *P,
                        Real tolerance = 0.001) const;

  void Eig(VectorBase<Real> *s, MatrixBase<Real> *P = NULL) const;

  void TopEigs(VectorBase<Real> *s, MatrixBase<Real> *P,
               MatrixIndexT lanczos_dim = 0) const;


  Real MaxAbsEig() const;

  void PrintEigs(const char *name) {
    Vector<Real> s((*this).NumRows());
    Matrix<Real> P((*this).NumRows(), (*this).NumCols());
    SymPosSemiDefEig(&s, &P);
    KALDI_LOG << "PrintEigs: " << name << ": " << s;
  }

  bool IsPosDef() const;  // returns true if Cholesky succeeds.
  void AddSp(const Real alpha, const SpMatrix<Real> &Ma) {
    this->AddPacked(alpha, Ma);
  }

  Real LogPosDefDet() const;

  Real LogDet(Real *det_sign = NULL) const;

  template<typename OtherReal>
  void AddVec2(const Real alpha, const VectorBase<OtherReal> &v);

  void AddVecVec(const Real alpha, const VectorBase<Real> &v,
                 const VectorBase<Real> &w);

  void AddVec2Sp(const Real alpha, const VectorBase<Real> &v,
                 const SpMatrix<Real> &S, const Real beta);

  template<typename OtherReal>
  void AddDiagVec(const Real alpha, const VectorBase<OtherReal> &v);

  void AddMat2(const Real alpha, const MatrixBase<Real> &M,
               MatrixTransposeType transM, const Real beta);

  void AddMat2Sp(const Real alpha, const MatrixBase<Real> &M,
                 MatrixTransposeType transM, const SpMatrix<Real> &A,
                 const Real beta = 0.0);

  void AddSmat2Sp(const Real alpha, const MatrixBase<Real> &M,
                  MatrixTransposeType transM, const SpMatrix<Real> &A,
                  const Real beta = 0.0);

  void AddTp2Sp(const Real alpha, const TpMatrix<Real> &T,
                MatrixTransposeType transM, const SpMatrix<Real> &A,
                const Real beta = 0.0);

  void AddTp2(const Real alpha, const TpMatrix<Real> &T,
              MatrixTransposeType transM, const Real beta = 0.0);

  void AddMat2Vec(const Real alpha, const MatrixBase<Real> &M,
                  MatrixTransposeType transM, const VectorBase<Real> &v,
                  const Real beta = 0.0);


  int ApplyFloor(const SpMatrix<Real> &Floor, Real alpha = 1.0,
                 bool verbose = false);

  int ApplyFloor(Real floor);

  bool IsDiagonal(Real cutoff = 1.0e-05) const;
  bool IsUnit(Real cutoff = 1.0e-05) const;
  bool IsZero(Real cutoff = 1.0e-05) const;
  bool IsTridiagonal(Real cutoff = 1.0e-05) const;

  Real FrobeniusNorm() const;

  bool ApproxEqual(const SpMatrix<Real> &other, float tol = 0.01) const;

  // LimitCond:
  // Limits the condition of symmetric positive semidefinite matrix to
  // a specified value
  // by flooring all eigenvalues to a positive number which is some multiple
  // of the largest one (or zero if there are no positive eigenvalues).
  // Takes the condition number we are willing to accept, and floors
  // eigenvalues to the largest eigenvalue divided by this.
  //  Returns #eigs floored or already equal to the floor.
  // Throws exception if input is not positive definite.
  // returns #floored.
  MatrixIndexT LimitCond(Real maxCond = 1.0e+5, bool invert = false);

  // as LimitCond but all done in double precision. // returns #floored.
  MatrixIndexT LimitCondDouble(Real maxCond = 1.0e+5, bool invert = false) {
    SpMatrix<double> dmat(*this);
    MatrixIndexT ans = dmat.LimitCond(maxCond, invert);
    (*this).CopyFromSp(dmat);
    return ans;
  }
  Real Trace() const;

  void Tridiagonalize(MatrixBase<Real> *Q);

  void Qr(MatrixBase<Real> *Q);

 private:
 void EigInternal(VectorBase<Real> *s, MatrixBase<Real> *P,
                   Real tolerance, int recurse) const;
};




float TraceSpSp(const SpMatrix<float> &A, const SpMatrix<float> &B);
double TraceSpSp(const SpMatrix<double> &A, const SpMatrix<double> &B);


template<typename Real>
inline bool ApproxEqual(const SpMatrix<Real> &A,
                        const SpMatrix<Real> &B, Real tol = 0.01) {
  return  A.ApproxEqual(B, tol);
}

template<typename Real>
inline void AssertEqual(const SpMatrix<Real> &A,
                        const SpMatrix<Real> &B, Real tol = 0.01) {
  KALDI_ASSERT(ApproxEqual(A, B, tol));
}



template<typename Real, typename OtherReal>
Real TraceSpSp(const SpMatrix<Real> &A, const SpMatrix<OtherReal> &B);



// TraceSpSpLower is the same as Trace(A B) except the lower-diagonal elements
// are counted only once not twice as they should be.  It is useful in certain
// optimizations.
template<typename Real>
Real TraceSpSpLower(const SpMatrix<Real> &A, const SpMatrix<Real> &B);


template<typename Real>
Real TraceSpMat(const SpMatrix<Real> &A, const MatrixBase<Real> &B);

template<typename Real>
Real TraceMatSpMat(const MatrixBase<Real> &A, MatrixTransposeType transA,
                   const SpMatrix<Real> &B, const MatrixBase<Real> &C,
                   MatrixTransposeType transC);

template<typename Real>
Real TraceMatSpMatSp(const MatrixBase<Real> &A, MatrixTransposeType transA,
                     const SpMatrix<Real> &B, const MatrixBase<Real> &C,
                     MatrixTransposeType transC, const SpMatrix<Real> &D);

template<typename Real>
Real VecSpVec(const VectorBase<Real> &v1, const SpMatrix<Real> &M,
               const VectorBase<Real> &v2);





struct SolverOptions {
  BaseFloat K; // maximum condition number
  BaseFloat eps;
  std::string name;
  bool optimize_delta;
  bool diagonal_precondition;
  bool print_debug_output;
  explicit SolverOptions(const std::string &name):
      K(1.0e+4), eps(1.0e-40), name(name),
      optimize_delta(true), diagonal_precondition(false),
      print_debug_output(true) { }
  SolverOptions(): K(1.0e+4), eps(1.0e-40), name("[unknown]"),
                   optimize_delta(true), diagonal_precondition(false),
                   print_debug_output(true) { }
  void Check() const;
};



template<typename Real>
Real SolveQuadraticProblem(const SpMatrix<Real> &H,
                           const VectorBase<Real> &g,
                           const SolverOptions &opts,
                           VectorBase<Real> *x);



template<typename Real>
Real SolveQuadraticMatrixProblem(const SpMatrix<Real> &Q,
                                 const MatrixBase<Real> &Y,
                                 const SpMatrix<Real> &P,
                                 const SolverOptions &opts,
                                 MatrixBase<Real> *M);

template<typename Real>
Real SolveDoubleQuadraticMatrixProblem(const MatrixBase<Real> &G,
                                       const SpMatrix<Real> &P1,
                                       const SpMatrix<Real> &P2,
                                       const SpMatrix<Real> &Q1,
                                       const SpMatrix<Real> &Q2,
                                       const SolverOptions &opts,
                                       MatrixBase<Real> *M);



}  // namespace kaldi


// Including the implementation (now actually just includes some
// template specializations).
#include "matrix/sp-matrix-inl.h"


#endif  // KALDI_MATRIX_SP_MATRIX_H_