Scaling.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr
00005 //
00006 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 #ifndef EIGEN_ITERSCALING_H
00011 #define EIGEN_ITERSCALING_H
00012 
00013 namespace Eigen {
00014 
00047 template<typename _MatrixType>
00048 class IterScaling
00049 {
00050   public:
00051     typedef _MatrixType MatrixType; 
00052     typedef typename MatrixType::Scalar Scalar;
00053     typedef typename MatrixType::Index Index;
00054     
00055   public:
00056     IterScaling() { init(); }
00057     
00058     IterScaling(const MatrixType& matrix)
00059     {
00060       init();
00061       compute(matrix);
00062     }
00063     
00064     ~IterScaling() { }
00065     
00073     void compute (const MatrixType& mat)
00074     {
00075       using std::abs;
00076       int m = mat.rows(); 
00077       int n = mat.cols();
00078       eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");
00079       m_left.resize(m); 
00080       m_right.resize(n);
00081       m_left.setOnes();
00082       m_right.setOnes();
00083       m_matrix = mat;
00084       VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
00085       Dr.resize(m); Dc.resize(n);
00086       DrRes.resize(m); DcRes.resize(n);
00087       double EpsRow = 1.0, EpsCol = 1.0;
00088       int its = 0; 
00089       do
00090       { // Iterate until the infinite norm of each row and column is approximately 1
00091         // Get the maximum value in each row and column
00092         Dr.setZero(); Dc.setZero();
00093         for (int k=0; k<m_matrix.outerSize(); ++k)
00094         {
00095           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00096           {
00097             if ( Dr(it.row()) < abs(it.value()) )
00098               Dr(it.row()) = abs(it.value());
00099             
00100             if ( Dc(it.col()) < abs(it.value()) )
00101               Dc(it.col()) = abs(it.value());
00102           }
00103         }
00104         for (int i = 0; i < m; ++i) 
00105         {
00106           Dr(i) = std::sqrt(Dr(i));
00107           Dc(i) = std::sqrt(Dc(i));
00108         }
00109         // Save the scaling factors 
00110         for (int i = 0; i < m; ++i) 
00111         {
00112           m_left(i) /= Dr(i);
00113           m_right(i) /= Dc(i);
00114         }
00115         // Scale the rows and the columns of the matrix
00116         DrRes.setZero(); DcRes.setZero(); 
00117         for (int k=0; k<m_matrix.outerSize(); ++k)
00118         {
00119           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00120           {
00121             it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
00122             // Accumulate the norms of the row and column vectors   
00123             if ( DrRes(it.row()) < abs(it.value()) )
00124               DrRes(it.row()) = abs(it.value());
00125             
00126             if ( DcRes(it.col()) < abs(it.value()) )
00127               DcRes(it.col()) = abs(it.value());
00128           }
00129         }  
00130         DrRes.array() = (1-DrRes.array()).abs();
00131         EpsRow = DrRes.maxCoeff();
00132         DcRes.array() = (1-DcRes.array()).abs();
00133         EpsCol = DcRes.maxCoeff();
00134         its++;
00135       }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
00136       m_isInitialized = true;
00137     }
00143     void computeRef (MatrixType& mat)
00144     {
00145       compute (mat);
00146       mat = m_matrix;
00147     }
00150     VectorXd& LeftScaling()
00151     {
00152       return m_left;
00153     }
00154     
00157     VectorXd& RightScaling()
00158     {
00159       return m_right;
00160     }
00161     
00164     void setTolerance(double tol)
00165     {
00166       m_tol = tol; 
00167     }
00168       
00169   protected:
00170     
00171     void init()
00172     {
00173       m_tol = 1e-10;
00174       m_maxits = 5;
00175       m_isInitialized = false;
00176     }
00177     
00178     MatrixType m_matrix;
00179     mutable ComputationInfo m_info; 
00180     bool m_isInitialized; 
00181     VectorXd m_left; // Left scaling vector
00182     VectorXd m_right; // m_right scaling vector
00183     double m_tol; 
00184     int m_maxits; // Maximum number of iterations allowed
00185 };
00186 }
00187 #endif
 All Classes Functions Variables Typedefs Enumerator