diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1119845d38..3daac3a6b0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -27,6 +27,7 @@
 - PR #3112: Speed test_array
 - PR #3111: Adding Cython to Code Coverage
 - PR #3129:  Update notebooks README
+- PR #3135: Add QuasiNewton tests
 - PR #3040: Improved Array Conversion with CumlArrayDescriptor and Decorators
 - PR #3134: Improving the Deprecation Message Formatting in Documentation
 - PR #3137: Reorganize Pytest Config and Add Quick Run Option
diff --git a/python/cuml/test/test_qn.py b/python/cuml/test/test_qn.py
index 0e88da4393..e90728d7d5 100644
--- a/python/cuml/test/test_qn.py
+++ b/python/cuml/test/test_qn.py
@@ -18,6 +18,9 @@
 import cupy as cp
 
 from cuml.solvers import QN as cuQN
+from cuml.preprocessing.model_selection import train_test_split
+from cuml.datasets.classification import make_classification
+from cuml.metrics import accuracy_score
 
 
 # todo: add util functions to better compare against precomputed solutions
@@ -34,22 +37,37 @@ def test_qn(loss, dtype, penalty, l1_strength, l2_strength, fit_intercept):
 
     tol = 1e-6
 
-    X = np.array(precomputed_X, dtype=dtype)
-
-    if loss == 'sigmoid':
-        y = np.array(precomputed_y_log, dtype=dtype)
-    else:
-        y = np.array(precomputed_y_multi, dtype=dtype)
-
     qn = cuQN(loss=loss, fit_intercept=fit_intercept, l1_strength=l1_strength,
               l2_strength=l2_strength, tol=1e-8, output_type="cupy")
 
-    qn.fit(X, y)
+    if loss == 'softmax':
+        X, y = make_classification(n_samples=5000,
+                                   n_informative=10,
+                                   n_features=20,
+                                   n_classes=4,
+                                   dtype=dtype)
+
+        X_train, X_test, y_train, y_test = train_test_split(X.astype(dtype),
+                                                            y.astype(dtype),
+                                                            stratify=True)
+        most_class = cp.unique(y)[cp.argmax(cp.bincount(y))]
+
+        baseline_preds = cp.array([most_class] * y_test.shape[0], dtype=dtype)
+        baseline_score = accuracy_score(y_test, baseline_preds)
 
-    print(qn.objective)
-    print(qn.coef_)
+        y_pred = qn.fit(X_train, y_train).predict(X_test)
+        cuml_score = accuracy_score(y_test, y_pred)
+
+        assert(cuml_score > baseline_score)
+        assert(cuml_score >= 0.50)
+
+    elif loss == 'sigmoid':
+        X = np.array(precomputed_X, dtype=dtype)
+        y = np.array(precomputed_y_log, dtype=dtype)
+        qn.fit(X, y)
+        print(qn.objective)
+        print(qn.coef_)
 
-    if loss == 'sigmoid':
         if penalty == 'none' and l1_strength == 0.0 and l2_strength == 0.0:
             if fit_intercept:
                 assert (qn.objective - 0.40263831615448) < tol
@@ -198,10 +216,6 @@ def test_qn(loss, dtype, penalty, l1_strength, l2_strength, fit_intercept):
 
                 print()
 
-    elif loss == 'softmax':
-        pytest.skip("Better initial conditions for softmax tests are "
-                    "in progress.")
-
     # todo add tests for softmax dtype=np.float64
     # elasticnet for this points converged to different solution
     # if loss == 'softmax':