path: root/kaldi_io/src/kaldi/matrix/sp-matrix.h

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


/// \addtogroup matrix_group
/// @{
template<typename Real> class SpMatrix;


/**
 * @brief Packed symetric matrix class
*/
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>() {}

  /// Copy constructor from CUDA version of SpMatrix
  /// This is defined in ../cudamatrix/cu-sp-matrix.h
  
  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

  /// Shallow swap.
  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
  }

  using PackedMatrix<Real>::operator =;
  using PackedMatrix<Real>::Scale;

  /// matrix inverse.
  /// if inverse_needed = false, will fill matrix with garbage.
  /// (only useful if logdet wanted).
  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);

  /// Returns maximum ratio of singular values.
  inline Real Cond() const {
    Matrix<Real> tmp(*this);
    return tmp.Cond();
  }

  /// Takes matrix to a fraction power via Svd.
  /// Will throw exception if matrix is not positive semidefinite
  /// (to within a tolerance)
  void ApplyPow(Real exponent);

  /// This is the version of SVD that we implement for symmetric positive
  /// definite matrices.  This exists for historical reasons; right now its
  /// internal implementation is the same as Eig().  It computes the eigenvalue
  /// decomposition (*this) = P * diag(s) * P^T with P orthogonal.  Will throw
  /// exception if input is not positive semidefinite to within a tolerance.
  void SymPosSemiDefEig(VectorBase<Real> *s, MatrixBase<Real> *P,
                        Real tolerance = 0.001) const;

  /// Solves the symmetric eigenvalue problem: at end we should have (*this) = P
  /// * diag(s) * P^T.  We solve the problem using the symmetric QR method.
  /// P may be NULL.
  /// Implemented in qr.cc.
  /// If you need the eigenvalues sorted, the function SortSvd declared in
  /// kaldi-matrix is suitable.
  void Eig(VectorBase<Real> *s, MatrixBase<Real> *P = NULL) const;
  
  /// This function gives you, approximately, the largest eigenvalues of the
  /// symmetric matrix and the corresponding eigenvectors.  (largest meaning,
  /// further from zero).  It does this by doing a SVD within the Krylov
  /// subspace generated by this matrix and a random vector.  This is
  /// a form of the Lanczos method with complete reorthogonalization, followed
  /// by SVD within a smaller dimension ("lanczos_dim").
  ///
  /// If *this is m by m, s should be of dimension n and P should be of
  /// dimension m by n, with n <= m.  The *columns* of P are the approximate
  /// eigenvectors; P * diag(s) * P^T would be a low-rank reconstruction of
  /// *this.  The columns of P will be orthogonal, and the elements of s will be
  /// the eigenvalues of *this projected into that subspace, but beyond that
  /// there are no exact guarantees.  (This is because the convergence of this
  /// method is statistical).  Note: it only makes sense to use this
  /// method if you are in very high dimension and n is substantially smaller
  /// than m: for example, if you want the 100 top eigenvalues of a 10k by 10k
  /// matrix.  This function calls Rand() to initialize the lanczos
  /// iterations and also for restarting.
  /// If lanczos_dim is zero, it will default to the greater of:
  /// s->Dim() + 50 or s->Dim() + s->Dim()/2, but not more than this->Dim().
  /// If lanczos_dim == this->Dim(), you might as well just call the function
  /// Eig() since the result will be the same, and Eig() would be faster; the
  /// whole point of this function is to reduce the dimension of the SVD
  /// computation.
  void TopEigs(VectorBase<Real> *s, MatrixBase<Real> *P,
               MatrixIndexT lanczos_dim = 0) const;


  
  /// Takes log of the matrix (does eigenvalue decomposition then takes
  /// log of eigenvalues and reconstructs).  Will throw of not +ve definite.
  void Log();


  // Takes exponential of the matrix (equivalent to doing eigenvalue
  // decomposition then taking exp of eigenvalues and reconstructing).
  void Exp();

  /// Returns the maximum of the absolute values of any of the
  /// eigenvalues.
  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);
  }

  /// Computes log determinant but only for +ve-def matrices
  /// (it uses Cholesky).
  /// If matrix is not +ve-def, it will throw an exception
  /// was LogPDDeterminant()
  Real LogPosDefDet() const;

  Real LogDet(Real *det_sign = NULL) const;

  /// rank-one update, this <-- this + alpha v v'
  template<typename OtherReal>
  void AddVec2(const Real alpha, const VectorBase<OtherReal> &v);

  /// rank-two update, this <-- this + alpha (v w' + w v').
  void AddVecVec(const Real alpha, const VectorBase<Real> &v,
                 const VectorBase<Real> &w);

  /// Does *this = beta * *thi + alpha * diag(v) * S * diag(v)
  void AddVec2Sp(const Real alpha, const VectorBase<Real> &v,
                 const SpMatrix<Real> &S, const Real beta);
  
  /// diagonal update, this <-- this + diag(v)
  template<typename OtherReal>
  void AddDiagVec(const Real alpha, const VectorBase<OtherReal> &v);

  /// rank-N update:
  /// if (transM == kNoTrans)
  /// (*this) = beta*(*this) + alpha * M * M^T,
  /// or  (if transM == kTrans)
  ///  (*this) = beta*(*this) + alpha * M^T * M
  /// Note: beta used to default to 0.0.