I am worried this is a breaking change for users who already have code running with 1d Y shaped as (n_samples, 1). Would it work for the imputation meta estimator if we made the estimator record the shape of Y at fit time to reuse it at predict time?

    X = np.random.randn(10, 3)
    Y = np.random.randn(10, 1)

    pls = Estimator(n_components=Y.shape[1])
    pls.fit(X, Y)

    X_test = np.random.randn(10, 3)
    Y_pred = pls.predict(X_test)

    assert Y_pred.shape == (10, 1)

    pls = Estimator(n_components=Y.shape[1])
    pls.fit(X,  np.ravel(Y))

    X_test = np.random.randn(10, 3)
    Y_pred = pls.predict(X_test)

    assert Y_pred.shape == (10,)

I guess that'd be fine.

We might not need to explicitly record y.shape as it should be embedded already in self.y_weights_ and self.y_loadings_

But for what it's worth, we're not consistent at all w.r.t. the output shape when passing y with y.shape == (n_samples, 1). Here are a few estimators:

Estimator                     , output shape   , 'supports multioutput'
HistGradientBoostingClassifier, (100,)         , False
HistGradientBoostingRegressor , (100,)         , False
RandomForestClassifier        , (100,)         , True
RandomForestRegressor         , (100,)         , True
LogisticRegression            , (100,)         , False
LinearRegression              , (100, 1)       , True
SGDClassifier                 , (100,)         , False
SGDRegressor                  , (100,)         , False
SVC                           , (100,)         , False
PLSRegression                 , (100, 1)       , True

Those that don't support multioutput will warn with "DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel()."
It's expected that these ones will output (n_samples,), but for example the RandomForest estimators support multioutput and still ravel the predictions, unlike PLSRegression and LinearRegression.

For ref, the our regressors checks ensure that y_pred.shape == y.shape and our classifiers checks ensure that y_pred.shape == (n_samples,), but because of the data that these checks use, neither of them can detect such discrepancies.