Scaling.h
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@inria.fr
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_SCALING_H
26 #define EIGEN_SCALING_H
27 
59 using std::abs;
60 using namespace Eigen;
61 template<typename _MatrixType>
62 class Scaling
63 {
64  public:
65  typedef _MatrixType MatrixType;
66  typedef typename MatrixType::Scalar Scalar;
67  typedef typename MatrixType::Index Index;
68 
69  public:
70  Scaling() { init(); }
71 
72  Scaling(const MatrixType& matrix)
73  {
74  init();
75  compute(matrix);
76  }
77 
78  ~Scaling() { }
79 
87  void compute (const MatrixType& mat)
88  {
89  int m = mat.rows();
90  int n = mat.cols();
91  assert((m>0 && m == n) && "Please give a non - empty matrix");
92  m_left.resize(m);
93  m_right.resize(n);
94  m_left.setOnes();
95  m_right.setOnes();
96  m_matrix = mat;
97  VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
98  Dr.resize(m); Dc.resize(n);
99  DrRes.resize(m); DcRes.resize(n);
100  double EpsRow = 1.0, EpsCol = 1.0;
101  int its = 0;
102  do
103  { // Iterate until the infinite norm of each row and column is approximately 1
104  // Get the maximum value in each row and column
105  Dr.setZero(); Dc.setZero();
106  for (int k=0; k<m_matrix.outerSize(); ++k)
107  {
108  for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
109  {
110  if ( Dr(it.row()) < abs(it.value()) )
111  Dr(it.row()) = abs(it.value());
112 
113  if ( Dc(it.col()) < abs(it.value()) )
114  Dc(it.col()) = abs(it.value());
115  }
116  }
117  for (int i = 0; i < m; ++i)
118  {
119  Dr(i) = std::sqrt(Dr(i));
120  Dc(i) = std::sqrt(Dc(i));
121  }
122  // Save the scaling factors
123  for (int i = 0; i < m; ++i)
124  {
125  m_left(i) /= Dr(i);
126  m_right(i) /= Dc(i);
127  }
128  // Scale the rows and the columns of the matrix
129  DrRes.setZero(); DcRes.setZero();
130  for (int k=0; k<m_matrix.outerSize(); ++k)
131  {
132  for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
133  {
134  it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
135  // Accumulate the norms of the row and column vectors
136  if ( DrRes(it.row()) < abs(it.value()) )
137  DrRes(it.row()) = abs(it.value());
138 
139  if ( DcRes(it.col()) < abs(it.value()) )
140  DcRes(it.col()) = abs(it.value());
141  }
142  }
143  DrRes.array() = (1-DrRes.array()).abs();
144  EpsRow = DrRes.maxCoeff();
145  DcRes.array() = (1-DcRes.array()).abs();
146  EpsCol = DcRes.maxCoeff();
147  its++;
148  }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
149  m_isInitialized = true;
150  }
156  void computeRef (MatrixType& mat)
157  {
158  compute (mat);
159  mat = m_matrix;
160  }
163  VectorXd& LeftScaling()
164  {
165  return m_left;
166  }
167 
170  VectorXd& RightScaling()
171  {
172  return m_right;
173  }
174 
177  void setTolerance(double tol)
178  {
179  m_tol = tol;
180  }
181 
182  protected:
183 
184  void init()
185  {
186  m_tol = 1e-10;
187  m_maxits = 5;
188  m_isInitialized = false;
189  }
190 
191  MatrixType m_matrix;
192  mutable ComputationInfo m_info;
193  bool m_isInitialized;
194  VectorXd m_left; // Left scaling vector
195  VectorXd m_right; // m_right scaling vector
196  double m_tol;
197  int m_maxits; // Maximum number of iterations allowed
198 };
199 
200 #endif