diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py
index fa7b1556e4e..61b0575cca2 100644
--- a/sklearn/linear_model/logistic.py
+++ b/sklearn/linear_model/logistic.py
@@ -22,7 +22,7 @@ from ..svm.base import _fit_liblinear
 from ..utils import check_array, check_consistent_length, compute_class_weight
 from ..utils import check_random_state
 from ..utils.extmath import (logsumexp, log_logistic, safe_sparse_dot,
-                             squared_norm)
+                             softmax, squared_norm)
 from ..utils.optimize import newton_cg
 from ..utils.validation import (as_float_array, DataConversionWarning,
                                 check_X_y, NotFittedError)
@@ -1111,11 +1111,7 @@ class LogisticRegression(BaseEstimator, LinearClassifierMixin,
         if calculate_ovr:
             return super(LogisticRegression, self)._predict_proba_lr(X)
         else:
-            prob = self.decision_function(X)
-            np.exp(prob, prob)
-            sum_prob = np.sum(prob, axis=1).reshape((-1, 1))
-            prob /= sum_prob
-            return prob
+            return softmax(self.decision_function(X), copy=False)
 
     def predict_log_proba(self, X):
         """Log of probability estimates.
diff --git a/sklearn/utils/extmath.py b/sklearn/utils/extmath.py
index ff0ea1a0658..717b30926ee 100644
--- a/sklearn/utils/extmath.py
+++ b/sklearn/utils/extmath.py
@@ -615,6 +615,40 @@ def log_logistic(X, out=None):
     return out
 
 
+def softmax(X, copy=True):
+    """
+    Calculate the softmax function.
+
+    The softmax function is calculated by
+    np.exp(X) / np.sum(np.exp(X), axis=1)
+
+    This will cause overflow when large values are exponentiated.
+    Hence the largest value in each row is subtracted from each data
+    point to prevent this.
+
+    Parameters
+    ----------
+    X: array-like, shape (M, N)
+        Argument to the logistic function
+
+    copy: bool, optional
+        Copy X or not.
+
+    Returns
+    -------
+    out: array, shape (M, N)
+        Softmax function evaluated at every point in x
+    """
+    if copy:
+        X = np.copy(X)
+    max_prob = np.max(X, axis=1).reshape((-1, 1))
+    X -= max_prob
+    np.exp(X, X)
+    sum_prob = np.sum(X, axis=1).reshape((-1, 1))
+    X /= sum_prob
+    return X
+
+
 def safe_min(X):
     """Returns the minimum value of a dense or a CSR/CSC matrix.
 
diff --git a/sklearn/utils/tests/test_extmath.py b/sklearn/utils/tests/test_extmath.py
index b3f91b29568..0ca5dd1e07c 100644
--- a/sklearn/utils/tests/test_extmath.py
+++ b/sklearn/utils/tests/test_extmath.py
@@ -29,6 +29,7 @@ from sklearn.utils.extmath import fast_dot, _fast_dot
 from sklearn.utils.extmath import svd_flip
 from sklearn.utils.extmath import _batch_mean_variance_update
 from sklearn.utils.extmath import _deterministic_vector_sign_flip
+from sklearn.utils.extmath import softmax
 from sklearn.datasets.samples_generator import make_low_rank_matrix
 
 
@@ -465,3 +466,11 @@ def test_vector_sign_flip():
     assert_array_equal(max_abs_rows, max_rows)
     signs = np.sign(data[range(data.shape[0]), max_abs_rows])
     assert_array_equal(data, data_flipped * signs[:, np.newaxis])
+
+
+def test_softmax():
+    rng = np.random.RandomState(0)
+    X = rng.randn(3, 5)
+    exp_X = np.exp(X)
+    sum_exp_X = np.sum(exp_X, axis=1).reshape((-1, 1))
+    assert_array_almost_equal(softmax(X), exp_X / sum_exp_X)