diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 0b7d03b95..427eec25f 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -47,11 +47,8 @@ namespace storm {
                 /**
                  * Initialize results vectors. vLowerOld, vUpperOld and wUpper[k=N] are initialized with zeros.
                  */
-                UnifPlusVectors(uint64_t steps, uint64_t noStates) : numberOfStates(noStates), resLowerOld(numberOfStates, storm::utility::zero<ValueType>()), resLowerNew(numberOfStates, -1), resUpperOld(numberOfStates, storm::utility::zero<ValueType>()), resUpperNew(numberOfStates, storm::utility::zero<ValueType>()) {
-                    // For wUpper we have to keep track of all previous results
-                    wUpper = std::vector<std::vector<ValueType>>(steps+1, std::vector<ValueType>(numberOfStates, -1));
-                    // Initialize entries for step N with zeros
-                    std::fill(wUpper[steps].begin(), wUpper[steps].end(), storm::utility::zero<ValueType>());
+                UnifPlusVectors(uint64_t steps, uint64_t noStates) : numberOfStates(noStates), steps(steps), resLowerOld(numberOfStates, storm::utility::zero<ValueType>()), resLowerNew(numberOfStates, -1), resUpper(numberOfStates, storm::utility::zero<ValueType>()), wUpperOld(numberOfStates, storm::utility::zero<ValueType>()), wUpperNew(numberOfStates, -1) {
+                    // Intentionally left empty
                 }
 
                 /**
@@ -60,23 +57,24 @@ namespace storm {
                 void prepareNewIteration() {
                     resLowerOld.swap(resLowerNew);
                     std::fill(resLowerNew.begin(), resLowerNew.end(), -1);
-                    resUpperOld.swap(resUpperNew);
-                    std::fill(resUpperNew.begin(), resUpperNew.end(), storm::utility::zero<ValueType>());
+                    wUpperOld.swap(wUpperNew);
+                    std::fill(wUpperNew.begin(), wUpperNew.end(), -1);
                 }
 
                 uint64_t numberOfStates;
+                uint64_t steps;
                 std::vector<ValueType> resLowerOld;
                 std::vector<ValueType> resLowerNew;
-                std::vector<ValueType> resUpperOld;
-                std::vector<ValueType> resUpperNew;
-                std::vector<std::vector<ValueType>> wUpper;
+                std::vector<ValueType> resUpper;
+                std::vector<ValueType> wUpperOld;
+                std::vector<ValueType> wUpperNew;
             };
             
             template<typename ValueType>
             void calculateUnifPlusVector(Environment const& env, uint64_t k, uint64_t state, bool calcLower, ValueType lambda, uint64_t numberOfProbabilisticChoices, std::vector<std::vector<ValueType>> const & relativeReachability, OptimizationDirection dir, UnifPlusVectors<ValueType>& unifVectors, storm::storage::SparseMatrix<ValueType> const& fullTransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, std::unique_ptr<storm::solver::MinMaxLinearEquationSolver<ValueType>> const& solver, storm::utility::numerical::FoxGlynnResult<ValueType> const& poisson, bool cycleFree) {
                 // Set reference to acutal vector
-                std::vector<ValueType>& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpper[k+1];
-                std::vector<ValueType>& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpper[k];
+                std::vector<ValueType>& resVectorOld = calcLower ? unifVectors.resLowerOld : unifVectors.wUpperOld;
+                std::vector<ValueType>& resVectorNew = calcLower ? unifVectors.resLowerNew : unifVectors.wUpperNew;
 
                 if (resVectorNew[state] != -1) {
                     // Result already calculated.
@@ -84,7 +82,7 @@ namespace storm {
                 }
                 
                 auto numberOfStates = fullTransitionMatrix.getRowGroupCount();
-                uint64_t N = unifVectors.wUpper.size() - 1;
+                uint64_t N = unifVectors.steps;
                 auto const& rowGroupIndices = fullTransitionMatrix.getRowGroupIndices();
                 ValueType res;
                 
@@ -389,26 +387,20 @@ namespace storm {
                             // Calculate results for lower bound and wUpper
                             calculateUnifPlusVector(env, k, state, true, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree);
                             calculateUnifPlusVector(env, k, state, false, lambda, numberOfProbabilisticChoices, relativeReachabilities, dir, unifVectors, fullTransitionMatrix, markovianAndGoalStates, psiStates, solver, foxGlynnResult, cycleFree);
-                        }
-
-                        if (k == 0) {
                             // Calculate result for upper bound
-                            // resUpperNew was already initialized with zeros
-                            uint64_t left = std::max(foxGlynnResult.left, (uint64_t)(k));
-                            uint64_t right = std::min(foxGlynnResult.right, N-1);
-                            for (uint64_t state = 0; state < numberOfStates; ++state) {
-                                for (uint64_t i = left; i <= right; ++i) {
-                                    STORM_LOG_ASSERT(unifVectors.wUpper[N-1-(i-k)][state] != -1, "wUpper was not computed before.");
-                                    unifVectors.resUpperNew[state] += foxGlynnResult.weights[i - foxGlynnResult.left] * unifVectors.wUpper[N-1-(i-k)][state];
-                                }
+                            uint64_t index = N-1-k;
+                            if (index >= foxGlynnResult.left && index <= foxGlynnResult.right) {
+                                STORM_LOG_ASSERT(unifVectors.wUpperNew[state] != -1, "wUpper was not computed before.");
+                                unifVectors.resUpper[state] += foxGlynnResult.weights[index - foxGlynnResult.left] * unifVectors.wUpperNew[state];
                             }
                         }
+
                     }
 
                     // Only iterate over result vector, as the results can only get more precise.
                     maxNorm = storm::utility::zero<ValueType>();
                     for (uint64_t i = 0; i < numberOfStates; i++){
-                        ValueType diff = storm::utility::abs(unifVectors.resUpperNew[i] - unifVectors.resLowerNew[i]);
+                        ValueType diff = storm::utility::abs(unifVectors.resUpper[i] - unifVectors.resLowerNew[i]);
                         maxNorm = std::max(maxNorm, diff);
                     }