Index: trunk/npstat/nm/Matrix.hh =================================================================== --- trunk/npstat/nm/Matrix.hh (revision 707) +++ trunk/npstat/nm/Matrix.hh (revision 708) @@ -1,826 +1,826 @@ #ifndef NPSTAT_MATRIX_HH_ #define NPSTAT_MATRIX_HH_ /*! // \file Matrix.hh // // \brief Template matrix class // // Author: I. Volobouev // // November 2008 */ #include #include #include #include #include "geners/ClassId.hh" #include "npstat/nm/EigenMethod.hh" #include "npstat/nm/SvdMethod.hh" #include "npstat/nm/GeneralizedComplex.hh" namespace npstat { /** // A simple helper class for matrix manipulations. Depending on how much // space is provided with the "Len" parameter, the data will be placed // either on the stack or on the heap. Do not set "Len" to 0, the code // assumes that the stack space is available for at least one element. // // C storage convention is used for internal data. In the element access // operations, array bounds are not checked. // // A number of operations (solving linear systems, calculating eigenvalues // and eigenvectors, singular value decomposition, etc) are performed // by calling appropriate routines from LAPACK. This usually limits // the Numeric template parameter types for which these operations are // available to float and double. If the operation is unavailable for // the template parameter used, std::invalid_argument exception is raised. // // Note that for simple matrix operations this class is likely to be // slower than matrix classes based on expression templates (such as // those in boost uBLAS or in Blitz++). If speed is really important // for your calculations, consider using a dedicated matrix library. */ template class Matrix { template friend class Matrix; public: typedef Numeric value_type; /** // Default constructor creates an unitialized matrix // which can be assigned from other matrices */ Matrix(); /** // This constructor creates an unitialized matrix // which can be assigned element-by-element or from // another matrix with the same dimensions */ Matrix(unsigned nrows, unsigned ncols); /** // This constructor initializes the matrix as follows: // // initCode = 0: all elements are initialized to 0. // // initCode = 1: matrix must be square; diagonal elements are // initialized to 1, all other elements to 0. */ Matrix(unsigned nrows, unsigned ncols, int initCode); /** // This constructor initializes the matrix from the // given 1-d array using C storage conventions */ Matrix(unsigned nrows, unsigned ncols, const Numeric* data); /** Copy constructor */ Matrix(const Matrix&); /** Move constructor */ Matrix(Matrix&&); /** Converting copy constructor */ template Matrix(const Matrix&); /** // Constructor from a subrange of another matrix. The minimum row // and column numbers are included, maximum numbers are excluded. */ template Matrix(const Matrix&, unsigned rowMin, unsigned rowMax, unsigned columnMin, unsigned columnMax); ~Matrix(); /** // Assignment operator. The matrix on the left must either be // in an uninitialized state or have compatible dimensions with // the matrix on the right. */ Matrix& operator=(const Matrix&); /** Move assignment operator */ Matrix& operator=(Matrix&&); /** Converting assignment operator */ template Matrix& operator=(const Matrix&); /** Set from triplets (for compatibility with Eigen) */ template Matrix& setFromTriplets(Iterator first, Iterator last); /** Set all data from an external array */ template Matrix& setData(const Num2* data, unsigned dataLength); /** Set a row from an external array */ template - Matrix& setRow( unsigned row, const Num2* data, unsigned dataLength); + Matrix& setRow(unsigned row, const Num2* data, unsigned dataLength); /** Set a column from an external array */ template Matrix& setColumn(unsigned col, const Num2* data, unsigned dataLength); /** // Tag this matrix as diagonal. This may improve the speed // of certain operations. Of course, this should be done only // if you know for sure that this matrix is, indeed, diagonal. // Works for square matrices only. // // This method can also be used to tag matrices as non-diagonal. */ Matrix& tagAsDiagonal(bool b=true); /** // Fill the main diagonal and set all other elements to zero. // This operation tags the matrix as diagonal. */ template Matrix& diagFill(const Num2* data, unsigned dataLen); //@{ /** A simple inspection of matrix properties */ inline unsigned nRows() const {return nrows_;} inline unsigned nColumns() const {return ncols_;} inline unsigned long length() const {return static_cast(nrows_)*ncols_;} inline Numeric* data() const {return data_;} inline bool isSquare() const {return data_ && ncols_ == nrows_;} bool isSymmetric() const; bool isAntiSymmetric() const; bool isDiagonal() const; bool isZero() const; //@} /** // This method resets the object to an unintialized state // Can be applied in order to force the assignment operators to work. */ Matrix& uninitialize(); /** // Check if this matrix has the same number of rows and columns // as the other one. "false" will be returned in case any one of // the two matrices (or both) is not initialized. */ bool isCompatible(const Matrix& other) const; /** // This method changes the object dimensions. // All data is lost in the process. */ Matrix& resize(unsigned nrows, unsigned ncols); /** This method sets all elements to 0 */ Matrix& zeroOut(); /** // This method sets all diagonal elements to 0. // It can only be used with square matrices. */ Matrix& clearMainDiagonal(); /** // This method sets all non-diagonal elements to 0. // This operation is only valid for square matrices. // It tags the matrix as diagonal. */ Matrix& makeDiagonal(); /** This method sets all elements to the given value */ Matrix& constFill(Numeric c); /** // Method for transposing this matrix. Uses std::swap for square // matrices but might allocate memory from the heap if the matrix // is not square. */ Matrix& Tthis(); //@{ /** Compare two matrices for equality */ bool operator==(const Matrix&) const; bool operator!=(const Matrix&) const; //@} /** // Non-const access to the data (works like 2-d array in C). // No bounds checking. */ Numeric* operator[](unsigned); /** // Const access to the data (works like 2-d array in C). // No bounds checking. */ const Numeric* operator[](unsigned) const; /** Data modification method with bounds checking */ Matrix& set(unsigned row, unsigned column, Numeric value); /** Access by value without bounds checking */ Numeric operator()(unsigned row, unsigned column) const; /** Access by value with bounds checking */ Numeric at(unsigned row, unsigned column) const; /** Sum of the values in the given row */ Numeric rowSum(unsigned row) const; /** Sum of the values in the given column */ Numeric columnSum(unsigned column) const; /** Remove row and/or column. Indices out of range are ignored */ Matrix removeRowAndColumn(unsigned row, unsigned column) const; /** Number of non-zero elements */ unsigned nonZeros() const; /** // Replace every n x m square block in the matrix by its sum. // The number of rows must be divisible by n. The number of // columns must be divisible by m. */ void coarseSum(unsigned n, unsigned m, Matrix* result) const; /** // Replace every n x m square block in the matrix by its average // The number of rows must be divisible by n. The number of // columns must be divisible by m. */ void coarseAverage(unsigned n, unsigned m, Matrix* result) const; /** Unary plus */ Matrix operator+() const; /** Unary minus */ Matrix operator-() const; //@{ /** Binary algebraic operation with matrix or scalar */ Matrix operator*(const Matrix& r) const; Matrix operator*(Numeric r) const; Matrix operator/(Numeric r) const; Matrix operator+(const Matrix& r) const; Matrix operator-(const Matrix& r) const; //@} /** Hadamard product (a.k.a. Schur product) */ Matrix hadamardProduct(const Matrix& r) const; /** Hadamard ratio. All elements of the denominator must be non-zero */ Matrix hadamardRatio(const Matrix& denominator) const; //@{ /** // Binary algebraic operation which writes its result into // an existing matrix potentially avoiding memory allocation */ void times(const Matrix& r, Matrix* result) const; void times(Numeric r, Matrix* result) const; void over(Numeric r, Matrix* result) const; void plus(const Matrix& r, Matrix* result) const; void minus(const Matrix& r, Matrix* result) const; //@} //@{ /** In-place algebraic operations with matrices and scalars */ Matrix& operator*=(Numeric r); Matrix& operator/=(Numeric r); Matrix& operator+=(const Matrix& r); Matrix& operator-=(const Matrix& r); //@} //@{ /** Multiplication by a vector represented by a pointer and array size */ template Matrix timesVector(const Num2* data, unsigned dataLen) const; template void timesVector(const Num2* data, unsigned dataLen, Num3* result, unsigned resultLen) const; //@} /** // Multiplication by a vector represented by a pointer and array size. // This function is useful when only one element of the result is needed // (or elements are needed one at a time). */ template Numeric timesVector(unsigned rowNumber, const Num2* data, unsigned dataLen) const; //@{ /** Multiplication by a row on the left */ template Matrix rowMultiply(const Num2* data, unsigned dataLen) const; template void rowMultiply(const Num2* data, unsigned dataLen, Num3* result, unsigned resultLen) const; //@} /** // Multiplication by a row on the left. This function is useful when // only one element of the result is needed (or elements are needed // one at a time). */ template Numeric rowMultiply(unsigned columnNumber, const Num2* data, unsigned dataLen) const; /** Bilinear form with a vector (v^T M v) */ template Numeric bilinear(const Num2* data, unsigned dataLen) const; /** Bilinear form with a matrix (r^T M r) */ Matrix bilinear(const Matrix& r) const; /** Bilinear form with a transposed matrix (r M r^T) */ Matrix bilinearT(const Matrix& r) const; /** // Solution of a linear system of equations M*x = rhs. // The matrix must be square. "false" is returned in case // the matrix is degenerate. */ bool solveLinearSystem(const Numeric* rhs, unsigned lenRhs, Numeric* solution) const; /** // Solution of linear systems of equations M*X = RHS. // The matrix M (this one) must be square. "false" is // returned in case this matrix is degenerate. */ bool solveLinearSystems(const Matrix& RHS, Matrix* X) const; /** // Solution of an underdetermined linear system of equations M*x = rhs // in the minimum norm sense (M is this matrix). The norm is defined by // sqrt(x^T*V^(-1)*x). The result is given by A*rhs, where A is the // right inverse of M (i.e., M*A = I) minimizing the norm. V plays // the role of the inverse Gram matrix. It must be symmetric and // positive-definite. // // This method can also be used for solving constrained least squares // problems. For example, minimization of (y - y0)^T*V^(-1)*(y - y0) // under the condition M*y = rhs can be reduced to solving an // underdetermined linear system using x = y - y0 and // M*x = rhs' = rhs - M*y0. The application code can call this function // assuming that V is the covariance matrix of y and using rhs' as the // "rhs" argument. Subsequently, y can be calculated as x + y0. // The covariance matrix of y can be subsequently calculated by // V_y = (I - P)*V*(I - P)^T, where P = A*M. Operator (I - P) actually // projects the x (i.e., y - y0) space onto the space of solutions of // the equation M*x = 0. // // "false" is returned in case of failure. */ bool underdeterminedLinearSystem(const Numeric* rhs, unsigned lenRhs, const Matrix& V, Numeric* solution, unsigned lenSolution, Numeric* resultNormSquared = 0, Matrix* A = 0) const; /** // Solution of a linear overdetermined system of equations M*x = rhs // in the least squares sense. The "lenRhs" parameter should be equal // to the number of matrix rows, and it should exceed "lenSolution" // ("lenSolution" should be equal to the number of columns). "false" // is returned in case of failure. */ bool linearLeastSquares(const Numeric* rhs, unsigned lenRhs, Numeric* solution, unsigned lenSolution) const; /** // Solution of a linear system of equations M*x = rhs1 in the least // squares sense, subject to the constraint B*x = rhs2. Number of // equations (i.e., dimensionality of rhs1) plus the number of // constraints (i.e., dimensionality of rhs2) should exceed the // dimensionality of x. "false" is returned in case of failure. If // the chi-square of the result is desired, make the "resultChiSquare" // argument point to some valid location. Same goes for the result // covariance matrix and all subsequent arguments which can be used // to extract various intermediate results: // // unconstrainedSolution -- Array to store the unconstrained solution. // Must have at least "lenSol" elements. // // unconstrainedCovmat -- Covariance matrix of the unconstrained // solution. // // projectionMatrix -- Projection matrix used to obtain a part // of the constrained solution from the // unconstrained one. // // A -- Matrix multiplied by rhs2 to get the second // part of the constrained solution. // // Note that requesting the result covariance matrix or any other // subsequent output switches the code to a different and slower // "textbook" algorithm instead of an optimized algorithm from Lapack. */ bool constrainedLeastSquares(const Numeric* rhs1, unsigned lenRhs1, const Matrix& B, const Numeric* rhs2, unsigned lenRhs2, Numeric* solution, unsigned lenSol, Numeric* resultChiSquare = 0, Matrix* resultCovarianceMatrix = 0, Numeric* unconstrainedSolution = 0, Matrix* unconstrainedCovmat = 0, Matrix* projectionMatrix = 0, Matrix* A = 0) const; /** // Weighted least squares problem. The inverse covariance matrix // must be symmetric and positive definite. If the chi-square of the // result is desired, make the "resultChiSquare" argument point to // some valid location. Same goes for the result covariance matrix. // "false" is returned in case of failure. */ bool weightedLeastSquares(const Numeric* rhs, unsigned lenRhs, const Matrix& inverseCovarianceMatrix, Numeric* solution, unsigned lenSolution, Numeric* resultChiSquare = 0, Matrix* resultCovarianceMatrix = 0) const; /** Return transposed matrix */ Matrix T() const; /** Return the product of this matrix transposed with this, M^T*M */ Matrix TtimesThis() const; /** Return the product of this matrix with its transpose, M*M^T */ Matrix timesT() const; /** Return the product of this matrix with the transposed argument */ Matrix timesT(const Matrix& r) const; /** Return the product of this matrix transposed with the argument */ Matrix Ttimes(const Matrix& r) const; /** // "directSum" method constructs a block-diagonal matrix with // this matrix and the argument matrix inside the diagonal blocks. // Block which contains this matrix is the top left one. */ Matrix directSum(const Matrix& added) const; /** Construct a symmetric matrix from this one (symmetrize) */ Matrix symmetrize() const; /** Construct an antisymmetric matrix from this one (antisymmetrize) */ Matrix antiSymmetrize() const; /** // Matrix outer product. The number of rows of the result equals // the product of the numbers of rows of this matrix and the argument // matrix. The number of columns of the result is the product of the // numbers of columns. */ Matrix outer(const Matrix& r) const; /** Minimum value for any row/column */ Numeric minValue() const; /** Maximum value for any row/column */ Numeric maxValue() const; /** Maximum absolute value for any row/column */ Numeric maxAbsValue() const; /** Row and column of the smallest matrix element */ std::pair argmin() const; /** Row and column of the largest matrix element */ std::pair argmax() const; /** // Row and column of the matrix element with the largest // absolute value */ std::pair argmaxAbs() const; /** Frobenius norm of this matrix */ double frobeniusNorm() const; //@{ /** Calculate the trace (matrix must be square) */ Numeric tr() const; inline Numeric sp() const {return tr();} //@} //@{ /** // Calculate the trace of a product of this matrix with another. // This works faster if only the trace is of interest but not the // product itself. */ Numeric productTr(const Matrix& r) const; inline Numeric productSp(const Matrix& r) const {return productTr(r);} //@} /** Calculate the determinant (matrix must be square) */ Numeric det() const; /** // Inverse of a real symmetric positive-definite matrix. Use this // function (it has better numerical precision than symInv()) if // you know that your matrix is symmetric and positive definite. // For example, a non-singular covariance matrix calculated for // some dataset has these properties. */ Matrix symPDInv() const; /** // Invert real symmetric positive-definite matrix using // eigendecomposition. Use this function if you know that your // matrix is supposed to be symmetric and positive-definite but // might not be due to numerical round-offs. The arguments of // this method are as folows: // // tol -- Eigenvalues whose ratios to the largest eigenvalue // are below this parameter will be either truncated // or extended (per value of the next argument). // This parameter must lie in (0, 1). // // extend -- If "true", the inverse of small eigenvalues will // be replaced by the inverse of the largest eigenvalue // divided by the tolerance parameter. If "false", // inverse of such eigenvalues will be set to 0. // // m -- LAPACK method to use for calculating the // eigendecomposition. */ Matrix symPDEigenInv(double tol, bool extend=true, EigenMethod m=EIGEN_SIMPLE) const; /** // Inverse of a real symmetric matrix. If the matrix is not symmetric, // its symmetrized version will be used internally. Use this function // (it has better numerical precision than inv()) if you know that // your matrix is symmetric. */ Matrix symInv() const; /** Inverse of a square tridiagonal matrix */ Matrix triDiagInv() const; /** // Inverse of a general matrix (which must be square and non-singular) */ Matrix inv() const; /** // Left pseudoinverse of a real, full rank matrix. This // pseudoinverse can be calculated when the number of columns // does not exceed the number of rows. */ Matrix leftInv() const; /** // Right pseudoinverse of a real, full rank matrix. This // pseudoinverse can be calculated when the number of rows // does not exceed the number of columns. */ Matrix rightInv() const; /** // Eigenvalues and eigenvectors of a real tridiagonal symmetric // matrix. The "eigenvectors" pointer can be NULL if only eigenvalues // are needed. If it is not NULL, the _rows_ of that matrix will be // set to the computed eigenvectors. // // The "m" parameter specifies the LAPACK method to use for // calculating the eigendecomposition. */ void tdSymEigen(Numeric* eigenvalues, unsigned lenEigenvalues, Matrix* eigenvectors=0, EigenMethod m=EIGEN_RRR) const; /** // Eigenvalues and eigenvectors of a real symmetric matrix. // The "eigenvectors" pointer can be NULL if only eigenvalues // are needed. If it is not NULL, the _rows_ of that matrix // will be set to the computed eigenvectors. // // The "m" parameter specifies the LAPACK method to use for // calculating the eigendecomposition. */ void symEigen(Numeric* eigenvalues, unsigned lenEigenvalues, Matrix* eigenvectors=0, EigenMethod m=EIGEN_SIMPLE) const; /** // Eigenvalues and eigenvectors of a real general matrix. // Either "rightEigenvectors" or "leftEigenvectors" or both // can be NULL if they are not needed. // // If an eigenvalue is real, the corresponding row of the // output eigenvector matrix is set to the computed eigenvector. // If an eigenvalue is complex, all such eigenvalues come in // complex conjugate pairs. Corresponding eigenvectors make up // such a pair as well. The matrix row which corresponds to the // first eigenvalue in the pair is set to the real part of the pair // and the matrix row which corresponds to the second eigenvalue // in the pair is set to the imaginary part of the pair. For the // right eigenvectors, the first complex eigenvector in the pair // needs to be constructed by adding the imaginary part. The second // eigenvector has the imaginary part subtracted. For the left // eigenvectors the order is reversed: the first eigenvector needs // to be constructed by subtracting the imaginary part, and the // second by adding it. */ void genEigen(typename GeneralizedComplex::type *eigenvalues, unsigned lenEigenvalues, Matrix* rightEigenvectors = 0, Matrix* leftEigenvectors = 0) const; /** // Singular value decomposition of a matrix. The length of the // singular values buffer must be at least min(nrows, ncolumns). // // On exit, matrices U and V will contain singular vectors of // the current matrix, row-by-row. The decomposition is then // // *this = U^T * diag(singularValues, nrows, ncolumns) * V */ void svd(Numeric* singularValues, unsigned lenSingularValues, Matrix* U, Matrix* V, SvdMethod m = SVD_D_AND_C) const; /** // Calculate entropy-based effective rank. It is assumed that // the matrix is symmetric positive semidefinite. // // The "tol" parameter specifies the eigenvalue cutoff: // all eigenvalues equal to tol*maxEigenvalue or smaller // will be ignored. This parameter must be non-negative. // // The "m" parameter specifies the LAPACK method to use for // calculating eigenvalues. // // If "eigenSum" pointer is not NULL, *eigenSum will be filled // on exit with the sum of eigenvalues above the cutoff divided // by the largest eigenvalue. */ double symPSDefEffectiveRank(double tol = 0.0, EigenMethod m = EIGEN_D_AND_C, double* eigenSum = 0) const; /** // Calculate a function of this matrix. It is assumed that the // matrix is symmetric and real. The calculation is performed by // calculating the eignenbasis, evaluating the function of the // eigenvalues, and then transforming back to the original basis. */ template Matrix symFcn(const Functor& fcn) const; /** Power of a matrix. Matrix must be square. */ Matrix pow(unsigned degree) const; /** // The following function derives a correlation matrix // out of a covariance matrix. In the returned matrix // all diagonal elements will be set to 1. */ Matrix covarToCorr() const; /** // The following function derives a covariance matrix // out of a correlation matrix and of the covariances. */ Matrix corrToCovar(const double* variances, unsigned lenVariances) const; //@{ /** Method related to "geners" I/O */ inline gs::ClassId classId() const {return gs::ClassId(*this);} bool write(std::ostream& of) const; //@} static const char* classname(); static inline unsigned version() {return 1;} static void restore(const gs::ClassId& id, std::istream& in, Matrix* m); private: void initialize(unsigned nrows, unsigned ncols); void calcDiagonal(); void makeCoarse(unsigned n, unsigned m, Matrix* result) const; void invertDiagonal(Matrix* result) const; Numeric local_[Len]; Numeric* data_; unsigned nrows_; unsigned ncols_; unsigned len_; bool isDiagonal_; bool diagonalityKnown_; #ifdef SWIG public: inline std::vector mainDiagonal() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::mainDiagonal: uninitialized matrix"); const unsigned len = std::min(nrows_, ncols_); std::vector d(len); for (unsigned i=0; isetData(data, dataLen);} inline void setRow2(unsigned row, const Numeric* data, unsigned dataLen) {this->setRow(row, data, dataLen);} inline void setColumn2(unsigned col, const Numeric* data, unsigned dataLen) {this->setColumn(col, data, dataLen);} inline Matrix timesVector2(const Numeric* data, unsigned dataLen) const {return timesVector(data, dataLen);} inline Matrix rowMultiply2(const Numeric* data, unsigned dataLen) const {return rowMultiply(data, dataLen);} inline Numeric bilinear2(const Numeric* data, unsigned dataLen) const {return bilinear(data, dataLen);} inline std::vector solveLinearSystem2( const Numeric* data, unsigned dataLen) const { std::vector result(dataLen); if (!solveLinearSystem(data, dataLen, dataLen ? &result[0] : 0)) result.clear(); return result; } inline std::pair symPSDefEffectiveRank2( double tol = 0.0, EigenMethod m = EIGEN_D_AND_C) const { double eigenSum = 0; double erunk = this->symPSDefEffectiveRank(tol, m, &eigenSum); return std::pair(erunk, eigenSum); } inline std::vector symEigenValues() const { std::vector result(nrows_); symEigen(&result[0], nrows_, 0); return result; } inline Matrix rmul2(Numeric r) const {return *this * r;} // Perhaps, some day SWIG will be able to wrap the following. Currently, // its type deduction for GeneralizedComplex::type fails. // // inline std::vector::type> // genEigenValues() const // { // std::vector::type> r(nrows_); // genEigen(&r[0], nrows_); // return r; // } #endif // SWIG }; /** Utility for making square diagonal matrices from the given array */ template Matrix diag(const Numeric* data, unsigned dataLen); /** // Utility for making rectangular diagonal matrices from the given array. // The length of the array must be at least min(nrows, ncols). */ template Matrix diag(const Numeric* data, unsigned nrows, unsigned ncols); } #include template std::ostream& operator<<(std::ostream& os, const npstat::Matrix& m); template inline npstat::Matrix operator*(const Numeric& l, const npstat::Matrix& r) { return r*l; } #include "npstat/nm/Matrix.icc" #endif // NPSTAT_MATRIX_HH_