Home | History | Annotate | Download | only in IterativeSolvers 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam (at) inria.fr
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_ITERSCALING_H
     11 #define EIGEN_ITERSCALING_H
     12 
     13 namespace Eigen {
     14 
     15 /**
     16   * \ingroup IterativeSolvers_Module
     17   * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
     18   *
     19   * This class can be used as a preprocessing tool to accelerate the convergence of iterative methods
     20   *
     21   * This feature is  useful to limit the pivoting amount during LU/ILU factorization
     22   * The  scaling strategy as presented here preserves the symmetry of the problem
     23   * NOTE It is assumed that the matrix does not have empty row or column,
     24   *
     25   * Example with key steps
     26   * \code
     27   * VectorXd x(n), b(n);
     28   * SparseMatrix<double> A;
     29   * // fill A and b;
     30   * IterScaling<SparseMatrix<double> > scal;
     31   * // Compute the left and right scaling vectors. The matrix is equilibrated at output
     32   * scal.computeRef(A);
     33   * // Scale the right hand side
     34   * b = scal.LeftScaling().cwiseProduct(b);
     35   * // Now, solve the equilibrated linear system with any available solver
     36   *
     37   * // Scale back the computed solution
     38   * x = scal.RightScaling().cwiseProduct(x);
     39   * \endcode
     40   *
     41   * \tparam _MatrixType the type of the matrix. It should be a real square sparsematrix
     42   *
     43   * References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552
     44   *
     45   * \sa \ref IncompleteLUT
     46   */
     47 template<typename _MatrixType>
     48 class IterScaling
     49 {
     50   public:
     51     typedef _MatrixType MatrixType;
     52     typedef typename MatrixType::Scalar Scalar;
     53     typedef typename MatrixType::Index Index;
     54 
     55   public:
     56     IterScaling() { init(); }
     57 
     58     IterScaling(const MatrixType& matrix)
     59     {
     60       init();
     61       compute(matrix);
     62     }
     63 
     64     ~IterScaling() { }
     65 
     66     /**
     67      * Compute the left and right diagonal matrices to scale the input matrix @p mat
     68      *
     69      * FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal.
     70      *
     71      * \sa LeftScaling() RightScaling()
     72      */
     73     void compute (const MatrixType& mat)
     74     {
     75       using std::abs;
     76       int m = mat.rows();
     77       int n = mat.cols();
     78       eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");
     79       m_left.resize(m);
     80       m_right.resize(n);
     81       m_left.setOnes();
     82       m_right.setOnes();
     83       m_matrix = mat;
     84       VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
     85       Dr.resize(m); Dc.resize(n);
     86       DrRes.resize(m); DcRes.resize(n);
     87       double EpsRow = 1.0, EpsCol = 1.0;
     88       int its = 0;
     89       do
     90       { // Iterate until the infinite norm of each row and column is approximately 1
     91         // Get the maximum value in each row and column
     92         Dr.setZero(); Dc.setZero();
     93         for (int k=0; k<m_matrix.outerSize(); ++k)
     94         {
     95           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
     96           {
     97             if ( Dr(it.row()) < abs(it.value()) )
     98               Dr(it.row()) = abs(it.value());
     99 
    100             if ( Dc(it.col()) < abs(it.value()) )
    101               Dc(it.col()) = abs(it.value());
    102           }
    103         }
    104         for (int i = 0; i < m; ++i)
    105         {
    106           Dr(i) = std::sqrt(Dr(i));
    107           Dc(i) = std::sqrt(Dc(i));
    108         }
    109         // Save the scaling factors
    110         for (int i = 0; i < m; ++i)
    111         {
    112           m_left(i) /= Dr(i);
    113           m_right(i) /= Dc(i);
    114         }
    115         // Scale the rows and the columns of the matrix
    116         DrRes.setZero(); DcRes.setZero();
    117         for (int k=0; k<m_matrix.outerSize(); ++k)
    118         {
    119           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
    120           {
    121             it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
    122             // Accumulate the norms of the row and column vectors
    123             if ( DrRes(it.row()) < abs(it.value()) )
    124               DrRes(it.row()) = abs(it.value());
    125 
    126             if ( DcRes(it.col()) < abs(it.value()) )
    127               DcRes(it.col()) = abs(it.value());
    128           }
    129         }
    130         DrRes.array() = (1-DrRes.array()).abs();
    131         EpsRow = DrRes.maxCoeff();
    132         DcRes.array() = (1-DcRes.array()).abs();
    133         EpsCol = DcRes.maxCoeff();
    134         its++;
    135       }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
    136       m_isInitialized = true;
    137     }
    138     /** Compute the left and right vectors to scale the vectors
    139      * the input matrix is scaled with the computed vectors at output
    140      *
    141      * \sa compute()
    142      */
    143     void computeRef (MatrixType& mat)
    144     {
    145       compute (mat);
    146       mat = m_matrix;
    147     }
    148     /** Get the vector to scale the rows of the matrix
    149      */
    150     VectorXd& LeftScaling()
    151     {
    152       return m_left;
    153     }
    154 
    155     /** Get the vector to scale the columns of the matrix
    156      */
    157     VectorXd& RightScaling()
    158     {
    159       return m_right;
    160     }
    161 
    162     /** Set the tolerance for the convergence of the iterative scaling algorithm
    163      */
    164     void setTolerance(double tol)
    165     {
    166       m_tol = tol;
    167     }
    168 
    169   protected:
    170 
    171     void init()
    172     {
    173       m_tol = 1e-10;
    174       m_maxits = 5;
    175       m_isInitialized = false;
    176     }
    177 
    178     MatrixType m_matrix;
    179     mutable ComputationInfo m_info;
    180     bool m_isInitialized;
    181     VectorXd m_left; // Left scaling vector
    182     VectorXd m_right; // m_right scaling vector
    183     double m_tol;
    184     int m_maxits; // Maximum number of iterations allowed
    185 };
    186 }
    187 #endif
    188