diff --git a/urbansim/urbanchoice/mnl.py b/urbansim/urbanchoice/mnl.py
index f9ed73c8..463eb252 100644
--- a/urbansim/urbanchoice/mnl.py
+++ b/urbansim/urbanchoice/mnl.py
@@ -35,6 +35,9 @@ def mnl_probs(data, beta, numalts):
         raise Exception("Number of alternatives is zero")
     utilities.reshape(numalts, utilities.size() // numalts)
 
+    # https://stats.stackexchange.com/questions/304758/softmax-overflow
+    utilities = utilities.subtract(utilities.max(0))
+
     exponentiated_utility = utilities.exp(inplace=True)
     if clamp:
         exponentiated_utility.inftoval(1e20)
@@ -245,6 +248,10 @@ def mnl_estimate(data, chosen, numalts, GPU=False, coeffrange=(-3, 3),
                                                    approx_grad=False,
                                                    bounds=bounds
                                                    )
+
+    if bfgs_result[2]['warnflag'] > 0:
+        logger.warn("mnl did not converge correctly: %s",  bfgs_result)
+
     beta = bfgs_result[0]
     stderr = mnl_loglik(
         beta, data, chosen, numalts, weights, stderr=1, lcgrad=lcgrad)
diff --git a/urbansim/urbanchoice/pmat.py b/urbansim/urbanchoice/pmat.py
index d25148b2..96e1e889 100644
--- a/urbansim/urbanchoice/pmat.py
+++ b/urbansim/urbanchoice/pmat.py
@@ -79,6 +79,12 @@ def cumsum(self, axis):
         # elif self.typ == 'cuda':
         #  return PMAT(misc.cumsum(self.mat,axis=axis))
 
+    def max(self, axis):
+        if self.typ == 'numpy':
+            return PMAT(np.max(self.mat, axis=axis))
+        elif self.typ == 'cuda':
+            return PMAT(self.mat.max(axis=axis))
+
     def argmax(self, axis):
         if self.typ == 'numpy':
             return PMAT(np.argmax(self.mat, axis=axis))