Scaling.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_ITERSCALING_H
#define EIGEN_ITERSCALING_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

template <typename MatrixType_>
class IterScaling {
 public:
  typedef MatrixType_ MatrixType;
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::Index Index;

 public:
  IterScaling() { init(); }

  IterScaling(const MatrixType& matrix) {
    init();
    compute(matrix);
  }

  ~IterScaling() {}

  void compute(const MatrixType& mat) {
    using std::abs;
    int m = mat.rows();
    int n = mat.cols();
    eigen_assert((m > 0 && m == n) && "Please give a non - empty matrix");
    m_left.resize(m);
    m_right.resize(n);
    m_left.setOnes();
    m_right.setOnes();
    m_matrix = mat;
    VectorXd Dr, Dc, DrRes, DcRes;  // Temporary Left and right scaling vectors
    Dr.resize(m);
    Dc.resize(n);
    DrRes.resize(m);
    DcRes.resize(n);
    double EpsRow = 1.0, EpsCol = 1.0;
    int its = 0;
    do {  // Iterate until the infinite norm of each row and column is approximately 1
      // Get the maximum value in each row and column
      Dr.setZero();
      Dc.setZero();
      for (int k = 0; k < m_matrix.outerSize(); ++k) {
        for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) {
          if (Dr(it.row()) < abs(it.value())) Dr(it.row()) = abs(it.value());

          if (Dc(it.col()) < abs(it.value())) Dc(it.col()) = abs(it.value());
        }
      }
      for (int i = 0; i < m; ++i) {
        Dr(i) = std::sqrt(Dr(i));
      }
      for (int i = 0; i < n; ++i) {
        Dc(i) = std::sqrt(Dc(i));
      }
      // Save the scaling factors
      for (int i = 0; i < m; ++i) {
        m_left(i) /= Dr(i);
      }
      for (int i = 0; i < n; ++i) {
        m_right(i) /= Dc(i);
      }
      // Scale the rows and the columns of the matrix
      DrRes.setZero();
      DcRes.setZero();
      for (int k = 0; k < m_matrix.outerSize(); ++k) {
        for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) {
          it.valueRef() = it.value() / (Dr(it.row()) * Dc(it.col()));
          // Accumulate the norms of the row and column vectors
          if (DrRes(it.row()) < abs(it.value())) DrRes(it.row()) = abs(it.value());

          if (DcRes(it.col()) < abs(it.value())) DcRes(it.col()) = abs(it.value());
        }
      }
      DrRes.array() = (1 - DrRes.array()).abs();
      EpsRow = DrRes.maxCoeff();
      DcRes.array() = (1 - DcRes.array()).abs();
      EpsCol = DcRes.maxCoeff();
      its++;
    } while ((EpsRow > m_tol || EpsCol > m_tol) && (its < m_maxits));
    m_isInitialized = true;
  }
  void computeRef(MatrixType& mat) {
    compute(mat);
    mat = m_matrix;
  }
  VectorXd& LeftScaling() { return m_left; }

  VectorXd& RightScaling() { return m_right; }

  void setTolerance(double tol) { m_tol = tol; }

 protected:
  void init() {
    m_tol = 1e-10;
    m_maxits = 5;
    m_isInitialized = false;
  }

  MatrixType m_matrix;
  mutable ComputationInfo m_info;
  bool m_isInitialized;
  VectorXd m_left;   // Left scaling vector
  VectorXd m_right;  // m_right scaling vector
  double m_tol;
  int m_maxits;  // Maximum number of iterations allowed
};
}  // namespace Eigen
#endif