diff --git a/src/storm/solver/Multiplier.h b/src/storm/solver/Multiplier.h
index a94934ae0..5f9036022 100644
--- a/src/storm/solver/Multiplier.h
+++ b/src/storm/solver/Multiplier.h
@@ -144,8 +144,6 @@ namespace storm {
             ~MultiplierFactory() = default;
 
             std::unique_ptr<Multiplier<ValueType>> create(Environment const& env, storm::storage::SparseMatrix<ValueType> const& matrix);
-            
-            
         };
         
     }
diff --git a/src/storm/storage/SparseMatrix.cpp b/src/storm/storage/SparseMatrix.cpp
index d909b76b2..3f0f25828 100644
--- a/src/storm/storage/SparseMatrix.cpp
+++ b/src/storm/storage/SparseMatrix.cpp
@@ -1501,7 +1501,7 @@ namespace storm {
                 multiplyWithVectorParallel(vector, tmpVector);
                 result = std::move(tmpVector);
             } else {
-                tbb::parallel_for(tbb::blocked_range<index_type>(0, result.size(), 10), TbbMultAddFunctor<ValueType>(columnsAndValues, rowIndications, vector, result, summand));
+                tbb::parallel_for(tbb::blocked_range<index_type>(0, result.size(), 100), TbbMultAddFunctor<ValueType>(columnsAndValues, rowIndications, vector, result, summand));
             }
         }
 #endif
@@ -1586,6 +1586,7 @@ namespace storm {
         
         template<typename ValueType>
         void SparseMatrix<ValueType>::multiplyAndReduceForward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+            bool min = dir == OptimizationDirection::Minimize;
             auto elementIt = this->begin();
             auto rowGroupIt = rowGroupIndices.begin();
             auto rowIt = rowIndications.begin();
@@ -1597,13 +1598,15 @@ namespace storm {
             if (choices) {
                 choiceIt = choices->begin();
             }
+
+            // Variables for correctly tracking choices (only update if new choice is strictly better).
+            ValueType oldSelectedChoiceValue;
+            uint64_t selectedChoice;
             
+            uint64_t currentRow = 0;
             for (auto resultIt = result.begin(), resultIte = result.end(); resultIt != resultIte; ++resultIt, ++choiceIt, ++rowGroupIt) {
                 ValueType currentValue = storm::utility::zero<ValueType>();
-                if (choices) {
-                    *choiceIt = 0;
-                }
-                
+
                 // Only multiply and reduce if there is at least one row in the group.
                 if (*rowGroupIt < *(rowGroupIt + 1)) {
                     if (summand) {
@@ -1615,30 +1618,43 @@ namespace storm {
                         currentValue += elementIt->getValue() * vector[elementIt->getColumn()];
                     }
                     
+                    if (choices) {
+                        selectedChoice = 0;
+                        if (*choiceIt == 0) {
+                            oldSelectedChoiceValue = currentValue;
+                        }
+                    }
+                    
                     ++rowIt;
+                    ++currentRow;
                     
-                    for (; static_cast<uint_fast64_t>(std::distance(rowIndications.begin(), rowIt)) < *(rowGroupIt + 1); ++rowIt) {
+                    for (; currentRow < *(rowGroupIt + 1); ++rowIt, ++currentRow) {
                         ValueType newValue = summand ? *summandIt : storm::utility::zero<ValueType>();
                         for (auto elementIte = this->begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) {
                             newValue += elementIt->getValue() * vector[elementIt->getColumn()];
                         }
                         
-                        if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) {
+                        if (choices && currentRow == *choiceIt + *rowGroupIt) {
+                            oldSelectedChoiceValue = newValue;
+                        }
+                        
+                        if (min ? newValue < currentValue : newValue > currentValue) {
                             currentValue = newValue;
                             if (choices) {
-                                *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt;
+                                selectedChoice = currentRow - *rowGroupIt;
                             }
                         }
                         if (summand) {
                             ++summandIt;
                         }
                     }
-                } else if (choices) {
-                    *choiceIt = 0;
+
+                    // Finally write value to target vector.
+                    *resultIt = currentValue;
+                    if (choices && (min ? currentValue < oldSelectedChoiceValue : currentValue > oldSelectedChoiceValue)) {
+                        *choiceIt = selectedChoice;
+                    }
                 }
-                
-                // Finally write value to target vector.
-                *resultIt = currentValue;
             }
         }
 
@@ -1651,6 +1667,7 @@ namespace storm {
         
         template<typename ValueType>
         void SparseMatrix<ValueType>::multiplyAndReduceBackward(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
+            bool min = dir == OptimizationDirection::Minimize;
             auto elementIt = this->end() - 1;
             auto rowGroupIt = rowGroupIndices.end() - 2;
             auto rowIt = rowIndications.end() - 2;
@@ -1662,13 +1679,15 @@ namespace storm {
             if (choices) {
                 choiceIt = choices->end() - 1;
             }
-            
+
+            // Variables for correctly tracking choices (only update if new choice is strictly better).
+            ValueType oldSelectedChoiceValue;
+            uint64_t selectedChoice;
+
+            uint64_t currentRow = this->getRowCount() - 1;
             for (auto resultIt = result.end() - 1, resultIte = result.begin() - 1; resultIt != resultIte; --resultIt, --choiceIt, --rowGroupIt) {
-                if (choices) {
-                    *choiceIt = 0;
-                }
                 ValueType currentValue = storm::utility::zero<ValueType>();
-                
+
                 // Only multiply and reduce if there is at least one row in the group.
                 if (*rowGroupIt < *(rowGroupIt + 1)) {
                     if (summand) {
@@ -1680,27 +1699,38 @@ namespace storm {
                         currentValue += elementIt->getValue() * vector[elementIt->getColumn()];
                     }
                     if (choices) {
-                        *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt;
+                        selectedChoice = currentRow - *rowGroupIt;
+                        if (*choiceIt == selectedChoice) {
+                            oldSelectedChoiceValue = currentValue;
+                        }
                     }
                     --rowIt;
+                    --currentRow;
                     
-                    for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, ++i, --summandIt) {
+                    for (uint64_t i = *rowGroupIt + 1, end = *(rowGroupIt + 1); i < end; --rowIt, --currentRow, ++i, --summandIt) {
                         ValueType newValue = summand ? *summandIt : storm::utility::zero<ValueType>();
                         for (auto elementIte = this->begin() + *rowIt - 1; elementIt != elementIte; --elementIt) {
                             newValue += elementIt->getValue() * vector[elementIt->getColumn()];
                         }
                         
-                        if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) {
+                        if (choices && currentRow == *choiceIt + *rowGroupIt) {
+                            oldSelectedChoiceValue = newValue;
+                        }
+                        
+                        if (min ? newValue < currentValue : newValue > currentValue) {
                             currentValue = newValue;
                             if (choices) {
-                                *choiceIt = std::distance(rowIndications.begin(), rowIt) - *rowGroupIt;
+                                selectedChoice = currentRow - *rowGroupIt;
                             }
                         }
                     }
+
+                    // Finally write value to target vector.
+                    *resultIt = currentValue;
+                    if (choices && (min ? currentValue < oldSelectedChoiceValue : currentValue > oldSelectedChoiceValue)) {
+                        *choiceIt = selectedChoice;
+                    }
                 }
-                
-                // Finally write value to target vector.
-                *resultIt = currentValue;
             }
         }
         
@@ -1724,6 +1754,8 @@ namespace storm {
             }
             
             void operator()(tbb::blocked_range<index_type> const& range) const {
+                bool min = dir == storm::OptimizationDirection::Minimize;
+
                 auto groupIt = rowGroupIndices.begin() + range.begin();
                 auto groupIte = rowGroupIndices.begin() + range.end();
                 
@@ -1740,11 +1772,12 @@ namespace storm {
                 
                 auto resultIt = result.begin() + range.begin();
                 
+                // Variables for correctly tracking choices (only update if new choice is strictly better).
+                ValueType oldSelectedChoiceValue;
+                uint64_t selectedChoice;
+
+                uint64_t currentRow = *groupIt;
                 for (; groupIt != groupIte; ++groupIt, ++resultIt, ++choiceIt) {
-                    if (choices) {
-                        *choiceIt = 0;
-                    }
-                    
                     ValueType currentValue = storm::utility::zero<ValueType>();
                     
                     // Only multiply and reduce if there is at least one row in the group.
@@ -1758,25 +1791,40 @@ namespace storm {
                             currentValue += elementIt->getValue() * x[elementIt->getColumn()];
                         }
                         
+                        if (choices) {
+                            selectedChoice = 0;
+                            if (*choiceIt == 0) {
+                                oldSelectedChoiceValue = currentValue;
+                            }
+                        }
+                        
                         ++rowIt;
+                        ++currentRow;
                         
-                        for (; static_cast<uint_fast64_t>(std::distance(rowIndications.begin(), rowIt)) < *(groupIt + 1); ++rowIt, ++summandIt) {
+                        for (; currentRow < *(groupIt + 1); ++rowIt, ++currentRow, ++summandIt) {
                             ValueType newValue = summand ? *summandIt : storm::utility::zero<ValueType>();
                             for (auto elementIte = columnsAndEntries.begin() + *(rowIt + 1); elementIt != elementIte; ++elementIt) {
                                 newValue += elementIt->getValue() * x[elementIt->getColumn()];
                             }
                             
-                            if ((dir == OptimizationDirection::Minimize && newValue < currentValue) || (dir == OptimizationDirection::Maximize && newValue > currentValue)) {
+                            if (choices && currentRow == *choiceIt + *groupIt) {
+                                oldSelectedChoiceValue = newValue;
+                            }
+
+                            if (min ? newValue < currentValue : newValue > currentValue) {
                                 currentValue = newValue;
                                 if (choices) {
-                                    *choiceIt = std::distance(rowIndications.begin(), rowIt) - *groupIt;
+                                    selectedChoice = currentRow - *groupIt;
                                 }
                             }
                         }
+                        
+                        // Finally write value to target vector.
+                        *resultIt = currentValue;
+                        if (choices && (min ? currentValue < oldSelectedChoiceValue : currentValue > oldSelectedChoiceValue)) {
+                            *choiceIt = selectedChoice;
+                        }
                     }
-                    
-                    // Finally write value to target vector.
-                    *resultIt = currentValue;
                 }
             }
             
@@ -1793,7 +1841,7 @@ namespace storm {
         
         template<typename ValueType>
         void SparseMatrix<ValueType>::multiplyAndReduceParallel(OptimizationDirection const& dir, std::vector<uint64_t> const& rowGroupIndices, std::vector<ValueType> const& vector, std::vector<ValueType> const* summand, std::vector<ValueType>& result, std::vector<uint_fast64_t>* choices) const {
-            tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 10), TbbMultAddReduceFunctor<ValueType>(dir, rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices));
+            tbb::parallel_for(tbb::blocked_range<index_type>(0, rowGroupIndices.size() - 1, 100), TbbMultAddReduceFunctor<ValueType>(dir, rowGroupIndices, columnsAndValues, rowIndications, vector, result, summand, choices));
         }
         
 #ifdef STORM_HAVE_CARL