From 5440d164b2288e3701c5e6cc247631439f3487aa Mon Sep 17 00:00:00 2001 From: dehnert Date: Mon, 11 Sep 2017 13:27:48 +0200 Subject: [PATCH] started on Walker-Chae algorithm --- .../modules/NativeEquationSolverSettings.cpp | 4 +- .../modules/NativeEquationSolverSettings.h | 2 +- .../solver/NativeLinearEquationSolver.cpp | 308 +++++++++++++----- src/storm/solver/NativeLinearEquationSolver.h | 11 +- src/storm/storage/SparseMatrix.cpp | 31 +- src/storm/storage/SparseMatrix.h | 11 + src/storm/utility/vector.h | 17 + 7 files changed, 303 insertions(+), 81 deletions(-) diff --git a/src/storm/settings/modules/NativeEquationSolverSettings.cpp b/src/storm/settings/modules/NativeEquationSolverSettings.cpp index 4d5c8ce59..b08ac40c0 100644 --- a/src/storm/settings/modules/NativeEquationSolverSettings.cpp +++ b/src/storm/settings/modules/NativeEquationSolverSettings.cpp @@ -24,7 +24,7 @@ namespace storm { const std::string NativeEquationSolverSettings::absoluteOptionName = "absolute"; NativeEquationSolverSettings::NativeEquationSolverSettings() : ModuleSettings(moduleName) { - std::vector methods = { "jacobi", "gaussseidel", "sor" }; + std::vector methods = { "jacobi", "gaussseidel", "sor", "walkerchae" }; this->addOption(storm::settings::OptionBuilder(moduleName, techniqueOptionName, true, "The method to be used for solving linear equation systems with the native engine.").addArgument(storm::settings::ArgumentBuilder::createStringArgument("name", "The name of the method to use.").addValidatorString(ArgumentValidatorFactory::createMultipleChoiceValidator(methods)).setDefaultValueString("jacobi").build()).build()); this->addOption(storm::settings::OptionBuilder(moduleName, maximalIterationsOptionName, false, "The maximal number of iterations to perform before iterative solving is aborted.").setShortName(maximalIterationsOptionShortName).addArgument(storm::settings::ArgumentBuilder::createUnsignedIntegerArgument("count", "The maximal iteration count.").setDefaultValueUnsignedInteger(20000).build()).build()); @@ -52,6 +52,8 @@ namespace storm { return NativeEquationSolverSettings::LinearEquationMethod::GaussSeidel; } else if (linearEquationSystemTechniqueAsString == "sor") { return NativeEquationSolverSettings::LinearEquationMethod::SOR; + } else if (linearEquationSystemTechniqueAsString == "walkerchae") { + return NativeEquationSolverSettings::LinearEquationMethod::WalkerChae; } STORM_LOG_THROW(false, storm::exceptions::IllegalArgumentValueException, "Unknown solution technique '" << linearEquationSystemTechniqueAsString << "' selected."); } diff --git a/src/storm/settings/modules/NativeEquationSolverSettings.h b/src/storm/settings/modules/NativeEquationSolverSettings.h index 2699d3f94..f05ef6b25 100644 --- a/src/storm/settings/modules/NativeEquationSolverSettings.h +++ b/src/storm/settings/modules/NativeEquationSolverSettings.h @@ -13,7 +13,7 @@ namespace storm { class NativeEquationSolverSettings : public ModuleSettings { public: // An enumeration of all available methods for solving linear equations. - enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR }; + enum class LinearEquationMethod { Jacobi, GaussSeidel, SOR, WalkerChae }; // An enumeration of all available convergence criteria. enum class ConvergenceCriterion { Absolute, Relative }; diff --git a/src/storm/solver/NativeLinearEquationSolver.cpp b/src/storm/solver/NativeLinearEquationSolver.cpp index 170edc03c..b3c1d4f37 100644 --- a/src/storm/solver/NativeLinearEquationSolver.cpp +++ b/src/storm/solver/NativeLinearEquationSolver.cpp @@ -5,6 +5,7 @@ #include "storm/settings/SettingsManager.h" #include "storm/settings/modules/NativeEquationSolverSettings.h" +#include "storm/utility/constants.h" #include "storm/utility/vector.h" #include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/InvalidSettingsException.h" @@ -23,6 +24,8 @@ namespace storm { method = SolutionMethod::Jacobi; } else if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::SOR) { method = SolutionMethod::SOR; + } else if (methodAsSetting == storm::settings::modules::NativeEquationSolverSettings::LinearEquationMethod::WalkerChae) { + method = SolutionMethod::WalkerChae; } else { STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "The selected solution technique is invalid for this solver."); } @@ -107,101 +110,251 @@ namespace storm { clearCache(); } - template - bool NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { + bool NativeLinearEquationSolver::solveEquationsSOR(std::vector& x, std::vector const& b, ValueType const& omega) const { + STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Gauss-Seidel, SOR omega = " << omega << ")"); + if (!this->cachedRowVector) { this->cachedRowVector = std::make_unique>(getMatrixRowCount()); } - if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::GaussSeidel) { - // Define the omega used for SOR. - ValueType omega = this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one(); - - STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Gauss-Seidel, SOR omega = " << omega << ")"); - - // Set up additional environment variables. - uint_fast64_t iterationCount = 0; - bool converged = false; + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations()) { + A->performSuccessiveOverRelaxationStep(omega, x, b); - while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations()) { - A->performSuccessiveOverRelaxationStep(omega, x, b); - - // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*this->cachedRowVector, x, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); - - // If we did not yet converge, we need to backup the contents of x. - if (!converged) { - *this->cachedRowVector = x; - } - - // Increase iteration count so we can abort if convergence is too slow. - ++iterationCount; - } + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision(*this->cachedRowVector, x, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()) || (this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(x)); - if(!this->isCachingEnabled()) { - clearCache(); + // If we did not yet converge, we need to backup the contents of x. + if (!converged) { + *this->cachedRowVector = x; } - - if (converged) { - STORM_LOG_INFO("Iterative solver converged in " << iterationCount << " iterations."); - } else { - STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterations."); - } - - return converged; + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + if (!this->isCachingEnabled()) { + clearCache(); + } + + if (converged) { + STORM_LOG_INFO("Iterative solver converged in " << iterationCount << " iterations."); } else { - STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Jacobi)"); - - // Get a Jacobi decomposition of the matrix A. - if(!jacobiDecomposition) { - jacobiDecomposition = std::make_unique, std::vector>>(A->getJacobiDecomposition()); - } - storm::storage::SparseMatrix const& jacobiLU = jacobiDecomposition->first; - std::vector const& jacobiD = jacobiDecomposition->second; + STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterations."); + } + + return converged; + } + + template + bool NativeLinearEquationSolver::solveEquationsJacobi(std::vector& x, std::vector const& b) const { + STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (Jacobi)"); + + if (!this->cachedRowVector) { + this->cachedRowVector = std::make_unique>(getMatrixRowCount()); + } + + // Get a Jacobi decomposition of the matrix A. + if (!jacobiDecomposition) { + jacobiDecomposition = std::make_unique, std::vector>>(A->getJacobiDecomposition()); + } + storm::storage::SparseMatrix const& jacobiLU = jacobiDecomposition->first; + std::vector const& jacobiD = jacobiDecomposition->second; + + std::vector* currentX = &x; + std::vector* nextX = this->cachedRowVector.get(); + + // Set up additional environment variables. + uint_fast64_t iterationCount = 0; + bool converged = false; + + while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { + // Compute D^-1 * (b - LU * x) and store result in nextX. + jacobiLU.multiplyWithVector(*currentX, *nextX); + storm::utility::vector::subtractVectors(b, *nextX, *nextX); + storm::utility::vector::multiplyVectorsPointwise(jacobiD, *nextX, *nextX); - std::vector* currentX = &x; - std::vector* nextX = this->cachedRowVector.get(); + // Now check if the process already converged within our precision. + converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()); - // Set up additional environment variables. - uint_fast64_t iterationCount = 0; - bool converged = false; + // Swap the two pointers as a preparation for the next iteration. + std::swap(nextX, currentX); - while (!converged && iterationCount < this->getSettings().getMaximalNumberOfIterations() && !(this->hasCustomTerminationCondition() && this->getTerminationCondition().terminateNow(*currentX))) { - // Compute D^-1 * (b - LU * x) and store result in nextX. - jacobiLU.multiplyWithVector(*currentX, *nextX); - storm::utility::vector::subtractVectors(b, *nextX, *nextX); - storm::utility::vector::multiplyVectorsPointwise(jacobiD, *nextX, *nextX); - - // Now check if the process already converged within our precision. - converged = storm::utility::vector::equalModuloPrecision(*currentX, *nextX, static_cast(this->getSettings().getPrecision()), this->getSettings().getRelativeTerminationCriterion()); - - // Swap the two pointers as a preparation for the next iteration. - std::swap(nextX, currentX); - - // Increase iteration count so we can abort if convergence is too slow. - ++iterationCount; - } - - // If the last iteration did not write to the original x we have to swap the contents, because the - // output has to be written to the input parameter x. - if (currentX == this->cachedRowVector.get()) { - std::swap(x, *currentX); + // Increase iteration count so we can abort if convergence is too slow. + ++iterationCount; + } + + // If the last iteration did not write to the original x we have to swap the contents, because the + // output has to be written to the input parameter x. + if (currentX == this->cachedRowVector.get()) { + std::swap(x, *currentX); + } + + if (!this->isCachingEnabled()) { + clearCache(); + } + + if (converged) { + STORM_LOG_INFO("Iterative solver converged in " << iterationCount << " iterations."); + } else { + STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterations."); + } + + return converged; + } + + template + void NativeLinearEquationSolver::computeWalkerChaeMatrix() const { + storm::storage::BitVector columnsWithNegativeEntries(this->A->getColumnCount()); + ValueType zero = storm::utility::zero(); + for (auto const& e : *this->A) { + if (e.getValue() < zero) { + columnsWithNegativeEntries.set(e.getColumn()); } - - if (!this->isCachingEnabled()) { - clearCache(); + } + std::vector columnsWithNegativeEntriesBefore = columnsWithNegativeEntries.getNumberOfSetBitsBeforeIndices(); + + // We now build an extended equation system matrix that only has non-negative coefficients. + storm::storage::SparseMatrixBuilder builder; + + uint64_t row = 0; + for (; row < this->A->getRowCount(); ++row) { + for (auto const& entry : this->A->getRow(row)) { + if (entry.getValue() < zero) { + builder.addNextValue(row, this->A->getRowCount() + columnsWithNegativeEntriesBefore[entry.getColumn()], -entry.getValue()); + } else { + builder.addNextValue(row, entry.getColumn(), entry.getValue()); + } } + } + ValueType one = storm::utility::one(); + for (auto column : columnsWithNegativeEntries) { + builder.addNextValue(row, column, one); + builder.addNextValue(row, this->A->getRowCount() + columnsWithNegativeEntriesBefore[column], one); + ++row; + } + + walkerChaeMatrix = std::make_unique>(builder.build()); + } + + template + bool NativeLinearEquationSolver::solveEquationsWalkerChae(std::vector& x, std::vector const& b) const { + STORM_LOG_INFO("Solving linear equation system (" << x.size() << " rows) with NativeLinearEquationSolver (WalkerChae)"); + + // (1) Compute an equivalent equation system that has only non-negative coefficients. + if (!walkerChaeMatrix) { + std::cout << *this->A << std::endl; + computeWalkerChaeMatrix(); + std::cout << *walkerChaeMatrix << std::endl; + } - if (converged) { - STORM_LOG_INFO("Iterative solver converged in " << iterationCount << " iterations."); - } else { - STORM_LOG_WARN("Iterative solver did not converge in " << iterationCount << " iterations."); - } + // (2) Enlarge the vectors x and b to account for additional variables. + x.resize(walkerChaeMatrix->getRowCount()); + + if (walkerChaeMatrix->getRowCount() > this->A->getRowCount() && !walkerChaeB) { + walkerChaeB = std::make_unique>(b); + walkerChaeB->resize(x.size()); + } + + // Choose a value for t in the algorithm. + ValueType t = storm::utility::convertNumber(1000); + + // Precompute some data. + std::vector columnSums(x.size()); + for (auto const& e : *walkerChaeMatrix) { + STORM_LOG_ASSERT(e.getValue() >= storm::utility::zero(), "Expecting only non-negative entries in WalkerChae matrix."); + columnSums[e.getColumn()] += e.getValue(); + } + + // Square the error bound, so we can use it to check for convergence. We take the squared error, because we + // do not want to compute the root in the 2-norm computation. + ValueType squaredErrorBound = storm::utility::pow(this->getSettings().getPrecision(), 2); + + // Create a vector that always holds Ax. + std::vector currentAx(x.size()); + walkerChaeMatrix->multiplyWithVector(x, currentAx); + + // Create an auxiliary vector that intermediately stores the result of the Walker-Chae step. + std::vector tmpX(x.size()); + + // Set up references to the x-vectors used in the iteration loop. + std::vector* currentX = &x; + std::vector* nextX = &tmpX; + + // Prepare a function that adds t to its input. + auto addT = [t] (ValueType const& value) { return value + t; }; + + // (3) Perform iterations until convergence. + bool converged = false; + uint64_t iterations = 0; + while (!converged && iterations < this->getSettings().getMaximalNumberOfIterations()) { + // Perform one Walker-Chae step. + A->performWalkerChaeStep(*currentX, columnSums, *walkerChaeB, currentAx, *nextX); - return iterationCount < this->getSettings().getMaximalNumberOfIterations(); + // Compute new Ax. + A->multiplyWithVector(*nextX, currentAx); + + // Check for convergence. + converged = storm::utility::vector::computeSquaredNorm2Difference(currentAx, *walkerChaeB); + + // If the method did not yet converge, we need to update the value of Ax. + if (!converged) { + // TODO: scale matrix diagonal entries with t and add them to *walkerChaeB. + + // Add t to all entries of x. + storm::utility::vector::applyPointwise(x, x, addT); + } + + std::swap(currentX, nextX); + + // Increase iteration count so we can abort if convergence is too slow. + ++iterations; } + + // If the last iteration did not write to the original x we have to swap the contents, because the + // output has to be written to the input parameter x. + if (currentX == &tmpX) { + std::swap(x, *currentX); + } + + if (!this->isCachingEnabled()) { + clearCache(); + } + + if (converged) { + STORM_LOG_INFO("Iterative solver converged in " << iterations << " iterations."); + } else { + STORM_LOG_WARN("Iterative solver did not converge in " << iterations << " iterations."); + } + + // Resize the solution to the right size. + x.resize(this->A->getRowCount()); + + // Finalize solution vector. + storm::utility::vector::applyPointwise(x, x, [&t,iterations] (ValueType const& value) { return value - iterations * t; } ); + + if (!this->isCachingEnabled()) { + clearCache(); + } + return converged; + } + + template + bool NativeLinearEquationSolver::solveEquations(std::vector& x, std::vector const& b) const { + if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR || this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::GaussSeidel) { + return this->solveEquationsSOR(x, b, this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::SOR ? this->getSettings().getOmega() : storm::utility::one()); + } else if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::Jacobi) { + return this->solveEquationsJacobi(x, b); + } else if (this->getSettings().getSolutionMethod() == NativeLinearEquationSolverSettings::SolutionMethod::WalkerChae) { + return this->solveEquationsWalkerChae(x, b); + } + + STORM_LOG_THROW(false, storm::exceptions::InvalidSettingsException, "Unknown solving technique."); + return false; } template @@ -271,6 +424,7 @@ namespace storm { template void NativeLinearEquationSolver::clearCache() const { jacobiDecomposition.reset(); + walkerChaeMatrix.reset(); LinearEquationSolver::clearCache(); } diff --git a/src/storm/solver/NativeLinearEquationSolver.h b/src/storm/solver/NativeLinearEquationSolver.h index 11cb39bd1..ebe180389 100644 --- a/src/storm/solver/NativeLinearEquationSolver.h +++ b/src/storm/solver/NativeLinearEquationSolver.h @@ -12,7 +12,7 @@ namespace storm { class NativeLinearEquationSolverSettings { public: enum class SolutionMethod { - Jacobi, GaussSeidel, SOR + Jacobi, GaussSeidel, SOR, WalkerChae }; NativeLinearEquationSolverSettings(); @@ -65,6 +65,12 @@ namespace storm { virtual uint64_t getMatrixRowCount() const override; virtual uint64_t getMatrixColumnCount() const override; + virtual bool solveEquationsSOR(std::vector& x, std::vector const& b, ValueType const& omega) const; + virtual bool solveEquationsJacobi(std::vector& x, std::vector const& b) const; + virtual bool solveEquationsWalkerChae(std::vector& x, std::vector const& b) const; + + void computeWalkerChaeMatrix() const; + // If the solver takes posession of the matrix, we store the moved matrix in this member, so it gets deleted // when the solver is destructed. std::unique_ptr> localA; @@ -78,6 +84,9 @@ namespace storm { // cached auxiliary data mutable std::unique_ptr, std::vector>> jacobiDecomposition; + + mutable std::unique_ptr> walkerChaeMatrix; + mutable std::unique_ptr> walkerChaeB; }; template diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp index 39c99a8e0..812e1214f 100644 --- a/src/storm/storage/SparseMatrix.cpp +++ b/src/storm/storage/SparseMatrix.cpp @@ -18,6 +18,7 @@ #include "storm/exceptions/InvalidStateException.h" #include "storm/exceptions/NotImplementedException.h" +#include "storm/exceptions/NotSupportedException.h" #include "storm/exceptions/InvalidArgumentException.h" #include "storm/exceptions/OutOfRangeException.h" @@ -1455,10 +1456,38 @@ namespace storm { #ifdef STORM_HAVE_CARL template<> void SparseMatrix::performSuccessiveOverRelaxationStep(Interval, std::vector&, std::vector const&) const { - STORM_LOG_THROW(false, storm::exceptions::NotImplementedException, "This operation is not supported."); + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); } #endif + + template + void SparseMatrix::performWalkerChaeStep(std::vector const& x, std::vector const& columnSums, std::vector const& b, std::vector const& ax, std::vector& result) const { + const_iterator it = this->begin(); + const_iterator ite; + std::vector::const_iterator rowIterator = rowIndications.begin(); + typename std::vector::iterator resultIterator = result.begin(); + typename std::vector::iterator resultIteratorEnd = result.end(); + typename std::vector::const_iterator xIterator = x.begin(); + typename std::vector::const_iterator columnSumsIterator = columnSums.begin(); + + for (; resultIterator != resultIteratorEnd; ++rowIterator, ++resultIterator, ++xIterator, ++columnSumsIterator) { + *resultIterator = storm::utility::zero(); + + for (ite = this->begin() + *(rowIterator + 1); it != ite; ++it) { + *resultIterator += it->getValue() * (b[it->getColumn()] / ax[it->getColumn()]); + } + + *resultIterator *= (*xIterator / *columnSumsIterator); + } + } +#ifdef STORM_HAVE_CARL + template<> + void SparseMatrix::performWalkerChaeStep(std::vector const& x, std::vector const& columnSums, std::vector const& b, std::vector const& ax, std::vector& result) const { + STORM_LOG_THROW(false, storm::exceptions::NotSupportedException, "This operation is not supported."); + } +#endif + template void SparseMatrix::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector const& rowGroupIndices, std::vector const& vector, std::vector const* summand, std::vector& result, std::vector* choices) const { auto elementIt = this->begin(); diff --git a/src/storm/storage/SparseMatrix.h b/src/storm/storage/SparseMatrix.h index 1045848f4..e5eb9ae9f 100644 --- a/src/storm/storage/SparseMatrix.h +++ b/src/storm/storage/SparseMatrix.h @@ -847,6 +847,17 @@ namespace storm { * @param b The 'right-hand side' of the problem. */ void performSuccessiveOverRelaxationStep(ValueType omega, std::vector& x, std::vector const& b) const; + + /*! + * Performs one step of the Walker-Chae technique. + * + * @param x The current solution vector. + * @param columnSums The sums of the entries of the individual columns. + * @param b The 'right-hand side' of the problem. + * @param ax A vector resulting from multiplying the current matrix with the vector x. + * @param result The vector to which to write the result. + */ + void performWalkerChaeStep(std::vector const& x, std::vector const& columnSums, std::vector const& b, std::vector const& ax, std::vector& result) const; /*! * Computes the sum of the entries in a given row. diff --git a/src/storm/utility/vector.h b/src/storm/utility/vector.h index 89fbdf1bf..b894ba4bf 100644 --- a/src/storm/utility/vector.h +++ b/src/storm/utility/vector.h @@ -822,6 +822,23 @@ namespace storm { return true; } + template + T computeSquaredNorm2Difference(std::vector const& b1, std::vector const& b2) { + STORM_LOG_ASSERT(b1.size() == b2.size(), "Vector sizes mismatch."); + + T result = storm::utility::zero(); + + auto b1It = b1.begin(); + auto b1Ite = b1.end(); + auto b2It = b2.begin(); + + for (; b1It != b1Ite; ++b1It, ++b2It) { + result += storm::utility::pow(*b1It - *b2It, 2); + } + + return result; + } + /*! * Takes the input vector and ensures that all entries conform to the bounds. */