Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 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 | // 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 |
Austin Schuh | 189376f | 2018-12-20 22:11:15 +1100 | [diff] [blame] | 12 | |
| 13 | namespace Eigen { |
| 14 | |
Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 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 | */ |
Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 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 | { |
Austin Schuh | 189376f | 2018-12-20 22:11:15 +1100 | [diff] [blame] | 75 | using std::abs; |
Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 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)); |
Austin Schuh | c55b017 | 2022-02-20 17:52:35 -0800 | [diff] [blame^] | 107 | } |
| 108 | for (int i = 0; i < n; ++i) |
| 109 | { |
Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 110 | Dc(i) = std::sqrt(Dc(i)); |
| 111 | } |
| 112 | // Save the scaling factors |
| 113 | for (int i = 0; i < m; ++i) |
| 114 | { |
| 115 | m_left(i) /= Dr(i); |
Austin Schuh | c55b017 | 2022-02-20 17:52:35 -0800 | [diff] [blame^] | 116 | } |
| 117 | for (int i = 0; i < n; ++i) |
| 118 | { |
Brian Silverman | 72890c2 | 2015-09-19 14:37:37 -0400 | [diff] [blame] | 119 | m_right(i) /= Dc(i); |
| 120 | } |
| 121 | // Scale the rows and the columns of the matrix |
| 122 | DrRes.setZero(); DcRes.setZero(); |
| 123 | for (int k=0; k<m_matrix.outerSize(); ++k) |
| 124 | { |
| 125 | for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) |
| 126 | { |
| 127 | it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) ); |
| 128 | // Accumulate the norms of the row and column vectors |
| 129 | if ( DrRes(it.row()) < abs(it.value()) ) |
| 130 | DrRes(it.row()) = abs(it.value()); |
| 131 | |
| 132 | if ( DcRes(it.col()) < abs(it.value()) ) |
| 133 | DcRes(it.col()) = abs(it.value()); |
| 134 | } |
| 135 | } |
| 136 | DrRes.array() = (1-DrRes.array()).abs(); |
| 137 | EpsRow = DrRes.maxCoeff(); |
| 138 | DcRes.array() = (1-DcRes.array()).abs(); |
| 139 | EpsCol = DcRes.maxCoeff(); |
| 140 | its++; |
| 141 | }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) ); |
| 142 | m_isInitialized = true; |
| 143 | } |
| 144 | /** Compute the left and right vectors to scale the vectors |
| 145 | * the input matrix is scaled with the computed vectors at output |
| 146 | * |
| 147 | * \sa compute() |
| 148 | */ |
| 149 | void computeRef (MatrixType& mat) |
| 150 | { |
| 151 | compute (mat); |
| 152 | mat = m_matrix; |
| 153 | } |
| 154 | /** Get the vector to scale the rows of the matrix |
| 155 | */ |
| 156 | VectorXd& LeftScaling() |
| 157 | { |
| 158 | return m_left; |
| 159 | } |
| 160 | |
| 161 | /** Get the vector to scale the columns of the matrix |
| 162 | */ |
| 163 | VectorXd& RightScaling() |
| 164 | { |
| 165 | return m_right; |
| 166 | } |
| 167 | |
| 168 | /** Set the tolerance for the convergence of the iterative scaling algorithm |
| 169 | */ |
| 170 | void setTolerance(double tol) |
| 171 | { |
| 172 | m_tol = tol; |
| 173 | } |
| 174 | |
| 175 | protected: |
| 176 | |
| 177 | void init() |
| 178 | { |
| 179 | m_tol = 1e-10; |
| 180 | m_maxits = 5; |
| 181 | m_isInitialized = false; |
| 182 | } |
| 183 | |
| 184 | MatrixType m_matrix; |
| 185 | mutable ComputationInfo m_info; |
| 186 | bool m_isInitialized; |
| 187 | VectorXd m_left; // Left scaling vector |
| 188 | VectorXd m_right; // m_right scaling vector |
| 189 | double m_tol; |
| 190 | int m_maxits; // Maximum number of iterations allowed |
| 191 | }; |
| 192 | } |
| 193 | #endif |