diff --git a/.gitmodules b/.gitmodules
index a1a876d8b..1489ea541 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -20,3 +20,6 @@
 [submodule "src/3rd_party/simple-websocket-server"]
 	path = src/3rd_party/simple-websocket-server
 	url = https://github.com/marian-nmt/Simple-WebSocket-Server
+[submodule "src/3rd_party/cub"]
+	path = src/3rd_party/cub
+	url = https://github.com/NVIDIA/cub
diff --git a/CHANGELOG.md b/CHANGELOG.md
index b0d05e954..8da848a7f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 ## [Unreleased]
 
 ### Added
+- Add new operator to compute path score for lemmas for a factored vocabulary
 - Local/global sharding with MPI training via `--sharding local`
 - fp16 support for factors.
 - Correct training with fp16 via `--fp16`. 
diff --git a/src/3rd_party/cub b/src/3rd_party/cub
new file mode 160000
index 000000000..52d58a889
--- /dev/null
+++ b/src/3rd_party/cub
@@ -0,0 +1 @@
+Subproject commit 52d58a88904da39c374e44a6a8ae0e4dcca5b71a
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b47663b4e..d0c8043cf 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_definitions(-DCUB_IGNORE_DEPRECATED_CPP_DIALECT=1)
+add_definitions(-DTHRUST_IGNORE_DEPRECATED_CPP_DIALECT=1)
 add_subdirectory(3rd_party)
 
 include_directories(.)
diff --git a/src/data/factored_vocab.h b/src/data/factored_vocab.h
index 215e92f09..39b893b63 100755
--- a/src/data/factored_vocab.h
+++ b/src/data/factored_vocab.h
@@ -3,6 +3,11 @@
 // and via dynamic_cast to FactoredVocab for factored-specific things used by
 // the Embedding and Output layers.
 
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+
 #pragma once
 
 #include "common/definitions.h"
@@ -69,6 +74,7 @@ class FactoredVocab : public IVocab {
   const std::string& getFactorGroupPrefix(size_t groupIndex) const { return groupPrefixes_[groupIndex]; } // for diagnostics only
   const std::string& getFactorName(size_t groupIndex, size_t factorIndex) const { return factorVocab_[(WordIndex)(factorIndex + groupRanges_[groupIndex].first)]; }
   std::string decodeForDiagnostics(const Words& sentence) const;
+  const std::vector<std::vector<bool>>& getLemmaHasFactorGroupVector() const {return lemmaHasFactorGroup_;};
 
   static constexpr size_t FACTOR_NOT_APPLICABLE = (SIZE_MAX - 1);
   static constexpr size_t FACTOR_NOT_SPECIFIED  = (SIZE_MAX - 2);
diff --git a/src/graph/expression_graph.h b/src/graph/expression_graph.h
index adc0aeae9..3d5d7beb9 100644
--- a/src/graph/expression_graph.h
+++ b/src/graph/expression_graph.h
@@ -1,3 +1,8 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+
 #pragma once
 
 #include "common/config.h"
@@ -36,6 +41,7 @@ class Tensors {
   typedef std::unordered_map<size_t, std::vector<WExpr>> WeakMemory;
   typedef std::unordered_map<size_t, std::vector<Expr>> Memory;
 
+  std::map<std::string, Expr> memoizationMap_;
   Ptr<WeakMemory> shortterm_;  // holds all nodes for a graph
   Ptr<Memory> longterm_;  // holds memoized nodes
 
@@ -109,6 +115,19 @@ class Tensors {
     return nullptr;
   }
 
+  void rememberByName(const std::string& name, Expr e) {
+    ABORT_IF(e == nullptr, "Expression must be non-null");
+    ABORT_IF(e->type() == "param", "Not intended for graph parameters");
+    memoizationMap_[name] = e;
+    findOrRemember(e);
+  }
+
+  Expr findByName(const std::string&name) {
+    if(memoizationMap_.count(name))
+      return findOrRemember(memoizationMap_[name]);
+    return nullptr;
+  }
+
   void clear() {
     tensors_->clear();
     shortterm_->clear();
@@ -709,6 +728,23 @@ class ExpressionGraph : public std::enable_shared_from_this<ExpressionGraph> {
    */
   Ptr<TensorAllocator> getTensorAllocator() { return tensors_->getTensorAllocator(); }
 
+
+  /**
+    * @TODO Add comment
+    *
+    */
+  void rememberByName(const std::string& name, Expr e) {
+    tensors_->rememberByName(name, e);
+  }
+
+  /**
+    * @TODO Add comment
+    *
+    */
+  Expr findByName(const std::string&name) {
+    return tensors_->findByName(name);
+  }
+
   /** Clear everything apart from parameters and memoized nodes */
   void clear() {
     count_ = 0;
diff --git a/src/graph/expression_operators.cpp b/src/graph/expression_operators.cpp
index f354caabc..6abb85a5e 100644
--- a/src/graph/expression_operators.cpp
+++ b/src/graph/expression_operators.cpp
@@ -1,3 +1,7 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */ 
 #include "graph/expression_operators.h"
 #include "layers/constructors.h"
 
@@ -683,6 +687,22 @@ Expr unlikelihood(Expr logits, Expr indices) {
   return -log(gather(1.f - softmax(logits), /*axis=*/-1, indicesWithLayout));
 }
 
+Expr addFactorMaxes(Expr lemmaHasFactorGroup, std::vector<Expr> groupLosses, Expr hypIndices, size_t group0Start) {
+  if(groupLosses.size() == 1) {
+    return groupLosses[0];
+  }
+
+  int numLemmas = groupLosses[0]->shape()[-1];
+
+  std::vector<Expr> nodes({lemmaHasFactorGroup});
+  if (hypIndices) {
+    nodes.push_back(hypIndices);
+  }
+  nodes.insert(nodes.end(), groupLosses.begin(), groupLosses.end());
+  bool hasShortList = hypIndices != nullptr;
+  return Expression<AddFactorMaxesOp>(nodes, hasShortList, group0Start, numLemmas);
+}
+
 Expr plus(const std::vector<Expr>& nodes) {
   ABORT_IF(nodes.size() > 1, "Not implemented");
   return nodes[0];
diff --git a/src/graph/expression_operators.h b/src/graph/expression_operators.h
index ca0739e44..0e0034fbd 100755
--- a/src/graph/expression_operators.h
+++ b/src/graph/expression_operators.h
@@ -1,3 +1,8 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+ 
 #pragma once
 #include "graph/expression_graph.h"
 #include "graph/node_initializers.h"
@@ -610,6 +615,15 @@ Expr atleast_3d(Expr a);
  */
 Expr atleast_4d(Expr a);
 
+/**
+  * @TODO Add Comment
+  *
+  */
+Expr addFactorMaxes(Expr lemmaHasFactorGroup, std::vector<Expr> groupLosses, Expr hypIndices, size_t group0Start);
+
+// create a constant of shape a->shape() and initialize with init
+// @TODO: add a && version, to avoid a ref count. NodeInitializers are typically temps.
+// @TODO: and/or make this a template on init
 /**
  * Converts input to an expression with a least n-dimension dimensions.
  * @param a     Expression
diff --git a/src/graph/node_operators_binary.h b/src/graph/node_operators_binary.h
index 261885ec4..474ead783 100644
--- a/src/graph/node_operators_binary.h
+++ b/src/graph/node_operators_binary.h
@@ -1,3 +1,7 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
 #pragma once
 
 #include <thread>
@@ -1281,6 +1285,73 @@ struct LayerNormalizationOp : public NaryNodeOp {
   float eps_;
 };
 
+struct AddFactorMaxesOp : public NaryNodeOp {
+size_t groupStart_;
+size_t numLemmas_;
+bool hasShortlist_;
+public:
+  AddFactorMaxesOp(const std::vector<Expr>& nodes, bool hasShortlist, size_t groupStart, size_t numLemmas)
+      : NaryNodeOp(nodes, getShape(nodes, hasShortlist), commonType(std::vector<Expr>(nodes.begin() + 1 + (int)hasShortlist, nodes.end())) ) {
+    groupStart_ = groupStart;
+    numLemmas_ = numLemmas;
+    hasShortlist_ = hasShortlist;
+  }
+
+  Shape getShape(const std::vector<Expr>& nodes, bool hasShortlist) {
+    ABORT_IF(nodes.empty(), "No child nodes given");
+    int start = hasShortlist? 2 : 1;
+    return nodes[start]->shape();
+  }
+
+  NodeOps forwardOps() override {
+    int start = hasShortlist_? 2 : 1;
+    std::vector<Tensor> losses;
+    for(int i = start; i < children().size(); ++i) {
+      losses.push_back(child(i)->val());
+    }
+    return {NodeOp(
+                    AddFactorMaxes(val_,
+                                   graph()->allocator(),
+                                   child(0)->val(), // lemmaHasFactorGroupTensor
+                                   hasShortlist_? child(1)->val() : nullptr, // indices
+                                   losses,
+                                   groupStart_, numLemmas_))};
+  }
+
+  NodeOps backwardOps() override {
+    ABORT("Not Implemented for Training");
+  }
+
+  const std::string type() override { return "AddFactorMaxesOp"; }
+
+  virtual size_t hash() override {
+    size_t seed = NaryNodeOp::hash();
+    util::hash_combine(seed, hasShortlist_);
+    util::hash_combine(seed, groupStart_);
+    util::hash_combine(seed, numLemmas_);
+    return seed;
+  }
+
+  virtual bool equal(Expr node) override {
+    if(!NaryNodeOp::equal(node))
+      return false;
+    auto cnode = std::dynamic_pointer_cast<AddFactorMaxesOp>(node);
+    if(!cnode)
+      return false;
+    if(hasShortlist_ != cnode->hasShortlist_)
+      return false;
+    if(groupStart_ != cnode->groupStart_)
+      return false;
+    if(numLemmas_ != cnode->numLemmas_)
+      return false;
+    return true;
+  }
+
+private:
+  friend class SerializationHelpers; // @TODO: use the same name for this as SqrtNodeOp
+};
+
+
 struct HighwayNodeOp : public NaryNodeOp {
   HighwayNodeOp(const std::vector<Expr>& nodes) : NaryNodeOp(nodes) {}
 
diff --git a/src/layers/generic.cpp b/src/layers/generic.cpp
index d44f40206..5c9a9ccdb 100755
--- a/src/layers/generic.cpp
+++ b/src/layers/generic.cpp
@@ -1,3 +1,9 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+
+#include <algorithm>
 #include "marian.h"
 
 #include "layers/generic.h"
@@ -80,15 +86,25 @@ namespace marian {
     //  - lemma: add all maxes of applicable factors
     if (groupIndex > 0) {
       sel = sel - max(sel, -1);
-    }
-    else {
+    } else {
       auto numGroups = getNumFactorGroups();
-      for (size_t g = 1; g < numGroups; g++) {
-        auto factorMaxima = max(logits_[g]->loss(), -1); // we cast since loss is likely ce-loss which has type float32
-        auto factorMasks = constant(getFactorMasks(g, shortlist ? shortlist->indices() : std::vector<WordIndex>()));
-        sel = sel + cast(factorMaxima, sel->value_type()) * cast(factorMasks, sel->value_type()); // those lemmas that don't have a factor get multiplied with 0
+      if(numGroups > 1 && graph()->isInference() && graph()->getBackend()->getDeviceId().type == DeviceType::gpu) {
+        Expr shortlistIndices = shortlist? indices(shortlist->indices()) : nullptr;
+        Expr lemmaHasFactorGroupTensor = getLemmaHasFactorGroupTensor();
+        std::vector<Expr> groupLosses(logits_.size());
+        std::transform(logits_.begin(), logits_.end(), groupLosses.begin(), [](const Ptr<RationalLoss>& loss) -> Expr {return loss->loss();});
+
+        sel = addFactorMaxes(lemmaHasFactorGroupTensor, groupLosses, 
+                             shortlistIndices, factoredVocab_->getGroupRange(0).first);
+      } else {
+        for (size_t g = 1; g < numGroups; g++) {
+          auto factorMaxima = max(logits_[g]->loss(), -1);
+          Expr factorMasks = constant(getFactorMasks(g, shortlist ? shortlist->indices() : std::vector<WordIndex>()));
+          sel = sel + factorMaxima * factorMasks; // those lemmas that don't have a factor get multiplied with 0
+        }
       }
     }
+    
 
     // if selIdx are given, then we must reshuffle accordingly
     if (!hypIndices.empty()) // use the same function that shuffles decoder state
@@ -190,6 +206,43 @@ namespace marian {
     return res;
   }
 
+  Expr Logits::getLemmaHasFactorGroupTensor() const {
+    auto g = graph();
+    const std::string name = "lemmaHasFactorGroup";
+    auto lemmaHasFactorGroupTensor = g->findByName(name);
+
+    if(!lemmaHasFactorGroupTensor) {
+      // We want to make this a graph param so the GPU can use this to implement lemmaHasFactorGroup to avoid unneeded memcpys.
+      const auto lemmaHasFactorGroup = factoredVocab_->getLemmaHasFactorGroupVector();
+      int dimLemma = (int)lemmaHasFactorGroup.size();
+      int dimFactorGroup = (int)lemmaHasFactorGroup[0].size();
+
+      for(const auto&  lemma : lemmaHasFactorGroup) {
+        // Paranoid check - Instead of aborting I think we can pad here and copy the array
+        ABORT_IF(lemma.size() != dimFactorGroup, "All groups must be the same size");
+      }
+
+      auto initFunc = [lemmaHasFactorGroup, dimLemma, dimFactorGroup] (Tensor t) {
+        std::vector<int8_t> flattened(dimLemma * dimFactorGroup);
+        int row = 0;
+        for(const auto& v : lemmaHasFactorGroup) {
+          int offset = row * dimFactorGroup;
+          for(int i = 0; i < v.size(); ++i) {
+            flattened[offset + i] = static_cast<int8_t>(v[i]);
+          }
+          ++row;
+        }
+        t->set(flattened);
+      };
+
+      lemmaHasFactorGroupTensor = graph()->constant({dimLemma, dimFactorGroup}, inits::fromLambda(initFunc), Type::int8);
+      lemmaHasFactorGroupTensor->setMemoize(true);
+      g->rememberByName(name, lemmaHasFactorGroupTensor);
+    }
+    ABORT_IF(!lemmaHasFactorGroupTensor, "Lemma has factor group tensor undefined?");
+    return lemmaHasFactorGroupTensor;
+  }
+
   Logits Logits::applyUnaryFunction(const std::function<Expr(Expr)>& f) const { // clone this but apply f to all loss values
     std::vector<Ptr<RationalLoss>> newLogits;
     for (const auto& l : logits_)
diff --git a/src/layers/generic.h b/src/layers/generic.h
index f47bb45e2..71a1fa54a 100755
--- a/src/layers/generic.h
+++ b/src/layers/generic.h
@@ -1,3 +1,8 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+ 
 #pragma once
 
 #include "marian.h"
@@ -142,6 +147,8 @@ class Logits {
     template<typename T> Expr constant(const std::vector<T>& data) const { return constant(Shape{(int)data.size()}, data); } // same as constant() but assuming vector
     Expr indices(const std::vector<uint32_t>& data) const { return graph()->indices(data); } // actually the same as constant(data) for this data type
     std::vector<float> getFactorMasks(size_t factorGroup, const std::vector<WordIndex>& indices) const;
+    Expr getLemmaHasFactorGroupTensor() const;
+    Expr addFactorMaxesHelper(Expr indices=nullptr) const;
 private:
     // members
     // @TODO: we don't use the RationalLoss component anymore, can be removed again, and replaced just by the Expr
diff --git a/src/tensors/cpu/tensor_operators.cpp b/src/tensors/cpu/tensor_operators.cpp
index 1191a2bec..705c1798f 100755
--- a/src/tensors/cpu/tensor_operators.cpp
+++ b/src/tensors/cpu/tensor_operators.cpp
@@ -3,6 +3,11 @@
  *   SPDX-License-Identifier: MIT
  */
 
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+
 #include "tensors/tensor_operators.h"
 #include "tensors/cpu/backend.h"
 #include "tensors/allocator.h"
@@ -24,6 +29,16 @@ void IsNaN(const Tensor /*in*/, Ptr<Allocator> /*allocator*/, bool& /*isNaN*/, b
   ABORT("Not implemented");
 }
 
+void AddFactorMaxes(Tensor /*out*/,
+                    Ptr<Allocator> /*allocator*/,
+                    const Tensor /*lemmaHasFactorGroupTensor*/,
+                    const Tensor /*indices*/,
+                    const std::vector<marian::Tensor>& /*groupLosses*/,
+                    size_t /*groupStart*/,
+                    size_t /*numLemmas*/) {
+  ABORT("AddFactorMaxes not implemented on CPU");
+}
+
 template <bool add, typename To, typename From>
 void CopyCastTo(To* out, const From* in, int length) {
   for(int i = 0; i < length; ++i)
diff --git a/src/tensors/gpu/tensor_operators.cu b/src/tensors/gpu/tensor_operators.cu
index 97f0cdfe0..286e97c59 100644
--- a/src/tensors/gpu/tensor_operators.cu
+++ b/src/tensors/gpu/tensor_operators.cu
@@ -1,3 +1,10 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+
+#include <limits>
+#include "common/logging.h"
 #include "common/types.h"
 #include "tensors/tensor_operators.h"
 
@@ -9,6 +16,20 @@
 
 #include "tensors/gpu/add_all.h"
 
+#if COMPILE_FP16
+#include <cuda_fp16.h>
+__device__ __forceinline__ half max(const half a, const half b) {
+  return a > b ? a : b;
+}
+#endif
+
+
+#if CUDA_VERSION >= 11000
+#include <cub/cub.cuh>
+#else
+#include "cub/cub/cub.cuh"
+#endif
+
 namespace marian {
 
 namespace gpu {
@@ -834,7 +855,7 @@ void LogSoftmax(Tensor out, Tensor in) {
   size_t m = out->shape().elements() / out->shape().back();
   size_t k = out->shape().back();
 
-  int blocks = std::min(MAX_BLOCKS, (int)m);
+  int blocks = (int)m;
   int threads = std::min(MAX_THREADS, (int)k);
   int shared = sizeof(float) * threads; // use float32 as accumulation type
 
@@ -3032,5 +3053,190 @@ void PoolingWithMaskingBackward(Tensor adj,
                                            width,
                                            lastWidth);
 }
+
+template <typename T>
+struct ptrInnerDimPair {
+  T* ptr;
+  int innerDim;
+};
+
+template <typename T>
+__global__ void gAddFactorMaxesPhase1(T* maxes, // [#blocks, numGroups - 1] 
+                                      ptrInnerDimPair<T>* lossAndLastDimSize, 
+                                      size_t numGroups, size_t sizeWithoutInnerDim, 
+                                      T minimal) {
+
+  typedef cub::WarpReduce<T> WarpReduce;
+  __shared__ typename WarpReduce::TempStorage temp_storage[32]; // Max thread size (1024) divided by warp size (32)
+
+
+  // We exploit the fact that the factor groups (except for the lemmas) tends to be small. Therefore, we perform
+  // all of the reductions across factor groups in parallel by splitting the work across warps.
+  const int lane = threadIdx.x % warpSize;
+  const int wid = threadIdx.x / warpSize;
+  const int numWarps = blockDim.x / warpSize; 
+
+  // Each block computes the maxes for each secondary factor within a given row. For example, if we have 6 factor groups,
+  // (so lemma[fg0] + 5 secondary groups) block 0 will compute the max in row 0 for fg1 to fg6. These is currently performed
+  // in parallel by assigning fg1 to warp 0. fg2 to warp 1 etc. If there are not enough warps to do all the factor groups then 
+  // multiple batched iterations are needed to calculate the factor maxes.
+  for(int warpFactorGroup = wid + 1; warpFactorGroup < numGroups; warpFactorGroup+=numWarps) {
+    T* groupLosses = lossAndLastDimSize[warpFactorGroup].ptr;
+    int groupLossesInnerDim = lossAndLastDimSize[warpFactorGroup].innerDim;
+    int groupLossSize = groupLossesInnerDim * sizeWithoutInnerDim;
+
+    T factorMaxima = minimal;
+    // We start by each warp computing offset to the row its block is responsible for in the loss
+
+    // This is per warp since groupLossesInnerDim is warp specific and the groupLossesPtr is also warp specific.
+    const int blockRowStart = blockIdx.x * groupLossesInnerDim;
+    for(int lossCol = lane; lossCol < groupLossesInnerDim; lossCol += warpSize) {
+      int warpOffset = blockRowStart + lossCol;
+      if (warpOffset < groupLossSize) {
+        factorMaxima = max(factorMaxima, groupLosses[warpOffset]);
+      }
+    }
+
+    // Each warp reduces the accumulated values. Then each each writes its final value to global for phase 2 of this op
+    factorMaxima = WarpReduce(temp_storage[wid]).Reduce(factorMaxima, cub::Max());
+    if (lane == 0) {
+      T* factorMaximasInBlock = maxes + blockIdx.x * (numGroups - 1);
+      factorMaximasInBlock[warpFactorGroup - 1] = factorMaxima;
+    }
+  }
+}
+
+
+  template <typename T>
+__global__ void gAddFactorMaxesPhase2(T* out, const int8_t* const lemmaHasFactorGroup, 
+                                      const IndexType* const indices, 
+                                      const T* lemmaLosses, 
+                                      size_t numGroups, 
+                                      const int sizeWithoutInnerDim, 
+                                      const T* maxes,
+                                      size_t groupStart, 
+                                      int lemmaHasFactorGroupWidth, 
+                                      size_t numLemmas) {
+
+  // Each thread is given an element in the output to construct. First we compute every thread's offset from the 
+  // base pointer along with the row and column indicies for the output cell the thread is responsible for constructing.
+  const int offsetFromBasePtr = blockIdx.x * blockDim.x + threadIdx.x;
+  const int rowForThread = offsetFromBasePtr / numLemmas;
+  const int colForThread = offsetFromBasePtr % numLemmas;
+
+  // We keep the maxes in shared memory since they are read repeatedly by each thread in a block
+  const T* globalFactorMaxesRow = maxes + rowForThread * (numGroups - 1);
+
+  // Compute a write guard for the output array
+  const int outputSize = numLemmas * sizeWithoutInnerDim;
+
+  // Now we start accumulating the thread's sum by initializing to the lemma value. The out tensor and lemmaLosses
+  // tensor are the same size so the output guard is used to prevent invalid reads.
+  if (offsetFromBasePtr < outputSize) {
+    T threadSum = lemmaLosses[offsetFromBasePtr];
+
+    // Each block accumulates its output locally then writes the final result to global. We use __ldg since the 
+    // data being read does not change throughout the life of the kernel so the constant cache can be used.
+    // This is more of a compiler hint than a necessity.
+    for(int g = 1; g < (int)numGroups; ++g) {
+      int lemma = indices? indices[colForThread] - groupStart : colForThread;
+      T factorMask = static_cast<T>(__ldg(&lemmaHasFactorGroup[lemma * lemmaHasFactorGroupWidth + g]));
+      T factorMaxima = globalFactorMaxesRow[g-1];
+      threadSum += (factorMask * factorMaxima);
+    }
+
+    // Finally, the accumulated sum is written back out to global memory.
+    out[offsetFromBasePtr] = threadSum;
+  }
+  
+}
+
+void AddFactorMaxes(Tensor out,
+                    Ptr<Allocator> allocator,
+                    const Tensor lemmaHasFactorGroupTensor,
+                    const Tensor indices,
+                    const std::vector<marian::Tensor>& groupLosses,
+                    size_t groupStart,
+                    size_t numLemmas) {
+
+  cudaSetDevice(out->getDeviceId().no);
+
+  ABORT_IF(lemmaHasFactorGroupTensor->type() != Type::int8, "lemmaHasFactor group tensor wrong type");
+  ABORT_IF(out->shape()[-1] != (int) numLemmas, "Output shape {} or numLemmas {} incorrect", out->shape(), numLemmas);
+
+  const int sizeWithoutInnerDim = out->shape().elements() / out->shape()[-1];
+  constexpr int warpSize = 32;
+  const int blocksPhase1 = sizeWithoutInnerDim;
+  const int threadsPhase1 = std::min(MAX_THREADS, (int)(warpSize * (groupLosses.size() - 1)) );
+
+  const int totalPhase2Threads = out->shape().elements();
+  const int threadsPhase2 = std::min(MAX_THREADS, out->shape()[-1]);
+  const int blocksPhase2 = (totalPhase2Threads + threadsPhase2 - 1) / threadsPhase2;
+
+  IPtr<MemoryPiece> mp_ptrs;
+  IPtr<MemoryPiece> factorMaxes;
+  if (out->type() == Type::float32) {
+    mp_ptrs = allocator->alloc<ptrInnerDimPair<float>>(groupLosses.size()); 
+    factorMaxes = allocator->alloc<float>(sizeWithoutInnerDim * (groupLosses.size() - 1));
+    ptrInnerDimPair<float>* dest = mp_ptrs->data<ptrInnerDimPair<float>>();
+
+    std::vector<ptrInnerDimPair<float>> lossPtrs;
+    for(const auto& t : groupLosses) {
+      ptrInnerDimPair<float> pair = {t->data<float>(), t->shape()[-1]};
+      lossPtrs.push_back(pair);
+    }
+
+    // Async just so that call is issued in per thread default stream
+    CUDA_CHECK(cudaMemcpyAsync(dest, lossPtrs.data(), lossPtrs.size() * sizeof(ptrInnerDimPair<float>), cudaMemcpyHostToDevice, 0));
+    gAddFactorMaxesPhase1<<<blocksPhase1, threadsPhase1>>>(factorMaxes->data<float>(), // [#rows, numGroups - 1] 
+                                                           dest, 
+                                                           groupLosses.size(), 
+                                                           sizeWithoutInnerDim, 
+                                                           -std::numeric_limits<float>::infinity());
+
+    gAddFactorMaxesPhase2<float><<<blocksPhase2, threadsPhase2>>>(out->data<float>(), 
+                                                                  lemmaHasFactorGroupTensor->data<int8_t>(), 
+                                                                  indices? indices->data<IndexType>(): nullptr,
+                                                                  groupLosses[0]->data<float>(),
+                                                                  groupLosses.size(), sizeWithoutInnerDim, 
+                                                                  factorMaxes->data<float>(), groupStart, 
+                                                                  lemmaHasFactorGroupTensor->shape()[1], 
+                                                                  numLemmas);
+  #if COMPILE_FP16
+  } else if (out->type() == Type::float16) {
+    mp_ptrs = allocator->alloc<ptrInnerDimPair<half>>(groupLosses.size()); 
+    factorMaxes = allocator->alloc<half>(sizeWithoutInnerDim * (groupLosses.size() - 1));
+    ptrInnerDimPair<half>* dest = mp_ptrs->data<ptrInnerDimPair<half>>();
+
+    std::vector<ptrInnerDimPair<half>> lossPtrs;
+    for(const auto& t : groupLosses) {
+      ptrInnerDimPair<half> pair = {t->data<half>(), t->shape()[-1]};
+      lossPtrs.push_back(pair);
+    }
+
+    // Async just so that call is issued in per thread default stream
+    CUDA_CHECK(cudaMemcpyAsync(dest, lossPtrs.data(), lossPtrs.size() * sizeof(ptrInnerDimPair<half>), cudaMemcpyHostToDevice, 0));
+    gAddFactorMaxesPhase1<<<blocksPhase1, threadsPhase1>>>(factorMaxes->data<half>(), // [#rows, numGroups - 1] 
+                                                           dest, 
+                                                           groupLosses.size(), 
+                                                           sizeWithoutInnerDim, 
+                                                           __float2half(-std::numeric_limits<float>::infinity()) );
+
+    gAddFactorMaxesPhase2<half><<<blocksPhase2, threadsPhase2>>>(out->data<half>(), 
+                                                                 lemmaHasFactorGroupTensor->data<int8_t>(), 
+                                                                 indices? indices->data<IndexType>(): nullptr,
+                                                                 groupLosses[0]->data<half>(),
+                                                                 groupLosses.size(), sizeWithoutInnerDim, 
+                                                                 factorMaxes->data<half>(), groupStart, 
+                                                                 lemmaHasFactorGroupTensor->shape()[1], 
+                                                                 numLemmas);
+  #endif
+  } else {
+    ABORT("AddFactorMaxes not implemented for type {}", out->type());
+  }
+  allocator->free(mp_ptrs);
+  allocator->free(factorMaxes);
+}
+
 }  // namespace gpu
 }  // namespace marian
diff --git a/src/tensors/tensor_operators.h b/src/tensors/tensor_operators.h
index 83bce8194..3299f55c1 100644
--- a/src/tensors/tensor_operators.h
+++ b/src/tensors/tensor_operators.h
@@ -1,3 +1,8 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+ 
 #pragma once
 
 #include "common/definitions.h"
@@ -123,6 +128,8 @@ DISPATCH4(ShiftGrad, marian::Tensor, marian::Tensor, marian::Shape, bool)
 
 DISPATCH3(Concatenate, marian::Tensor, const std::vector<marian::Tensor>&, int)
 
+DISPATCH7(AddFactorMaxes, marian::Tensor, Ptr<Allocator>, const marian::Tensor, const marian::Tensor, const std::vector<marian::Tensor>&, size_t, size_t)
+
 // clang-format on
 
 // Bernoulli(tensor, 0.5f, 2.f, -1.f) generates a tensor composed of 50% of 1 and 50% of -1.
diff --git a/src/translator/beam_search.cpp b/src/translator/beam_search.cpp
index 5c1989a68..7ffa03b03 100755
--- a/src/translator/beam_search.cpp
+++ b/src/translator/beam_search.cpp
@@ -1,3 +1,8 @@
+/* Part of this file was contributed by NVIDIA under license:
+ *   Copyright (C) 2020 NVIDIA Corporation
+ *   SPDX-License-Identifier: MIT
+ */
+
 #include "translator/beam_search.h"
 
 #include "data/factored_vocab.h"
@@ -396,7 +401,8 @@ Histories BeamSearch::search(Ptr<ExpressionGraph> graph, Ptr<data::CorpusBatch>
         }
         if(factorGroup == 0)
           currentDimBatch = (IndexType) batchIndices.size(); // keep batch size constant for all factor groups in a time step
-        prevPathScores = graph->constant({(int)maxBeamSize, 1, (int)currentDimBatch, 1}, inits::fromVector(prevScores));
+        // Avoid unnecessary memcpy on GPU  
+        if(anyCanExpand) prevPathScores = graph->constant({(int)maxBeamSize, 1, (int)currentDimBatch, 1}, inits::fromVector(prevScores));
       }
       if (!anyCanExpand) // all words cannot expand this factor: skip
         continue;