src/IterativeSolvers/Scaling.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_ITERSCALING_H
0011 #define EIGEN_ITERSCALING_H
0012
0013 namespace Eigen {
0014
0015 /**
0016   * \ingroup IterativeSolvers_Module
0017   * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
0018   *
0019   * This class can be used as a preprocessing tool to accelerate the convergence of iterative methods
0020   *
0021   * This feature is  useful to limit the pivoting amount during LU/ILU factorization
0022   * The  scaling strategy as presented here preserves the symmetry of the problem
0023   * NOTE It is assumed that the matrix does not have empty row or column,
0024   *
0025   * Example with key steps
0026   * \code
0027   * VectorXd x(n), b(n);
0028   * SparseMatrix<double> A;
0029   * // fill A and b;
0030   * IterScaling<SparseMatrix<double> > scal;
0031   * // Compute the left and right scaling vectors. The matrix is equilibrated at output
0032   * scal.computeRef(A);
0033   * // Scale the right hand side
0034   * b = scal.LeftScaling().cwiseProduct(b);
0035   * // Now, solve the equilibrated linear system with any available solver
0036   *
0037   * // Scale back the computed solution
0038   * x = scal.RightScaling().cwiseProduct(x);
0039   * \endcode
0040   *
0041   * \tparam _MatrixType the type of the matrix. It should be a real square sparsematrix
0042   *
0043   * References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552
0044   *
0045   * \sa \ref IncompleteLUT
0046   */
0047 template<typename _MatrixType>
0048 class IterScaling
0049 {
0050   public:
0051     typedef _MatrixType MatrixType;
0052     typedef typename MatrixType::Scalar Scalar;
0053     typedef typename MatrixType::Index Index;
0054
0055   public:
0056     IterScaling() { init(); }
0057
0058     IterScaling(const MatrixType& matrix)
0059     {
0060       init();
0061       compute(matrix);
0062     }
0063
0064     ~IterScaling() { }
0065
0066     /**
0067      * Compute the left and right diagonal matrices to scale the input matrix @p mat
0068      *
0069      * FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal.
0070      *
0071      * \sa LeftScaling() RightScaling()
0072      */
0073     void compute (const MatrixType& mat)
0074     {
0075       using std::abs;
0076       int m = mat.rows();
0077       int n = mat.cols();
0078       eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");
0079       m_left.resize(m);
0080       m_right.resize(n);
0081       m_left.setOnes();
0082       m_right.setOnes();
0083       m_matrix = mat;
0084       VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
0085       Dr.resize(m); Dc.resize(n);
0086       DrRes.resize(m); DcRes.resize(n);
0087       double EpsRow = 1.0, EpsCol = 1.0;
0088       int its = 0;
0089       do
0090       { // Iterate until the infinite norm of each row and column is approximately 1
0091         // Get the maximum value in each row and column
0092         Dr.setZero(); Dc.setZero();
0093         for (int k=0; k<m_matrix.outerSize(); ++k)
0094         {
0095           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
0096           {
0097             if ( Dr(it.row()) < abs(it.value()) )
0098               Dr(it.row()) = abs(it.value());
0099
0100             if ( Dc(it.col()) < abs(it.value()) )
0101               Dc(it.col()) = abs(it.value());
0102           }
0103         }
0104         for (int i = 0; i < m; ++i)
0105         {
0106           Dr(i) = std::sqrt(Dr(i));
0107         }
0108         for (int i = 0; i < n; ++i)
0109         {
0110           Dc(i) = std::sqrt(Dc(i));
0111         }
0112         // Save the scaling factors
0113         for (int i = 0; i < m; ++i)
0114         {
0115           m_left(i) /= Dr(i);
0116         }
0117         for (int i = 0; i < n; ++i)
0118         {
0119           m_right(i) /= Dc(i);
0120         }
0121         // Scale the rows and the columns of the matrix
0122         DrRes.setZero(); DcRes.setZero();
0123         for (int k=0; k<m_matrix.outerSize(); ++k)
0124         {
0125           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
0126           {
0127             it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
0128             // Accumulate the norms of the row and column vectors
0129             if ( DrRes(it.row()) < abs(it.value()) )
0130               DrRes(it.row()) = abs(it.value());
0131
0132             if ( DcRes(it.col()) < abs(it.value()) )
0133               DcRes(it.col()) = abs(it.value());
0134           }
0135         }
0136         DrRes.array() = (1-DrRes.array()).abs();
0137         EpsRow = DrRes.maxCoeff();
0138         DcRes.array() = (1-DcRes.array()).abs();
0139         EpsCol = DcRes.maxCoeff();
0140         its++;
0141       }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
0142       m_isInitialized = true;
0143     }
0144     /** Compute the left and right vectors to scale the vectors
0145      * the input matrix is scaled with the computed vectors at output
0146      *
0147      * \sa compute()
0148      */
0149     void computeRef (MatrixType& mat)
0150     {
0151       compute (mat);
0152       mat = m_matrix;
0153     }
0154     /** Get the vector to scale the rows of the matrix
0155      */
0156     VectorXd& LeftScaling()
0157     {
0158       return m_left;
0159     }
0160
0161     /** Get the vector to scale the columns of the matrix
0162      */
0163     VectorXd& RightScaling()
0164     {
0165       return m_right;
0166     }
0167
0168     /** Set the tolerance for the convergence of the iterative scaling algorithm
0169      */
0170     void setTolerance(double tol)
0171     {
0172       m_tol = tol;
0173     }
0174
0175   protected:
0176
0177     void init()
0178     {
0179       m_tol = 1e-10;
0180       m_maxits = 5;
0181       m_isInitialized = false;
0182     }
0183
0184     MatrixType m_matrix;
0185     mutable ComputationInfo m_info;
0186     bool m_isInitialized;
0187     VectorXd m_left; // Left scaling vector
0188     VectorXd m_right; // m_right scaling vector
0189     double m_tol;
0190     int m_maxits; // Maximum number of iterations allowed
0191 };
0192 }
0193 #endif