Repository: scikit-learn/scikit-learn
Pull request: API Accept 'auto' option for 'dual' parameter in LinearSVC and LinearSVR
Fixes #24651
closes #24731

<!--
copilot:summary
-->
### <samp>🤖 Generated by Copilot at cd87daa</samp>

Set `dual` to "auto" for `LinearSVC` and `LinearSVR` in some tests to avoid future warnings. This is a temporary fix for compatibility with the upcoming version of scikit-learn.

Problem statement:
Default value for solving an SVM in primal and dual should be determined automatically
### Describe the workflow you want to enable

Currently the default value for LinearSVC and LinearSVR with regards to solving the problem in dual is True. This means that for users that are unfamiliar with support-vector machines and do not know the importance of this parameter they could be using te framework inefficiently. 

I would propose for the default value to be None and a simple check to occur. If n_samples > n_features it's set to dual and otherwise it is set to primal. 

### Describe your proposed solution

This is a quick and easy fix and would only require a few lines of code in the fit method:

```
        if self.dual != None:
            pass
        else:
            if X.shape[0] < X.shape[1]:
                self.dual = True
            else:
                self.dual = False
```


### Describe alternatives you've considered, if relevant

_No response_

### Additional context

_No response_

Hints from issue discussion:
@agramfort I don't know enough about the potential side-effect of solving the SVM problem in the primal. Would having an `"auto"` option indeed be better to switch automatically from dual to primal depending on the shape of `X`?
indeed. Dual should be faster when n_samples < n_features. This would deserve a tiny benchmark though.
I would like to be involved in the benchmark and the eventual PR if it comes to that @agramfort.

Do I just simulate and arbitrary dataset with sklearn that has n > d using the current default (dual) and post the results in the comments?

I also expect that running the benchmark with dimensions 100000x20 is trivial in primal and will result in a memory issue on most machines in dual.
I think if you have the time then it would be great to contribute a benchmark and the PR @ChidiRnweke 
Summary of the benchmarks that I ran:

To put it shortly, for users that don't change the default parameters this is a big improvement for some datasets.

![barplot primal dual](https://user-images.githubusercontent.com/88949660/197419436-0b2d1f20-572b-4bb1-a327-b76dd11d1b01.png)

EDIT: Added a barplot instead of raw data to make the results readable.


It would be easier to have a bar plot to quickly compare the results :)
Instead of comparing `dual=True` and `dual="auto"`, I think that we want to compare `dual=True/False` and check that for `n_samples < n_features`, `dual=True` will always lead to a faster `fit` time.
I edited my comment to incorporate the feedback in terms of the benchmark and having a plot. Will be working on the changes on the adjustments to the pull request somewhere next week.
Is it stale?
It looks so. Do you want to take over @glevv ?

Commit message: lint

Patch to review:
diff --git a/doc/whats_new/v1.3.rst b/doc/whats_new/v1.3.rst
index e47c74a54edd6..480ec9e740031 100644
--- a/doc/whats_new/v1.3.rst
+++ b/doc/whats_new/v1.3.rst
@@ -443,6 +443,12 @@ Changelog
   The `sample_interval_` attribute is deprecated and will be removed in 1.5.
   :pr:`25190` by :user:`Vincent Maladière <Vincent-Maladiere>`.
 
+:mod:`sklearn.svm`
+..................
+
+    |Enhancement| `dual` parameter now accepts `auto` option for :class:`svm.LinearSVC` and :class:`svm.LinearSVR`.
+    :pr:`26093` by :user:`Gleb Levitski <glevv>`.
+
 :mod:`sklearn.tree`
 ...................
 
diff --git a/sklearn/calibration.py b/sklearn/calibration.py
index 36394a925309e..3eeb1b9fa2bdb 100644
--- a/sklearn/calibration.py
+++ b/sklearn/calibration.py
@@ -327,7 +327,7 @@ def fit(self, X, y, sample_weight=None, **fit_params):
         if estimator is None:
             # we want all classifiers that don't expose a random_state
             # to be deterministic (and we don't want to expose this one).
-            estimator = LinearSVC(random_state=0)
+            estimator = LinearSVC(dual="auto", random_state=0)
 
         self.calibrated_classifiers_ = []
         if self.cv == "prefit":
diff --git a/sklearn/ensemble/_stacking.py b/sklearn/ensemble/_stacking.py
index 81f8ef9c087a3..10f7a606f20c9 100644
--- a/sklearn/ensemble/_stacking.py
+++ b/sklearn/ensemble/_stacking.py
@@ -544,7 +544,7 @@ class StackingClassifier(ClassifierMixin, _BaseStacking):
     >>> estimators = [
     ...     ('rf', RandomForestClassifier(n_estimators=10, random_state=42)),
     ...     ('svr', make_pipeline(StandardScaler(),
-    ...                           LinearSVC(random_state=42)))
+    ...                           LinearSVC(dual="auto", random_state=42)))
     ... ]
     >>> clf = StackingClassifier(
     ...     estimators=estimators, final_estimator=LogisticRegression()
@@ -887,7 +887,7 @@ class StackingRegressor(RegressorMixin, _BaseStacking):
     >>> X, y = load_diabetes(return_X_y=True)
     >>> estimators = [
     ...     ('lr', RidgeCV()),
-    ...     ('svr', LinearSVR(random_state=42))
+    ...     ('svr', LinearSVR(dual="auto", random_state=42))
     ... ]
     >>> reg = StackingRegressor(
     ...     estimators=estimators,
diff --git a/sklearn/ensemble/tests/test_common.py b/sklearn/ensemble/tests/test_common.py
index 8dbbf0fae67f7..5bafe08881ae9 100644
--- a/sklearn/ensemble/tests/test_common.py
+++ b/sklearn/ensemble/tests/test_common.py
@@ -30,7 +30,7 @@
             StackingClassifier(
                 estimators=[
                     ("lr", LogisticRegression()),
-                    ("svm", LinearSVC()),
+                    ("svm", LinearSVC(dual="auto")),
                     ("rf", RandomForestClassifier(n_estimators=5, max_depth=3)),
                 ],
                 cv=2,
@@ -41,7 +41,7 @@
             VotingClassifier(
                 estimators=[
                     ("lr", LogisticRegression()),
-                    ("svm", LinearSVC()),
+                    ("svm", LinearSVC(dual="auto")),
                     ("rf", RandomForestClassifier(n_estimators=5, max_depth=3)),
                 ]
             ),
@@ -51,7 +51,7 @@
             StackingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR()),
+                    ("svm", LinearSVR(dual="auto")),
                     ("rf", RandomForestRegressor(n_estimators=5, max_depth=3)),
                 ],
                 cv=2,
@@ -62,7 +62,7 @@
             VotingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR()),
+                    ("svm", LinearSVR(dual="auto")),
                     ("rf", RandomForestRegressor(n_estimators=5, max_depth=3)),
                 ]
             ),
diff --git a/sklearn/ensemble/tests/test_stacking.py b/sklearn/ensemble/tests/test_stacking.py
index 956997156b455..2c04171fcd0f4 100644
--- a/sklearn/ensemble/tests/test_stacking.py
+++ b/sklearn/ensemble/tests/test_stacking.py
@@ -76,7 +76,7 @@ def test_stacking_classifier_iris(cv, final_estimator, passthrough):
     X_train, X_test, y_train, y_test = train_test_split(
         scale(X_iris), y_iris, stratify=y_iris, random_state=42
     )
-    estimators = [("lr", LogisticRegression()), ("svc", LinearSVC())]
+    estimators = [("lr", LogisticRegression()), ("svc", LinearSVC(dual="auto"))]
     clf = StackingClassifier(
         estimators=estimators,
         final_estimator=final_estimator,
@@ -128,7 +128,7 @@ def test_stacking_classifier_drop_column_binary_classification():
     assert X_trans.shape[1] == 2
 
     # LinearSVC does not implement 'predict_proba' and will not drop one column
-    estimators = [("lr", LogisticRegression()), ("svc", LinearSVC())]
+    estimators = [("lr", LogisticRegression()), ("svc", LinearSVC(dual="auto"))]
     clf.set_params(estimators=estimators)
 
     clf.fit(X_train, y_train)
@@ -142,10 +142,12 @@ def test_stacking_classifier_drop_estimator():
     X_train, X_test, y_train, _ = train_test_split(
         scale(X_iris), y_iris, stratify=y_iris, random_state=42
     )
-    estimators = [("lr", "drop"), ("svc", LinearSVC(random_state=0))]
+    estimators = [("lr", "drop"), ("svc", LinearSVC(dual="auto", random_state=0))]
     rf = RandomForestClassifier(n_estimators=10, random_state=42)
     clf = StackingClassifier(
-        estimators=[("svc", LinearSVC(random_state=0))], final_estimator=rf, cv=5
+        estimators=[("svc", LinearSVC(dual="auto", random_state=0))],
+        final_estimator=rf,
+        cv=5,
     )
     clf_drop = StackingClassifier(estimators=estimators, final_estimator=rf, cv=5)
 
@@ -162,10 +164,12 @@ def test_stacking_regressor_drop_estimator():
     X_train, X_test, y_train, _ = train_test_split(
         scale(X_diabetes), y_diabetes, random_state=42
     )
-    estimators = [("lr", "drop"), ("svr", LinearSVR(random_state=0))]
+    estimators = [("lr", "drop"), ("svr", LinearSVR(dual="auto", random_state=0))]
     rf = RandomForestRegressor(n_estimators=10, random_state=42)
     reg = StackingRegressor(
-        estimators=[("svr", LinearSVR(random_state=0))], final_estimator=rf, cv=5
+        estimators=[("svr", LinearSVR(dual="auto", random_state=0))],
+        final_estimator=rf,
+        cv=5,
     )
     reg_drop = StackingRegressor(estimators=estimators, final_estimator=rf, cv=5)
 
@@ -191,7 +195,7 @@ def test_stacking_regressor_diabetes(cv, final_estimator, predict_params, passth
     X_train, X_test, y_train, _ = train_test_split(
         scale(X_diabetes), y_diabetes, random_state=42
     )
-    estimators = [("lr", LinearRegression()), ("svr", LinearSVR())]
+    estimators = [("lr", LinearRegression()), ("svr", LinearSVR(dual="auto"))]
     reg = StackingRegressor(
         estimators=estimators,
         final_estimator=final_estimator,
@@ -227,7 +231,7 @@ def test_stacking_regressor_sparse_passthrough(fmt):
     X_train, X_test, y_train, _ = train_test_split(
         sparse.coo_matrix(scale(X_diabetes)).asformat(fmt), y_diabetes, random_state=42
     )
-    estimators = [("lr", LinearRegression()), ("svr", LinearSVR())]
+    estimators = [("lr", LinearRegression()), ("svr", LinearSVR(dual="auto"))]
     rf = RandomForestRegressor(n_estimators=10, random_state=42)
     clf = StackingRegressor(
         estimators=estimators, final_estimator=rf, cv=5, passthrough=True
@@ -245,7 +249,7 @@ def test_stacking_classifier_sparse_passthrough(fmt):
     X_train, X_test, y_train, _ = train_test_split(
         sparse.coo_matrix(scale(X_iris)).asformat(fmt), y_iris, random_state=42
     )
-    estimators = [("lr", LogisticRegression()), ("svc", LinearSVC())]
+    estimators = [("lr", LogisticRegression()), ("svc", LinearSVC(dual="auto"))]
     rf = RandomForestClassifier(n_estimators=10, random_state=42)
     clf = StackingClassifier(
         estimators=estimators, final_estimator=rf, cv=5, passthrough=True
@@ -318,7 +322,7 @@ def fit(self, X, y):
             {
                 "estimators": [
                     ("lr", LogisticRegression()),
-                    ("cor", LinearSVC(max_iter=50_000)),
+                    ("cor", LinearSVC(dual="auto", max_iter=50_000)),
                 ],
                 "final_estimator": NoWeightClassifier(),
             },
@@ -346,7 +350,10 @@ def test_stacking_classifier_error(y, params, type_err, msg_err):
         (
             y_diabetes,
             {
-                "estimators": [("lr", LinearRegression()), ("cor", LinearSVR())],
+                "estimators": [
+                    ("lr", LinearRegression()),
+                    ("cor", LinearSVR(dual="auto")),
+                ],
                 "final_estimator": NoWeightRegressor(),
             },
             TypeError,
@@ -367,7 +374,7 @@ def test_stacking_regressor_error(y, params, type_err, msg_err):
             StackingClassifier(
                 estimators=[
                     ("lr", LogisticRegression(random_state=0)),
-                    ("svm", LinearSVC(random_state=0)),
+                    ("svm", LinearSVC(dual="auto", random_state=0)),
                 ]
             ),
             X_iris[:100],
@@ -377,7 +384,7 @@ def test_stacking_regressor_error(y, params, type_err, msg_err):
             StackingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR(random_state=0)),
+                    ("svm", LinearSVR(dual="auto", random_state=0)),
                 ]
             ),
             X_diabetes,
@@ -411,7 +418,7 @@ def test_stacking_classifier_stratify_default():
     clf = StackingClassifier(
         estimators=[
             ("lr", LogisticRegression(max_iter=10_000)),
-            ("svm", LinearSVC(max_iter=10_000)),
+            ("svm", LinearSVC(dual="auto", max_iter=10_000)),
         ]
     )
     # since iris is not shuffled, a simple k-fold would not contain the
@@ -426,7 +433,7 @@ def test_stacking_classifier_stratify_default():
             StackingClassifier(
                 estimators=[
                     ("lr", LogisticRegression()),
-                    ("svm", LinearSVC(random_state=42)),
+                    ("svm", LinearSVC(dual="auto", random_state=42)),
                 ],
                 final_estimator=LogisticRegression(),
                 cv=KFold(shuffle=True, random_state=42),
@@ -437,7 +444,7 @@ def test_stacking_classifier_stratify_default():
             StackingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR(random_state=42)),
+                    ("svm", LinearSVR(dual="auto", random_state=42)),
                 ],
                 final_estimator=LinearRegression(),
                 cv=KFold(shuffle=True, random_state=42),
@@ -494,7 +501,7 @@ def test_stacking_classifier_sample_weight_fit_param():
             StackingClassifier(
                 estimators=[
                     ("lr", LogisticRegression()),
-                    ("svm", LinearSVC(random_state=42)),
+                    ("svm", LinearSVC(dual="auto", random_state=42)),
                 ],
                 final_estimator=LogisticRegression(),
             ),
@@ -504,7 +511,7 @@ def test_stacking_classifier_sample_weight_fit_param():
             StackingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR(random_state=42)),
+                    ("svm", LinearSVR(dual="auto", random_state=42)),
                 ],
                 final_estimator=LinearRegression(),
             ),
@@ -610,7 +617,7 @@ def test_stacking_prefit(Stacker, Estimator, stack_method, final_estimator, X, y
             StackingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR()),
+                    ("svm", LinearSVR(dual="auto")),
                 ],
                 cv="prefit",
             ),
@@ -776,7 +783,7 @@ def test_stacking_classifier_multilabel_auto_predict(stack_method, passthrough):
             StackingClassifier(
                 estimators=[
                     ("lr", LogisticRegression(random_state=0)),
-                    ("svm", LinearSVC(random_state=0)),
+                    ("svm", LinearSVC(dual="auto", random_state=0)),
                 ]
             ),
             iris.feature_names,
@@ -796,7 +803,7 @@ def test_stacking_classifier_multilabel_auto_predict(stack_method, passthrough):
                 estimators=[
                     ("lr", LogisticRegression(random_state=0)),
                     ("other", "drop"),
-                    ("svm", LinearSVC(random_state=0)),
+                    ("svm", LinearSVC(dual="auto", random_state=0)),
                 ]
             ),
             iris.feature_names,
@@ -811,7 +818,7 @@ def test_stacking_classifier_multilabel_auto_predict(stack_method, passthrough):
             StackingRegressor(
                 estimators=[
                     ("lr", LinearRegression()),
-                    ("svm", LinearSVR(random_state=0)),
+                    ("svm", LinearSVR(dual="auto", random_state=0)),
                 ]
             ),
             diabetes.feature_names,
diff --git a/sklearn/feature_extraction/tests/test_text.py b/sklearn/feature_extraction/tests/test_text.py
index d14a86d8f675e..80a42aaea5af0 100644
--- a/sklearn/feature_extraction/tests/test_text.py
+++ b/sklearn/feature_extraction/tests/test_text.py
@@ -934,7 +934,7 @@ def test_count_vectorizer_pipeline_grid_selection():
         data, target, test_size=0.2, random_state=0
     )
 
-    pipeline = Pipeline([("vect", CountVectorizer()), ("svc", LinearSVC())])
+    pipeline = Pipeline([("vect", CountVectorizer()), ("svc", LinearSVC(dual="auto"))])
 
     parameters = {
         "vect__ngram_range": [(1, 1), (1, 2)],
@@ -970,7 +970,7 @@ def test_vectorizer_pipeline_grid_selection():
         data, target, test_size=0.1, random_state=0
     )
 
-    pipeline = Pipeline([("vect", TfidfVectorizer()), ("svc", LinearSVC())])
+    pipeline = Pipeline([("vect", TfidfVectorizer()), ("svc", LinearSVC(dual="auto"))])
 
     parameters = {
         "vect__ngram_range": [(1, 1), (1, 2)],
@@ -1004,7 +1004,7 @@ def test_vectorizer_pipeline_cross_validation():
     # label junk food as -1, the others as +1
     target = [-1] * len(JUNK_FOOD_DOCS) + [1] * len(NOTJUNK_FOOD_DOCS)
 
-    pipeline = Pipeline([("vect", TfidfVectorizer()), ("svc", LinearSVC())])
+    pipeline = Pipeline([("vect", TfidfVectorizer()), ("svc", LinearSVC(dual="auto"))])
 
     cv_scores = cross_val_score(pipeline, data, target, cv=3)
     assert_array_equal(cv_scores, [1.0, 1.0, 1.0])
diff --git a/sklearn/feature_selection/tests/test_from_model.py b/sklearn/feature_selection/tests/test_from_model.py
index 8f4f0cdcc9308..7b408201bc7f5 100644
--- a/sklearn/feature_selection/tests/test_from_model.py
+++ b/sklearn/feature_selection/tests/test_from_model.py
@@ -405,7 +405,7 @@ def test_partial_fit():
 
 
 def test_calling_fit_reinitializes():
-    est = LinearSVC(random_state=0)
+    est = LinearSVC(dual="auto", random_state=0)
     transformer = SelectFromModel(estimator=est)
     transformer.fit(data, y)
     transformer.set_params(estimator__C=100)
diff --git a/sklearn/feature_selection/tests/test_rfe.py b/sklearn/feature_selection/tests/test_rfe.py
index fa7d5e826ad5f..fa7aeea19be6c 100644
--- a/sklearn/feature_selection/tests/test_rfe.py
+++ b/sklearn/feature_selection/tests/test_rfe.py
@@ -466,7 +466,7 @@ def test_rfe_wrapped_estimator(importance_getter, selector, expected_n_features)
     # Non-regression test for
     # https://github.com/scikit-learn/scikit-learn/issues/15312
     X, y = make_friedman1(n_samples=50, n_features=10, random_state=0)
-    estimator = LinearSVR(random_state=0)
+    estimator = LinearSVR(dual="auto", random_state=0)
 
     log_estimator = TransformedTargetRegressor(
         regressor=estimator, func=np.log, inverse_func=np.exp
@@ -488,7 +488,7 @@ def test_rfe_wrapped_estimator(importance_getter, selector, expected_n_features)
 @pytest.mark.parametrize("Selector", [RFE, RFECV])
 def test_rfe_importance_getter_validation(importance_getter, err_type, Selector):
     X, y = make_friedman1(n_samples=50, n_features=10, random_state=42)
-    estimator = LinearSVR()
+    estimator = LinearSVR(dual="auto")
     log_estimator = TransformedTargetRegressor(
         regressor=estimator, func=np.log, inverse_func=np.exp
     )
diff --git a/sklearn/kernel_approximation.py b/sklearn/kernel_approximation.py
index b79b4dcad457f..faa098e634937 100644
--- a/sklearn/kernel_approximation.py
+++ b/sklearn/kernel_approximation.py
@@ -965,13 +965,13 @@ class Nystroem(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator)
     >>> from sklearn.kernel_approximation import Nystroem
     >>> X, y = datasets.load_digits(n_class=9, return_X_y=True)
     >>> data = X / 16.
-    >>> clf = svm.LinearSVC()
+    >>> clf = svm.LinearSVC(dual="auto")
     >>> feature_map_nystroem = Nystroem(gamma=.2,
     ...                                 random_state=1,
     ...                                 n_components=300)
     >>> data_transformed = feature_map_nystroem.fit_transform(data)
     >>> clf.fit(data_transformed, y)
-    LinearSVC()
+    LinearSVC(dual='auto')
     >>> clf.score(data_transformed, y)
     0.9987...
     """
diff --git a/sklearn/metrics/_classification.py b/sklearn/metrics/_classification.py
index 5beb13df8069c..695bd3a74328f 100644
--- a/sklearn/metrics/_classification.py
+++ b/sklearn/metrics/_classification.py
@@ -2973,9 +2973,9 @@ def hinge_loss(y_true, pred_decision, *, labels=None, sample_weight=None):
     >>> from sklearn.metrics import hinge_loss
     >>> X = [[0], [1]]
     >>> y = [-1, 1]
-    >>> est = svm.LinearSVC(random_state=0)
+    >>> est = svm.LinearSVC(dual="auto", random_state=0)
     >>> est.fit(X, y)
-    LinearSVC(random_state=0)
+    LinearSVC(dual='auto', random_state=0)
     >>> pred_decision = est.decision_function([[-2], [3], [0.5]])
     >>> pred_decision
     array([-2.18...,  2.36...,  0.09...])
@@ -2988,9 +2988,9 @@ def hinge_loss(y_true, pred_decision, *, labels=None, sample_weight=None):
     >>> X = np.array([[0], [1], [2], [3]])
     >>> Y = np.array([0, 1, 2, 3])
     >>> labels = np.array([0, 1, 2, 3])
-    >>> est = svm.LinearSVC()
+    >>> est = svm.LinearSVC(dual="auto")
     >>> est.fit(X, Y)
-    LinearSVC()
+    LinearSVC(dual='auto')
     >>> pred_decision = est.decision_function([[-1], [2], [3]])
     >>> y_true = [0, 2, 3]
     >>> hinge_loss(y_true, pred_decision, labels=labels)
diff --git a/sklearn/metrics/tests/test_score_objects.py b/sklearn/metrics/tests/test_score_objects.py
index 7f3e804f68d46..bbdcc71350d9a 100644
--- a/sklearn/metrics/tests/test_score_objects.py
+++ b/sklearn/metrics/tests/test_score_objects.py
@@ -291,7 +291,7 @@ def test_check_scoring_and_check_multimetric_scoring(scoring):
     # To make sure the check_scoring is correctly applied to the constituent
     # scorers
 
-    estimator = LinearSVC(random_state=0)
+    estimator = LinearSVC(dual="auto", random_state=0)
     estimator.fit([[1], [2], [3]], [1, 1, 0])
 
     scorers = _check_multimetric_scoring(estimator, scoring)
@@ -353,11 +353,11 @@ def test_check_scoring_gridsearchcv():
     # test that check_scoring works on GridSearchCV and pipeline.
     # slightly redundant non-regression test.
 
-    grid = GridSearchCV(LinearSVC(), param_grid={"C": [0.1, 1]}, cv=3)
+    grid = GridSearchCV(LinearSVC(dual="auto"), param_grid={"C": [0.1, 1]}, cv=3)
     scorer = check_scoring(grid, scoring="f1")
     assert isinstance(scorer, _PredictScorer)
 
-    pipe = make_pipeline(LinearSVC())
+    pipe = make_pipeline(LinearSVC(dual="auto"))
     scorer = check_scoring(pipe, scoring="f1")
     assert isinstance(scorer, _PredictScorer)
 
@@ -405,7 +405,7 @@ def test_classification_binary_scores(scorer_name, metric):
     # binary classification.
     X, y = make_blobs(random_state=0, centers=2)
     X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     clf.fit(X_train, y_train)
 
     score = get_scorer(scorer_name)(clf, X_test, y_test)
@@ -455,7 +455,7 @@ def test_custom_scorer_pickling():
     # test that custom scorer can be pickled
     X, y = make_blobs(random_state=0, centers=2)
     X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     clf.fit(X_train, y_train)
 
     scorer = make_scorer(fbeta_score, beta=2)
@@ -566,7 +566,7 @@ def test_thresholded_scorers_multilabel_indicator_data():
     assert_almost_equal(score1, score2)
 
     # Multilabel decision function
-    clf = OneVsRestClassifier(LinearSVC(random_state=0))
+    clf = OneVsRestClassifier(LinearSVC(dual="auto", random_state=0))
     clf.fit(X_train, y_train)
     score1 = get_scorer("roc_auc")(clf, X_test, y_test)
     score2 = roc_auc_score(y_test, clf.decision_function(X_test))
diff --git a/sklearn/model_selection/tests/test_search.py b/sklearn/model_selection/tests/test_search.py
index 0aef56da2d3c8..a8ab33a334c85 100644
--- a/sklearn/model_selection/tests/test_search.py
+++ b/sklearn/model_selection/tests/test_search.py
@@ -249,10 +249,10 @@ def test_SearchCV_with_fit_params(SearchCV):
 @ignore_warnings
 def test_grid_search_no_score():
     # Test grid-search on classifier that has no score function.
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     X, y = make_blobs(random_state=0, centers=2)
     Cs = [0.1, 1, 10]
-    clf_no_score = LinearSVCNoScore(random_state=0)
+    clf_no_score = LinearSVCNoScore(dual="auto", random_state=0)
     grid_search = GridSearchCV(clf, {"C": Cs}, scoring="accuracy")
     grid_search.fit(X, y)
 
@@ -273,13 +273,13 @@ def test_grid_search_no_score():
 
 def test_grid_search_score_method():
     X, y = make_classification(n_samples=100, n_classes=2, flip_y=0.2, random_state=0)
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     grid = {"C": [0.1]}
 
     search_no_scoring = GridSearchCV(clf, grid, scoring=None).fit(X, y)
     search_accuracy = GridSearchCV(clf, grid, scoring="accuracy").fit(X, y)
     search_no_score_method_auc = GridSearchCV(
-        LinearSVCNoScore(), grid, scoring="roc_auc"
+        LinearSVCNoScore(dual="auto"), grid, scoring="roc_auc"
     ).fit(X, y)
     search_auc = GridSearchCV(clf, grid, scoring="roc_auc").fit(X, y)
 
@@ -307,7 +307,7 @@ def test_grid_search_groups():
     X, y = make_classification(n_samples=15, n_classes=2, random_state=0)
     groups = rng.randint(0, 3, 15)
 
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     grid = {"C": [1]}
 
     group_cvs = [
@@ -336,7 +336,7 @@ def test_classes__property():
     y = np.array([0] * 5 + [1] * 5)
     Cs = [0.1, 1, 10]
 
-    grid_search = GridSearchCV(LinearSVC(random_state=0), {"C": Cs})
+    grid_search = GridSearchCV(LinearSVC(dual="auto", random_state=0), {"C": Cs})
     grid_search.fit(X, y)
     assert_array_equal(grid_search.best_estimator_.classes_, grid_search.classes_)
 
@@ -346,11 +346,13 @@ def test_classes__property():
     assert not hasattr(grid_search, "classes_")
 
     # Test that the grid searcher has no classes_ attribute before it's fit
-    grid_search = GridSearchCV(LinearSVC(random_state=0), {"C": Cs})
+    grid_search = GridSearchCV(LinearSVC(dual="auto", random_state=0), {"C": Cs})
     assert not hasattr(grid_search, "classes_")
 
     # Test that the grid searcher has no classes_ attribute without a refit
-    grid_search = GridSearchCV(LinearSVC(random_state=0), {"C": Cs}, refit=False)
+    grid_search = GridSearchCV(
+        LinearSVC(dual="auto", random_state=0), {"C": Cs}, refit=False
+    )
     grid_search.fit(X, y)
     assert not hasattr(grid_search, "classes_")
 
@@ -410,7 +412,7 @@ def test_grid_search_error():
     # Test that grid search will capture errors on data with different length
     X_, y_ = make_classification(n_samples=200, n_features=100, random_state=0)
 
-    clf = LinearSVC()
+    clf = LinearSVC(dual="auto")
     cv = GridSearchCV(clf, {"C": [0.1, 1.0]})
     with pytest.raises(ValueError):
         cv.fit(X_[:180], y_)
@@ -483,14 +485,14 @@ def test_grid_search_sparse():
     # Test that grid search works with both dense and sparse matrices
     X_, y_ = make_classification(n_samples=200, n_features=100, random_state=0)
 
-    clf = LinearSVC()
+    clf = LinearSVC(dual="auto")
     cv = GridSearchCV(clf, {"C": [0.1, 1.0]})
     cv.fit(X_[:180], y_[:180])
     y_pred = cv.predict(X_[180:])
     C = cv.best_estimator_.C
 
     X_ = sp.csr_matrix(X_)
-    clf = LinearSVC()
+    clf = LinearSVC(dual="auto")
     cv = GridSearchCV(clf, {"C": [0.1, 1.0]})
     cv.fit(X_[:180].tocoo(), y_[:180])
     y_pred2 = cv.predict(X_[180:])
@@ -503,14 +505,14 @@ def test_grid_search_sparse():
 def test_grid_search_sparse_scoring():
     X_, y_ = make_classification(n_samples=200, n_features=100, random_state=0)
 
-    clf = LinearSVC()
+    clf = LinearSVC(dual="auto")
     cv = GridSearchCV(clf, {"C": [0.1, 1.0]}, scoring="f1")
     cv.fit(X_[:180], y_[:180])
     y_pred = cv.predict(X_[180:])
     C = cv.best_estimator_.C
 
     X_ = sp.csr_matrix(X_)
-    clf = LinearSVC()
+    clf = LinearSVC(dual="auto")
     cv = GridSearchCV(clf, {"C": [0.1, 1.0]}, scoring="f1")
     cv.fit(X_[:180], y_[:180])
     y_pred2 = cv.predict(X_[180:])
@@ -620,7 +622,7 @@ def refit_callable(cv_results):
         # clf.cv_results_.
         X, y = make_classification(n_samples=100, n_features=4, random_state=42)
         clf = GridSearchCV(
-            LinearSVC(random_state=42),
+            LinearSVC(dual="auto", random_state=42),
             {"C": [0.01, 0.1, 1]},
             scoring="precision",
             refit=True,
@@ -637,7 +639,7 @@ def refit_callable(cv_results):
 
     X, y = make_classification(n_samples=100, n_features=4, random_state=42)
     clf = GridSearchCV(
-        LinearSVC(random_state=42),
+        LinearSVC(dual="auto", random_state=42),
         {"C": [0.01, 0.1, 1]},
         scoring="precision",
         refit=refit_callable,
@@ -664,7 +666,7 @@ def refit_callable_invalid_type(cv_results):
     X, y = make_classification(n_samples=100, n_features=4, random_state=42)
 
     clf = GridSearchCV(
-        LinearSVC(random_state=42),
+        LinearSVC(dual="auto", random_state=42),
         {"C": [0.1, 1]},
         scoring="precision",
         refit=refit_callable_invalid_type,
@@ -690,7 +692,7 @@ def refit_callable_out_bound(cv_results):
     X, y = make_classification(n_samples=100, n_features=4, random_state=42)
 
     clf = search_cv(
-        LinearSVC(random_state=42),
+        LinearSVC(dual="auto", random_state=42),
         {"C": [0.1, 1]},
         scoring="precision",
         refit=refit_callable_out_bound,
@@ -716,7 +718,7 @@ def refit_callable(cv_results):
     X, y = make_classification(n_samples=100, n_features=4, random_state=42)
     scoring = {"Accuracy": make_scorer(accuracy_score), "prec": "precision"}
     clf = GridSearchCV(
-        LinearSVC(random_state=42),
+        LinearSVC(dual="auto", random_state=42),
         {"C": [0.01, 0.1, 1]},
         scoring=scoring,
         refit=refit_callable,
@@ -1365,7 +1367,7 @@ def test_search_cv_results_none_param():
 
 @ignore_warnings()
 def test_search_cv_timing():
-    svc = LinearSVC(random_state=0)
+    svc = LinearSVC(dual="auto", random_state=0)
 
     X = [
         [
@@ -1407,7 +1409,7 @@ def test_search_cv_timing():
 def test_grid_search_correct_score_results():
     # test that correct scores are used
     n_splits = 3
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     X, y = make_blobs(random_state=0, centers=2)
     Cs = [0.1, 1, 10]
     for score in ["f1", "roc_auc"]:
@@ -1746,7 +1748,7 @@ def test_stochastic_gradient_loss_param():
 def test_search_train_scores_set_to_false():
     X = np.arange(6).reshape(6, -1)
     y = [0, 0, 0, 1, 1, 1]
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
 
     gs = GridSearchCV(clf, param_grid={"C": [0.1, 0.2]}, cv=3)
     gs.fit(X, y)
@@ -1759,7 +1761,7 @@ def test_grid_search_cv_splits_consistency():
     X, y = make_classification(n_samples=n_samples, random_state=0)
 
     gs = GridSearchCV(
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         param_grid={"C": [0.1, 0.2, 0.3]},
         cv=OneTimeSplitter(n_splits=n_splits, n_samples=n_samples),
         return_train_score=True,
@@ -1767,7 +1769,7 @@ def test_grid_search_cv_splits_consistency():
     gs.fit(X, y)
 
     gs2 = GridSearchCV(
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         param_grid={"C": [0.1, 0.2, 0.3]},
         cv=KFold(n_splits=n_splits),
         return_train_score=True,
@@ -1780,7 +1782,7 @@ def test_grid_search_cv_splits_consistency():
         GeneratorType,
     )
     gs3 = GridSearchCV(
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         param_grid={"C": [0.1, 0.2, 0.3]},
         cv=KFold(n_splits=n_splits, shuffle=True, random_state=0).split(X, y),
         return_train_score=True,
@@ -1788,7 +1790,7 @@ def test_grid_search_cv_splits_consistency():
     gs3.fit(X, y)
 
     gs4 = GridSearchCV(
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         param_grid={"C": [0.1, 0.2, 0.3]},
         cv=KFold(n_splits=n_splits, shuffle=True, random_state=0),
         return_train_score=True,
@@ -1824,7 +1826,7 @@ def _pop_time_keys(cv_results):
 
     # Check consistency of folds across the parameters
     gs = GridSearchCV(
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         param_grid={"C": [0.1, 0.1, 0.2, 0.2]},
         cv=KFold(n_splits=n_splits, shuffle=True),
         return_train_score=True,
@@ -2042,7 +2044,7 @@ def custom_scorer(clf, X, y):
         return {"tn": cm[0, 0], "fp": cm[0, 1], "fn": cm[1, 0], "tp": cm[1, 1]}
 
     X, y = make_classification(n_samples=40, n_features=4, random_state=42)
-    est = LinearSVC(random_state=42)
+    est = LinearSVC(dual="auto", random_state=42)
     search = GridSearchCV(est, {"C": [0.1, 1]}, scoring=custom_scorer, refit="fp")
 
     search.fit(X, y)
@@ -2066,7 +2068,7 @@ def custom_scorer(est, X, y):
         }
 
     X, y = make_classification(n_samples=40, n_features=4, random_state=42)
-    est = LinearSVC(random_state=42)
+    est = LinearSVC(dual="auto", random_state=42)
     search_callable = GridSearchCV(
         est, {"C": [0.1, 1]}, scoring=custom_scorer, refit="recall"
     )
@@ -2089,7 +2091,7 @@ def custom_scorer(est, X, y):
         return recall_score(y, y_pred)
 
     X, y = make_classification(n_samples=40, n_features=4, random_state=42)
-    est = LinearSVC(random_state=42)
+    est = LinearSVC(dual="auto", random_state=42)
     search_callable = GridSearchCV(
         est, {"C": [0.1, 1]}, scoring=custom_scorer, refit=True
     )
@@ -2117,7 +2119,7 @@ def bad_scorer(est, X, y):
 
     X, y = make_classification(n_samples=40, n_features=4, random_state=42)
     clf = GridSearchCV(
-        LinearSVC(random_state=42),
+        LinearSVC(dual="auto", random_state=42),
         {"C": [0.1, 1]},
         scoring=bad_scorer,
         refit="good_name",
@@ -2399,7 +2401,7 @@ def test_search_cv_verbose_3(capsys, return_train_score):
     """Check that search cv with verbose>2 shows the score for single
     metrics. non-regression test for #19658."""
     X, y = make_classification(n_samples=100, n_classes=2, flip_y=0.2, random_state=0)
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     grid = {"C": [0.1]}
 
     GridSearchCV(
diff --git a/sklearn/model_selection/tests/test_successive_halving.py b/sklearn/model_selection/tests/test_successive_halving.py
index d0e3b9f78c918..035e20dc701d8 100644
--- a/sklearn/model_selection/tests/test_successive_halving.py
+++ b/sklearn/model_selection/tests/test_successive_halving.py
@@ -725,7 +725,7 @@ def test_groups_support(Est):
     X, y = make_classification(n_samples=50, n_classes=2, random_state=0)
     groups = rng.randint(0, 3, 50)
 
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     grid = {"C": [1]}
 
     group_cvs = [
diff --git a/sklearn/model_selection/tests/test_validation.py b/sklearn/model_selection/tests/test_validation.py
index 1ed42d555309a..f21c3eb326135 100644
--- a/sklearn/model_selection/tests/test_validation.py
+++ b/sklearn/model_selection/tests/test_validation.py
@@ -2383,7 +2383,7 @@ def custom_scorer(clf, X, y):
         return {"tn": cm[0, 0], "fp": cm[0, 1], "fn": cm[1, 0], "tp": cm[1, 1]}
 
     X, y = make_classification(n_samples=40, n_features=4, random_state=42)
-    est = LinearSVC(random_state=42)
+    est = LinearSVC(dual="auto", random_state=42)
     est.fit(X, y)
     cv_results = cross_validate(est, X, y, cv=5, scoring=custom_scorer)
 
diff --git a/sklearn/multiclass.py b/sklearn/multiclass.py
index f3f69a3c6bbf9..74684e608d3c1 100644
--- a/sklearn/multiclass.py
+++ b/sklearn/multiclass.py
@@ -641,7 +641,7 @@ class OneVsOneClassifier(MetaEstimatorMixin, ClassifierMixin, BaseEstimator):
     >>> X_train, X_test, y_train, y_test = train_test_split(
     ...     X, y, test_size=0.33, shuffle=True, random_state=0)
     >>> clf = OneVsOneClassifier(
-    ...     LinearSVC(random_state=0)).fit(X_train, y_train)
+    ...     LinearSVC(dual="auto", random_state=0)).fit(X_train, y_train)
     >>> clf.predict(X_test[:10])
     array([2, 1, 0, 2, 0, 2, 0, 1, 1, 1])
     """
diff --git a/sklearn/svm/_classes.py b/sklearn/svm/_classes.py
index 94ba958b5a905..24ae8a7650965 100644
--- a/sklearn/svm/_classes.py
+++ b/sklearn/svm/_classes.py
@@ -1,14 +1,35 @@
 from numbers import Integral, Real
+import warnings
 
 import numpy as np
 
-from ._base import _fit_liblinear, BaseSVC, BaseLibSVM
+from ._base import _fit_liblinear, _get_liblinear_solver_type, BaseSVC, BaseLibSVM
 from ..base import BaseEstimator, RegressorMixin, OutlierMixin
 from ..linear_model._base import LinearClassifierMixin, SparseCoefMixin, LinearModel
 from ..utils import deprecated
 from ..utils.validation import _num_samples
 from ..utils.multiclass import check_classification_targets
-from ..utils._param_validation import Interval, StrOptions
+from ..utils._param_validation import Interval, StrOptions, Hidden
+
+
+def _choose_dual_automatically(loss, penalty, multi_class, X):
+    """Choose dual parameter value wrt to loss, penalty and shape of the data"""
+    if X.shape[0] < X.shape[1]:
+        try:
+            # check if the combination of parameters supports
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, True)
+            return True
+        except ValueError:
+            return False
+    else:
+        try:
+            # check if the combination of parameters doesn't support
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, False)
+            return False
+        except ValueError:
+            return True
 
 
 class LinearSVC(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
@@ -37,9 +58,14 @@ class LinearSVC(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
         square of the hinge loss. The combination of ``penalty='l1'``
         and ``loss='hinge'`` is not supported.
 
-    dual : bool, default=True
+    dual : "auto" or bool, default=True
         Select the algorithm to either solve the dual or primal
         optimization problem. Prefer dual=False when n_samples > n_features.
+        `dual="auto"` will choose the best option considering `n_samples`,
+        `n_features`, `loss` and `penalty`.
+
+        .. versionchanged:: 1.3
+           The default value will change from `True` to `"auto"` in 1.5.
 
     tol : float, default=1e-4
         Tolerance for stopping criteria.
@@ -176,10 +202,10 @@ class LinearSVC(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
     >>> from sklearn.datasets import make_classification
     >>> X, y = make_classification(n_features=4, random_state=0)
     >>> clf = make_pipeline(StandardScaler(),
-    ...                     LinearSVC(random_state=0, tol=1e-5))
+    ...                     LinearSVC(dual="auto", random_state=0, tol=1e-5))
     >>> clf.fit(X, y)
     Pipeline(steps=[('standardscaler', StandardScaler()),
-                    ('linearsvc', LinearSVC(random_state=0, tol=1e-05))])
+                    ('linearsvc', LinearSVC(dual='auto', random_state=0, tol=1e-05))])
 
     >>> print(clf.named_steps['linearsvc'].coef_)
     [[0.141...   0.526... 0.679... 0.493...]]
@@ -193,7 +219,7 @@ class LinearSVC(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
     _parameter_constraints: dict = {
         "penalty": [StrOptions({"l1", "l2"})],
         "loss": [StrOptions({"hinge", "squared_hinge"})],
-        "dual": ["boolean"],
+        "dual": ["boolean", StrOptions({"auto"}), Hidden(StrOptions({"warn"}))],
         "tol": [Interval(Real, 0.0, None, closed="neither")],
         "C": [Interval(Real, 0.0, None, closed="neither")],
         "multi_class": [StrOptions({"ovr", "crammer_singer"})],
@@ -210,7 +236,7 @@ def __init__(
         penalty="l2",
         loss="squared_hinge",
         *,
-        dual=True,
+        dual="warn",
         tol=1e-4,
         C=1.0,
         multi_class="ovr",
@@ -271,6 +297,22 @@ def fit(self, X, y, sample_weight=None):
         check_classification_targets(y)
         self.classes_ = np.unique(y)
 
+        if self.dual == "auto":
+            self._dual = _choose_dual_automatically(
+                self.loss, self.penalty, self.multi_class, X
+            )
+        elif self.dual == "warn":
+            warnings.warn(
+                (
+                    "The default value of `dual` will change from `True` to `'auto'` in"
+                    " 1.5. Set the value of `dual` explicitly to suppress the warning."
+                ),
+                FutureWarning,
+            )
+            self._dual = True
+        else:
+            self._dual = self.dual
+
         self.coef_, self.intercept_, n_iter_ = _fit_liblinear(
             X,
             y,
@@ -279,7 +321,7 @@ def fit(self, X, y, sample_weight=None):
             self.intercept_scaling,
             self.class_weight,
             self.penalty,
-            self.dual,
+            self._dual,
             self.verbose,
             self.max_iter,
             self.tol,
@@ -362,9 +404,14 @@ class LinearSVR(RegressorMixin, LinearModel):
         To lessen the effect of regularization on synthetic feature weight
         (and therefore on the intercept) intercept_scaling has to be increased.
 
-    dual : bool, default=True
+    dual : "auto" or bool, default=True
         Select the algorithm to either solve the dual or primal
         optimization problem. Prefer dual=False when n_samples > n_features.
+        `dual="auto"` will choose the best option considering `n_samples`,
+        `n_features` and `loss`.
+
+        .. versionchanged:: 1.3
+           The default value will change from `True` to `"auto"` in 1.5.
 
     verbose : int, default=0
         Enable verbose output. Note that this setting takes advantage of a
@@ -429,10 +476,10 @@ class LinearSVR(RegressorMixin, LinearModel):
     >>> from sklearn.datasets import make_regression
     >>> X, y = make_regression(n_features=4, random_state=0)
     >>> regr = make_pipeline(StandardScaler(),
-    ...                      LinearSVR(random_state=0, tol=1e-5))
+    ...                      LinearSVR(dual="auto", random_state=0, tol=1e-5))
     >>> regr.fit(X, y)
     Pipeline(steps=[('standardscaler', StandardScaler()),
-                    ('linearsvr', LinearSVR(random_state=0, tol=1e-05))])
+                    ('linearsvr', LinearSVR(dual='auto', random_state=0, tol=1e-05))])
 
     >>> print(regr.named_steps['linearsvr'].coef_)
     [18.582... 27.023... 44.357... 64.522...]
@@ -449,7 +496,7 @@ class LinearSVR(RegressorMixin, LinearModel):
         "loss": [StrOptions({"epsilon_insensitive", "squared_epsilon_insensitive"})],
         "fit_intercept": ["boolean"],
         "intercept_scaling": [Interval(Real, 0, None, closed="neither")],
-        "dual": ["boolean"],
+        "dual": ["boolean", StrOptions({"auto"}), Hidden(StrOptions({"warn"}))],
         "verbose": ["verbose"],
         "random_state": ["random_state"],
         "max_iter": [Interval(Integral, 0, None, closed="left")],
@@ -464,7 +511,7 @@ def __init__(
         loss="epsilon_insensitive",
         fit_intercept=True,
         intercept_scaling=1.0,
-        dual=True,
+        dual="warn",
         verbose=0,
         random_state=None,
         max_iter=1000,
@@ -515,6 +562,21 @@ def fit(self, X, y, sample_weight=None):
             accept_large_sparse=False,
         )
         penalty = "l2"  # SVR only accepts l2 penalty
+
+        if self.dual == "auto":
+            self._dual = _choose_dual_automatically(self.loss, penalty, "ovr", X)
+        elif self.dual == "warn":
+            warnings.warn(
+                (
+                    "The default value of `dual` will change from `True` to `'auto'` in"
+                    " 1.5. Set the value of `dual` explicitly to suppress the warning."
+                ),
+                FutureWarning,
+            )
+            self._dual = True
+        else:
+            self._dual = self.dual
+
         self.coef_, self.intercept_, n_iter_ = _fit_liblinear(
             X,
             y,
@@ -523,7 +585,7 @@ def fit(self, X, y, sample_weight=None):
             self.intercept_scaling,
             None,
             penalty,
-            self.dual,
+            self._dual,
             self.verbose,
             self.max_iter,
             self.tol,
diff --git a/sklearn/svm/tests/test_sparse.py b/sklearn/svm/tests/test_sparse.py
index 9f0102f08fa54..97c63b0597c48 100644
--- a/sklearn/svm/tests/test_sparse.py
+++ b/sklearn/svm/tests/test_sparse.py
@@ -230,8 +230,8 @@ def test_error():
 
 def test_linearsvc():
     # Similar to test_SVC
-    clf = svm.LinearSVC(random_state=0).fit(X, Y)
-    sp_clf = svm.LinearSVC(random_state=0).fit(X_sp, Y)
+    clf = svm.LinearSVC(dual="auto", random_state=0).fit(X, Y)
+    sp_clf = svm.LinearSVC(dual="auto", random_state=0).fit(X_sp, Y)
 
     assert sp_clf.fit_intercept
 
@@ -250,8 +250,10 @@ def test_linearsvc():
 def test_linearsvc_iris():
     # Test the sparse LinearSVC with the iris dataset
 
-    sp_clf = svm.LinearSVC(random_state=0).fit(iris.data, iris.target)
-    clf = svm.LinearSVC(random_state=0).fit(iris.data.toarray(), iris.target)
+    sp_clf = svm.LinearSVC(dual="auto", random_state=0).fit(iris.data, iris.target)
+    clf = svm.LinearSVC(dual="auto", random_state=0).fit(
+        iris.data.toarray(), iris.target
+    )
 
     assert clf.fit_intercept == sp_clf.fit_intercept
 
@@ -282,7 +284,7 @@ def test_weight():
     X_ = sparse.csr_matrix(X_)
     for clf in (
         linear_model.LogisticRegression(),
-        svm.LinearSVC(random_state=0),
+        svm.LinearSVC(dual="auto", random_state=0),
         svm.SVC(),
     ):
         clf.set_params(class_weight={0: 5})
diff --git a/sklearn/svm/tests/test_svm.py b/sklearn/svm/tests/test_svm.py
index ca23360f9d462..0e4cbeca464e6 100644
--- a/sklearn/svm/tests/test_svm.py
+++ b/sklearn/svm/tests/test_svm.py
@@ -15,6 +15,7 @@
 from scipy import sparse
 from sklearn import svm, linear_model, datasets, metrics, base
 from sklearn.svm import LinearSVC, OneClassSVM, SVR, NuSVR, LinearSVR
+from sklearn.svm._classes import _choose_dual_automatically
 from sklearn.model_selection import train_test_split
 from sklearn.datasets import make_classification, make_blobs
 from sklearn.metrics import f1_score
@@ -30,6 +31,13 @@
 # mypy error: Module 'sklearn.svm' has no attribute '_libsvm'
 from sklearn.svm import _libsvm  # type: ignore
 
+# TODO(1.5): Remove
+msg = (
+    r"The default value of `dual` will change from `True` to `'auto'` in 1.5. "
+    r"Set the value of `dual` explicitly to suppress the warning.:FutureWarning"
+)
+pytestmark = pytest.mark.filterwarnings("ignore:" + msg)
+
 # toy sample
 X = [[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]]
 Y = [1, 1, 1, 2, 2, 2]
@@ -1390,3 +1398,25 @@ def test_svm_class_weights_deprecation(Klass):
     )
     with pytest.warns(FutureWarning, match=re.escape(msg)):
         getattr(clf, "class_weight_")
+
+
+# TODO(1.5): Remove
+@pytest.mark.parametrize("Estimator", [LinearSVR, LinearSVC])
+def test_dual_auto_deprecation_warning(Estimator):
+    svm = Estimator()
+    msg = (
+        "The default value of `dual` will change from `True` to `'auto'` in"
+        " 1.5. Set the value of `dual` explicitly to suppress the warning."
+    )
+    with pytest.warns(FutureWarning, match=re.escape(msg)):
+        svm.fit(X, Y)
+
+
+@pytest.mark.parametrize("loss", ["squared_hinge", "squared_epsilon_insensitive"])
+def test_dual_auto(loss):
+    # OvR, L2, N > M (6,2)
+    dual = _choose_dual_automatically(loss, "l2", "ovr", np.asarray(X))
+    assert dual is False
+    # OvR, L2, N < M (2,6)
+    dual = _choose_dual_automatically(loss, "l2", "ovr", np.asarray(X).T)
+    assert dual is True
diff --git a/sklearn/tests/test_calibration.py b/sklearn/tests/test_calibration.py
index 679ec6ed52c15..a19785a60c308 100644
--- a/sklearn/tests/test_calibration.py
+++ b/sklearn/tests/test_calibration.py
@@ -152,7 +152,7 @@ def test_sample_weight(data, method, ensemble):
     X_train, y_train, sw_train = X[:n_samples], y[:n_samples], sample_weight[:n_samples]
     X_test = X[n_samples:]
 
-    estimator = LinearSVC(random_state=42)
+    estimator = LinearSVC(dual="auto", random_state=42)
     calibrated_clf = CalibratedClassifierCV(estimator, method=method, ensemble=ensemble)
     calibrated_clf.fit(X_train, y_train, sample_weight=sw_train)
     probs_with_sw = calibrated_clf.predict_proba(X_test)
@@ -173,7 +173,7 @@ def test_parallel_execution(data, method, ensemble):
     X, y = data
     X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=42)
 
-    estimator = make_pipeline(StandardScaler(), LinearSVC(random_state=42))
+    estimator = make_pipeline(StandardScaler(), LinearSVC(dual="auto", random_state=42))
 
     cal_clf_parallel = CalibratedClassifierCV(
         estimator, method=method, n_jobs=2, ensemble=ensemble
@@ -202,7 +202,7 @@ def multiclass_brier(y_true, proba_pred, n_classes):
 
     # Test calibration for multiclass with classifier that implements
     # only decision function.
-    clf = LinearSVC(random_state=7)
+    clf = LinearSVC(dual="auto", random_state=7)
     X, y = make_blobs(
         n_samples=500, n_features=100, random_state=seed, centers=10, cluster_std=15.0
     )
@@ -333,7 +333,7 @@ def test_calibration_ensemble_false(data, method):
     # Test that `ensemble=False` is the same as using predictions from
     # `cross_val_predict` to train calibrator.
     X, y = data
-    clf = LinearSVC(random_state=7)
+    clf = LinearSVC(dual="auto", random_state=7)
 
     cal_clf = CalibratedClassifierCV(clf, method=method, cv=3, ensemble=False)
     cal_clf.fit(X, y)
@@ -422,7 +422,7 @@ def test_calibration_prob_sum(ensemble):
     # issue #7796
     num_classes = 2
     X, y = make_classification(n_samples=10, n_features=5, n_classes=num_classes)
-    clf = LinearSVC(C=1.0, random_state=7)
+    clf = LinearSVC(dual="auto", C=1.0, random_state=7)
     clf_prob = CalibratedClassifierCV(
         clf, method="sigmoid", cv=LeaveOneOut(), ensemble=ensemble
     )
@@ -440,7 +440,7 @@ def test_calibration_less_classes(ensemble):
     # class label
     X = np.random.randn(10, 5)
     y = np.arange(10)
-    clf = LinearSVC(C=1.0, random_state=7)
+    clf = LinearSVC(dual="auto", C=1.0, random_state=7)
     cal_clf = CalibratedClassifierCV(
         clf, method="sigmoid", cv=LeaveOneOut(), ensemble=ensemble
     )
@@ -535,8 +535,8 @@ def test_calibration_dict_pipeline(dict_data, dict_data_pipeline):
 @pytest.mark.parametrize(
     "clf, cv",
     [
-        pytest.param(LinearSVC(C=1), 2),
-        pytest.param(LinearSVC(C=1), "prefit"),
+        pytest.param(LinearSVC(dual="auto", C=1), 2),
+        pytest.param(LinearSVC(dual="auto", C=1), "prefit"),
     ],
 )
 def test_calibration_attributes(clf, cv):
@@ -560,7 +560,7 @@ def test_calibration_inconsistent_prefit_n_features_in():
     # Check that `n_features_in_` from prefit base estimator
     # is consistent with training set
     X, y = make_classification(n_samples=10, n_features=5, n_classes=2, random_state=7)
-    clf = LinearSVC(C=1).fit(X, y)
+    clf = LinearSVC(dual="auto", C=1).fit(X, y)
     calib_clf = CalibratedClassifierCV(clf, cv="prefit")
 
     msg = "X has 3 features, but LinearSVC is expecting 5 features as input."
diff --git a/sklearn/tests/test_docstring_parameters.py b/sklearn/tests/test_docstring_parameters.py
index 707f60ec7fd90..9ae2196cfc20c 100644
--- a/sklearn/tests/test_docstring_parameters.py
+++ b/sklearn/tests/test_docstring_parameters.py
@@ -256,6 +256,10 @@ def test_fit_docstring_attributes(name, Estimator):
     if Estimator.__name__ in ("KMeans", "MiniBatchKMeans"):
         est.set_params(n_init="auto")
 
+    # TODO(1.4): TO BE REMOVED for 1.5 (avoid FutureWarning)
+    if Estimator.__name__ in ("LinearSVC", "LinearSVR"):
+        est.set_params(dual="auto")
+
     # TODO(1.4): TO BE REMOVED for 1.4 (avoid FutureWarning)
     if Estimator.__name__ in (
         "MultinomialNB",
diff --git a/sklearn/tests/test_multiclass.py b/sklearn/tests/test_multiclass.py
index 342d76d357263..472d1adadc050 100644
--- a/sklearn/tests/test_multiclass.py
+++ b/sklearn/tests/test_multiclass.py
@@ -53,7 +53,7 @@
 
 
 def test_ovr_exceptions():
-    ovr = OneVsRestClassifier(LinearSVC(random_state=0))
+    ovr = OneVsRestClassifier(LinearSVC(dual="auto", random_state=0))
 
     # test predicting without fitting
     with pytest.raises(NotFittedError):
@@ -82,11 +82,11 @@ def test_check_classification_targets():
 
 def test_ovr_fit_predict():
     # A classifier which implements decision_function.
-    ovr = OneVsRestClassifier(LinearSVC(random_state=0))
+    ovr = OneVsRestClassifier(LinearSVC(dual="auto", random_state=0))
     pred = ovr.fit(iris.data, iris.target).predict(iris.data)
     assert len(ovr.estimators_) == n_classes
 
-    clf = LinearSVC(random_state=0)
+    clf = LinearSVC(dual="auto", random_state=0)
     pred2 = clf.fit(iris.data, iris.target).predict(iris.data)
     assert np.mean(iris.target == pred) == np.mean(iris.target == pred2)
 
@@ -257,7 +257,7 @@ def test_ovr_multiclass():
 
     for base_clf in (
         MultinomialNB(),
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         LinearRegression(),
         Ridge(),
         ElasticNet(),
@@ -302,7 +302,7 @@ def conduct_test(base_clf, test_predict_proba=False):
         assert y_pred == 1
 
     for base_clf in (
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         LinearRegression(),
         Ridge(),
         ElasticNet(),
@@ -320,7 +320,7 @@ def test_ovr_multilabel():
 
     for base_clf in (
         MultinomialNB(),
-        LinearSVC(random_state=0),
+        LinearSVC(dual="auto", random_state=0),
         LinearRegression(),
         Ridge(),
         ElasticNet(),
@@ -458,7 +458,7 @@ def test_ovr_single_label_decision_function():
 
 
 def test_ovr_gridsearch():
-    ovr = OneVsRestClassifier(LinearSVC(random_state=0))
+    ovr = OneVsRestClassifier(LinearSVC(dual="auto", random_state=0))
     Cs = [0.1, 0.5, 0.8]
     cv = GridSearchCV(ovr, {"estimator__C": Cs})
     cv.fit(iris.data, iris.target)
@@ -479,7 +479,7 @@ def test_ovr_pipeline():
 
 
 def test_ovo_exceptions():
-    ovo = OneVsOneClassifier(LinearSVC(random_state=0))
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto", random_state=0))
     with pytest.raises(NotFittedError):
         ovo.predict([])
 
@@ -487,7 +487,7 @@ def test_ovo_exceptions():
 def test_ovo_fit_on_list():
     # Test that OneVsOne fitting works with a list of targets and yields the
     # same output as predict from an array
-    ovo = OneVsOneClassifier(LinearSVC(random_state=0))
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto", random_state=0))
     prediction_from_array = ovo.fit(iris.data, iris.target).predict(iris.data)
     iris_data_list = [list(a) for a in iris.data]
     prediction_from_list = ovo.fit(iris_data_list, list(iris.target)).predict(
@@ -498,7 +498,7 @@ def test_ovo_fit_on_list():
 
 def test_ovo_fit_predict():
     # A classifier which implements decision_function.
-    ovo = OneVsOneClassifier(LinearSVC(random_state=0))
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto", random_state=0))
     ovo.fit(iris.data, iris.target).predict(iris.data)
     assert len(ovo.estimators_) == n_classes * (n_classes - 1) / 2
 
@@ -564,7 +564,7 @@ def test_ovo_partial_fit_predict():
 def test_ovo_decision_function():
     n_samples = iris.data.shape[0]
 
-    ovo_clf = OneVsOneClassifier(LinearSVC(random_state=0))
+    ovo_clf = OneVsOneClassifier(LinearSVC(dual="auto", random_state=0))
     # first binary
     ovo_clf.fit(iris.data, iris.target == 0)
     decisions = ovo_clf.decision_function(iris.data)
@@ -609,7 +609,7 @@ def test_ovo_decision_function():
 
 
 def test_ovo_gridsearch():
-    ovo = OneVsOneClassifier(LinearSVC(random_state=0))
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto", random_state=0))
     Cs = [0.1, 0.5, 0.8]
     cv = GridSearchCV(ovo, {"estimator__C": Cs})
     cv.fit(iris.data, iris.target)
@@ -659,7 +659,7 @@ def test_ovo_string_y():
     X = np.eye(4)
     y = np.array(["a", "b", "c", "d"])
 
-    ovo = OneVsOneClassifier(LinearSVC())
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto"))
     ovo.fit(X, y)
     assert_array_equal(y, ovo.predict(X))
 
@@ -669,7 +669,7 @@ def test_ovo_one_class():
     X = np.eye(4)
     y = np.array(["a"] * 4)
 
-    ovo = OneVsOneClassifier(LinearSVC())
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto"))
     msg = "when only one class"
     with pytest.raises(ValueError, match=msg):
         ovo.fit(X, y)
@@ -680,21 +680,23 @@ def test_ovo_float_y():
     X = iris.data
     y = iris.data[:, 0]
 
-    ovo = OneVsOneClassifier(LinearSVC())
+    ovo = OneVsOneClassifier(LinearSVC(dual="auto"))
     msg = "Unknown label type"
     with pytest.raises(ValueError, match=msg):
         ovo.fit(X, y)
 
 
 def test_ecoc_exceptions():
-    ecoc = OutputCodeClassifier(LinearSVC(random_state=0))
+    ecoc = OutputCodeClassifier(LinearSVC(dual="auto", random_state=0))
     with pytest.raises(NotFittedError):
         ecoc.predict([])
 
 
 def test_ecoc_fit_predict():
     # A classifier which implements decision_function.
-    ecoc = OutputCodeClassifier(LinearSVC(random_state=0), code_size=2, random_state=0)
+    ecoc = OutputCodeClassifier(
+        LinearSVC(dual="auto", random_state=0), code_size=2, random_state=0
+    )
     ecoc.fit(iris.data, iris.target).predict(iris.data)
     assert len(ecoc.estimators_) == n_classes * 2
 
@@ -705,7 +707,7 @@ def test_ecoc_fit_predict():
 
 
 def test_ecoc_gridsearch():
-    ecoc = OutputCodeClassifier(LinearSVC(random_state=0), random_state=0)
+    ecoc = OutputCodeClassifier(LinearSVC(dual="auto", random_state=0), random_state=0)
     Cs = [0.1, 0.5, 0.8]
     cv = GridSearchCV(ecoc, {"estimator__C": Cs})
     cv.fit(iris.data, iris.target)
@@ -718,7 +720,7 @@ def test_ecoc_float_y():
     X = iris.data
     y = iris.data[:, 0]
 
-    ovo = OutputCodeClassifier(LinearSVC())
+    ovo = OutputCodeClassifier(LinearSVC(dual="auto"))
     msg = "Unknown label type"
     with pytest.raises(ValueError, match=msg):
         ovo.fit(X, y)
@@ -745,7 +747,7 @@ def test_ecoc_delegate_sparse_base_estimator():
         ecoc.predict(X_sp)
 
     # smoke test to check when sparse input should be supported
-    ecoc = OutputCodeClassifier(LinearSVC(random_state=0))
+    ecoc = OutputCodeClassifier(LinearSVC(dual="auto", random_state=0))
     ecoc.fit(X_sp, y).predict(X_sp)
     assert len(ecoc.estimators_) == 4
 
diff --git a/sklearn/tests/test_multioutput.py b/sklearn/tests/test_multioutput.py
index 371c516f197c5..c319ec08bb7b9 100644
--- a/sklearn/tests/test_multioutput.py
+++ b/sklearn/tests/test_multioutput.py
@@ -304,7 +304,7 @@ def test_multi_output_classification():
 
 def test_multiclass_multioutput_estimator():
     # test to check meta of meta estimators
-    svc = LinearSVC(random_state=0)
+    svc = LinearSVC(dual="auto", random_state=0)
     multi_class_svc = OneVsRestClassifier(svc)
     multi_target_svc = MultiOutputClassifier(multi_class_svc)
 
@@ -409,7 +409,7 @@ def test_multi_output_classification_partial_fit_sample_weights():
 def test_multi_output_exceptions():
     # NotFittedError when fit is not done but score, predict and
     # and predict_proba are called
-    moc = MultiOutputClassifier(LinearSVC(random_state=0))
+    moc = MultiOutputClassifier(LinearSVC(dual="auto", random_state=0))
     with pytest.raises(NotFittedError):
         moc.score(X, y)
 
@@ -445,7 +445,7 @@ def test_multi_output_delegate_predict_proba():
     assert hasattr(moc, "predict_proba")
 
     # A base estimator without `predict_proba` should raise an AttributeError
-    moc = MultiOutputClassifier(LinearSVC())
+    moc = MultiOutputClassifier(LinearSVC(dual="auto"))
     assert not hasattr(moc, "predict_proba")
     msg = "'LinearSVC' object has no attribute 'predict_proba'"
     with pytest.raises(AttributeError, match=msg):
@@ -471,7 +471,7 @@ def generate_multilabel_dataset_with_correlations():
 def test_classifier_chain_fit_and_predict_with_linear_svc():
     # Fit classifier chain and verify predict performance using LinearSVC
     X, Y = generate_multilabel_dataset_with_correlations()
-    classifier_chain = ClassifierChain(LinearSVC())
+    classifier_chain = ClassifierChain(LinearSVC(dual="auto"))
     classifier_chain.fit(X, Y)
 
     Y_pred = classifier_chain.predict(X)


Code review comments:
diff_hunk: @@ -1,14 +1,35 @@
 from numbers import Integral, Real
+import warnings
 
 import numpy as np
 
-from ._base import _fit_liblinear, BaseSVC, BaseLibSVM
+from ._base import _fit_liblinear, _get_liblinear_solver_type, BaseSVC, BaseLibSVM
 from ..base import BaseEstimator, RegressorMixin, OutlierMixin
 from ..linear_model._base import LinearClassifierMixin, SparseCoefMixin, LinearModel
 from ..utils import deprecated
 from ..utils.validation import _num_samples
 from ..utils.multiclass import check_classification_targets
-from ..utils._param_validation import Interval, StrOptions
+from ..utils._param_validation import Interval, StrOptions, Hidden
+
+
+def _choose_dual_automatically(loss, penalty, multi_class, X):
+    """Choose dual parameter value wrt to loss, penalty and shape of the data"""
+    if X.shape[0] < X.shape[1]:
+        try:
+            # check if the combination of parameters supports
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, True)
+            return True
+        except ValueError:
+            return False
+    else:
+        try:
+            # check if the combination of parameters doesn't support
+            # dual formulation
path: sklearn/svm/_classes.py
text: @user1:
```suggestion
```

diff_hunk: @@ -1,14 +1,35 @@
 from numbers import Integral, Real
+import warnings
 
 import numpy as np
 
-from ._base import _fit_liblinear, BaseSVC, BaseLibSVM
+from ._base import _fit_liblinear, _get_liblinear_solver_type, BaseSVC, BaseLibSVM
 from ..base import BaseEstimator, RegressorMixin, OutlierMixin
 from ..linear_model._base import LinearClassifierMixin, SparseCoefMixin, LinearModel
 from ..utils import deprecated
 from ..utils.validation import _num_samples
 from ..utils.multiclass import check_classification_targets
-from ..utils._param_validation import Interval, StrOptions
+from ..utils._param_validation import Interval, StrOptions, Hidden
+
+
+def _choose_dual_automatically(loss, penalty, multi_class, X):
+    """Choose dual parameter value wrt to loss, penalty and shape of the data"""
+    if X.shape[0] < X.shape[1]:
+        try:
+            # check if the combination of parameters supports
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, True)
+            return True
+        except ValueError:
path: sklearn/svm/_classes.py
text: @user1:
```suggestion
        except ValueError:  # dual not supported for the combination
```

diff_hunk: @@ -1,14 +1,35 @@
 from numbers import Integral, Real
+import warnings
 
 import numpy as np
 
-from ._base import _fit_liblinear, BaseSVC, BaseLibSVM
+from ._base import _fit_liblinear, _get_liblinear_solver_type, BaseSVC, BaseLibSVM
 from ..base import BaseEstimator, RegressorMixin, OutlierMixin
 from ..linear_model._base import LinearClassifierMixin, SparseCoefMixin, LinearModel
 from ..utils import deprecated
 from ..utils.validation import _num_samples
 from ..utils.multiclass import check_classification_targets
-from ..utils._param_validation import Interval, StrOptions
+from ..utils._param_validation import Interval, StrOptions, Hidden
+
+
+def _choose_dual_automatically(loss, penalty, multi_class, X):
+    """Choose dual parameter value wrt to loss, penalty and shape of the data"""
+    if X.shape[0] < X.shape[1]:
+        try:
+            # check if the combination of parameters supports
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, True)
+            return True
+        except ValueError:
+            return False
+    else:
+        try:
+            # check if the combination of parameters doesn't support
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, False)
+            return False
+        except ValueError:
path: sklearn/svm/_classes.py
text: @user1:
```suggestion
        except ValueError:  # primal not supported by the combination
```

diff_hunk: @@ -1,14 +1,35 @@
 from numbers import Integral, Real
+import warnings
 
 import numpy as np
 
-from ._base import _fit_liblinear, BaseSVC, BaseLibSVM
+from ._base import _fit_liblinear, _get_liblinear_solver_type, BaseSVC, BaseLibSVM
 from ..base import BaseEstimator, RegressorMixin, OutlierMixin
 from ..linear_model._base import LinearClassifierMixin, SparseCoefMixin, LinearModel
 from ..utils import deprecated
 from ..utils.validation import _num_samples
 from ..utils.multiclass import check_classification_targets
-from ..utils._param_validation import Interval, StrOptions
+from ..utils._param_validation import Interval, StrOptions, Hidden
+
+
+def _choose_dual_automatically(loss, penalty, multi_class, X):
+    """Choose dual parameter value wrt to loss, penalty and shape of the data"""
+    if X.shape[0] < X.shape[1]:
+        try:
+            # check if the combination of parameters supports
+            # dual formulation
+            _get_liblinear_solver_type(multi_class, penalty, loss, True)
+            return True
+        except ValueError:
+            return False
path: sklearn/svm/_classes.py
text: @user1:
Apparently, we don't have a test for this case. It could be nice to hand-craft a case that this one happen. Otherwise, it means that we don't need the `try ... except` here.