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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_SCALING_H
00026 #define EIGEN_SCALING_H
00027 
00059 using std::abs; 
00060 using namespace Eigen;
00061 template<typename _MatrixType>
00062 class Scaling
00063 {
00064   public:
00065     typedef _MatrixType MatrixType; 
00066     typedef typename MatrixType::Scalar Scalar;
00067     typedef typename MatrixType::Index Index;
00068     
00069   public:
00070     Scaling() { init(); }
00071     
00072     Scaling(const MatrixType& matrix)
00073     {
00074       init();
00075       compute(matrix);
00076     }
00077     
00078     ~Scaling() { }
00079     
00087     void compute (const MatrixType& mat)
00088     {
00089       int m = mat.rows(); 
00090       int n = mat.cols();
00091       assert((m>0 && m == n) && "Please give a non - empty matrix");
00092       m_left.resize(m); 
00093       m_right.resize(n);
00094       m_left.setOnes();
00095       m_right.setOnes();
00096       m_matrix = mat;
00097       VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
00098       Dr.resize(m); Dc.resize(n);
00099       DrRes.resize(m); DcRes.resize(n);
00100       double EpsRow = 1.0, EpsCol = 1.0;
00101       int its = 0; 
00102       do
00103       { // Iterate until the infinite norm of each row and column is approximately 1
00104         // Get the maximum value in each row and column
00105         Dr.setZero(); Dc.setZero();
00106         for (int k=0; k<m_matrix.outerSize(); ++k)
00107         {
00108           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00109           {
00110             if ( Dr(it.row()) < abs(it.value()) )
00111               Dr(it.row()) = abs(it.value());
00112             
00113             if ( Dc(it.col()) < abs(it.value()) )
00114               Dc(it.col()) = abs(it.value());
00115           }
00116         }
00117         for (int i = 0; i < m; ++i) 
00118         {
00119           Dr(i) = std::sqrt(Dr(i));
00120           Dc(i) = std::sqrt(Dc(i));
00121         }
00122         // Save the scaling factors 
00123         for (int i = 0; i < m; ++i) 
00124         {
00125           m_left(i) /= Dr(i);
00126           m_right(i) /= Dc(i);
00127         }
00128         // Scale the rows and the columns of the matrix
00129         DrRes.setZero(); DcRes.setZero(); 
00130         for (int k=0; k<m_matrix.outerSize(); ++k)
00131         {
00132           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00133           {
00134             it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
00135             // Accumulate the norms of the row and column vectors   
00136             if ( DrRes(it.row()) < abs(it.value()) )
00137               DrRes(it.row()) = abs(it.value());
00138             
00139             if ( DcRes(it.col()) < abs(it.value()) )
00140               DcRes(it.col()) = abs(it.value());
00141           }
00142         }  
00143         DrRes.array() = (1-DrRes.array()).abs();
00144         EpsRow = DrRes.maxCoeff();
00145         DcRes.array() = (1-DcRes.array()).abs();
00146         EpsCol = DcRes.maxCoeff();
00147         its++;
00148       }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
00149       m_isInitialized = true;
00150     }
00156     void computeRef (MatrixType& mat)
00157     {
00158       compute (mat);
00159       mat = m_matrix;
00160     }
00163     VectorXd& LeftScaling()
00164     {
00165       return m_left;
00166     }
00167     
00170     VectorXd& RightScaling()
00171     {
00172       return m_right;
00173     }
00174     
00177     void setTolerance(double tol)
00178     {
00179       m_tol = tol; 
00180     }
00181       
00182   protected:
00183     
00184     void init()
00185     {
00186       m_tol = 1e-10;
00187       m_maxits = 5;
00188       m_isInitialized = false;
00189     }
00190     
00191     MatrixType m_matrix;
00192     mutable ComputationInfo m_info; 
00193     bool m_isInitialized; 
00194     VectorXd m_left; // Left scaling vector
00195     VectorXd m_right; // m_right scaling vector
00196     double m_tol; 
00197     int m_maxits; // Maximum number of iterations allowed
00198 };
00199 
00200 #endif