BUG: fix stability of partial_fit pinv()

[tylerjereddy]

• Related to TST, CI: test_partial_fit is failing in Linux CI (with NumPy 2.5.0) #117, but does not really fix it. This just avoids calculating the Moore-Penrose pseudoinverse of a Hermitian matrix with NumPy because of the problematic vendored OpenBLAS, instead using the equivalent SciPy operation since stable releases of SciPy currently use an older OpenBLAS that doesn't have this stability issue.

• Based on the analysis in the above issue, it seems that we are safe to assume that when partial_fit receives an array via the pinv codepath, the array is always symmetric/Hermitian. Enforcing this assumption in the code seems to allow a more numerically stable outcome that allows the failing test for the ensemble estimator in the above issue to pass locally on this branch on x86_64 Linux with latest NumPy (2.5.1) wheel/updated OpenBLAS that ships with it.

• It would be better if we also had a regression test to justify the need for this change with the non-ensemble estimators, but we haven't succeeded in finding one. Nonetheless, this does appear to offer slightly more reliable numerical outcomes for pinv in general, given that we can guarantee symmetric matrices via the Gram matrix control flow (for partial fits).

• Emma had done some testing on the alternative suggestion in that issue, using lstsq instead, but did not have success in getting all tests passing.

• I don't believe we can make the same symmetry assumption for the full fit based on some rudimentary testing while on travel.

---

AI usage disclosure: an LLM was consulted on the likely safety of assuming hermitian=True as described in the analysis in the cognate issue.

[tylerjereddy]

Review of this PR was requested from @aganguly-lab a week or two ago. If you find a reproducer for the non-ensemble estimator that's great, but if you don't then please summarize your review findings/opinion on the changes and then I will try to circle back to see if I can help with that reproducer again.

[tylerjereddy]

Also, note that it looks like other research teams are starting to be able to reproduce similar findings/problems upstream in the last few days: OpenMathLib/OpenBLAS#5866 (comment)

[aganguly-lab]

I think this can be merged. I wrote a regression test that fails the partial_fit method of the non-ensemble classifier under the exact same conditions of Issue 117. That is, it fails on linux OS with numpy versions 2.5.0-2.5.2, and it passes when the changes of this merge request are applied. This test should pass on MacOS and for numpy versions 2.4.6 or below even without the proposed merge.

The test:

# Both fixtures were written by AI
@pytest.fixture
def fixtures_dir():
    """Returns the fixtures directory path"""
    return Path(__file__).parent / "fixtures"


@pytest.fixture
def Hermitian_data(fixtures_dir):
    X = np.load(fixtures_dir / "X_Hermitian.npy")
    y = np.load(fixtures_dir / "y_Hermitian.npy")
    return (X, y)


@pytest.mark.parametrize("hidden_layer_sizes", [(5, 5)])
@pytest.mark.parametrize("n_classes", [2])
@pytest.mark.parametrize("activation", ["identity"])
@pytest.mark.parametrize("weight_init", ["zeros"])
@pytest.mark.parametrize("alpha", [None])
@pytest.mark.parametrize("direct_links", [True])
@pytest.mark.parametrize("attr", ["coeff_"])
def test_GFDLClassifier_Pinv_Hermitian_True(
    Hermitian_data,
    hidden_layer_sizes,
    direct_links,
    n_classes,
    activation,
    weight_init,
    alpha,
    attr
):
    # Regression test for MR #120: recreates issue #117 in GFDLClassifier.
    # Generates data for which partial fit of the non-ensemble classifier fails.
    # Data is set to the design matrix of the Ensemble Classifier whose failed
    # test caused Issue #117. Weights are set to 0 simplify and stabilize generation
    # of the design matrix.
    # This should be a very rare issue.
    # The classifier fails without stabilizing np.linalg.pinv
    # using the hermitian = True option.
    # The test passes regardless on Mac or numpy version 2.4.6 or below.

    X, y = Hermitian_data

    kwargs = {
        "hidden_layer_sizes": hidden_layer_sizes,
        "activation": activation,
        "weight_scheme": weight_init,
        "seed": 0,
        "reg_alpha": alpha,
    }
    if direct_links is not None:
        kwargs["direct_links"] = direct_links

    ff_model = GFDLClassifier(**kwargs)
    pf_model = clone(ff_model)

    ff_model.fit(X, y)

    classes = np.unique(y)
    batch = 25
    for start in range(0, len(X), batch):
        end = min(start + batch, len(X))
        Xb = X[start:end]
        yb = y[start:end]
        if start == 0:
            pf_model.partial_fit(Xb, yb, classes=classes)
        else:
            pf_model.partial_fit(Xb, yb)

    assert_allclose(
        getattr(pf_model, attr), getattr(ff_model, attr), rtol=1e-5, atol=1e-3
        )

Instructions to run the test:

1. Add the following import to test_model.py: from pathlib import Path
2. Copy the code above to the end of the file.
3. Add a directory to the tests folder named fixtures.
4. Download the attached text files, and convert them into .npy files with the same name.
5. Save the .npy files in fixtures

Method: The test uses 0 weights and the identity activation function to reduce the RVFL Classifier to a linear classifier. The data (X,y) was generated such that design matrix of the resulting non-ensemble methods matched the design matrix from the failing test in Issue 117. I was unable to generate random synthetic data that replicated this error.

AI Use: AI was used for developing path logic (the two pytest fixtures and the file structure). The rest of the test was created without the use of AI.
y_Hermitian.txt
X_Hermitian.txt

[tylerjereddy]

@aganguly-lab thanks. I've pushed up a simplified version of the regression test here now, that you should check carefully.

I think this can be merged.

Are you sure? The CI is failing on x86_64 Linux, and I can reproduce locally in that scenario on this branch as well.

If I run python -m pytest "src/gfdl/tests/test_model.py::test_partial_fit_performance[GFDLClassifier-None-None-0.7631578947368421]" with NumPy 2.5.2 on this branch it seems to fail consistently. Did you check the full testsuite on x86_64 Linux against latest NumPy here?

On the other hand, if I revert the Hermitian flag changes that failing test starts passing again.

Does it make sense that our estimator performance drops on the Breast cancer standard dataset when using pinv at the default singular value cutoff when we assume symmetry? If so, we can just change that test as well, but if not, it may require more thought.

The traceback for that new failing test is visible in the CI here, and also below the fold.

Details

______________________________________ test_partial_fit_performance[GFDLClassifier-None-None-0.7631578947368421] _______________________________________

Classifier = <class 'gfdl.model.GFDLClassifier'>, ridge_alpha = None, rcond = None, expected = 0.7631578947368421

    @pytest.mark.parametrize(
        "Classifier, ridge_alpha, rcond, expected",
        [
            # NOTE: for Moore-Penrose, a large singular value
            # cutoff (rcond) or alpha regularization is required
            # to achieve reasonable accuracy
            # Without rtol or reg accuracy ~= 0.7632
            # With reg accuracy ~= 0.9649
            # With rtol accuracy ~= 0.9474
            (GFDLClassifier, 5, None, 0.9649122807017544),
            (GFDLClassifier, None, None, 0.7631578947368421),
            (GFDLClassifier, None, 1e-3, 0.9473684210526315),
    
            # NOTE: for Moore-Penrose, a large singular value
            # cutoff (rcond) or alpha regularization is required
            # to achieve reasonable accuracy
            # Without rtol or reg accuracy ~= 0.7456
            # With reg accuracy ~= 0.9737
            # With rtol accuracy ~= 0.9561
            (EnsembleGFDLClassifier, 5, None, 0.9736842105263158),
            (EnsembleGFDLClassifier, None, None, 0.7456140350877193),
            (EnsembleGFDLClassifier, None, 1e-3, 0.956140350877193),
            # NOTE: this behavior may be exacerbated by using shallower,
            # wider networks
            # With this dataset in particular we're achieving better accuracies
            # with smaller architectures, but for the sake of this test
            # we're using a shallow and wide architecture
        ],
    )
    def test_partial_fit_performance(Classifier, ridge_alpha, rcond, expected):
        X, y = load_breast_cancer(return_X_y=True)
    
        X_train, X_test, y_train, y_test = train_test_split(
            X, y, test_size=0.2, random_state=42, shuffle=True
        )
    
        scaler = StandardScaler()
        X_train = scaler.fit_transform(X_train)
        X_test = scaler.transform(X_test)
        classes = np.unique(y)
    
        model = Classifier(
            hidden_layer_sizes=[500]
            * 3,  # partial_fit is slow so smaller network for speed
            activation="tanh",
            weight_scheme="uniform",
            seed=0,
            reg_alpha=ridge_alpha,
            rtol=rcond,
        )
    
        batch = 20
        for start in range(0, X_train.shape[0], batch):
            end = min(start + batch, X_train.shape[0])
            if start == 0:
                model.partial_fit(X_train[start:end], y_train[start:end], classes=classes)
            else:
                model.partial_fit(X_train[start:end], y_train[start:end])
    
        actual = model.score(X_test, y_test)
        # RandomForestClassifier() with default params scores 0.973 here
        # RVFL with above params scores comparatively:
>       assert_allclose(actual, expected)
E       AssertionError: 
E       Not equal to tolerance rtol=1e-07, atol=0
E       
E       Mismatched elements: 1 / 1 (100%)
E       Max absolute difference among violations: 0.03508772
E       Max relative difference among violations: 0.04597701
E        ACTUAL: array(0.72807)
E        DESIRED: array(0.763158)

src/gfdl/tests/test_model.py:672: AssertionError

Note that reverting to NumPy 2.5.1 also allows the test to pass, which is why the CI was passing here originally. I haven't checked if there was another OpenBLAS bump in 2.5.2 (usually not in bug fix releases, but can happen sometimes...).

[aganguly-lab]

Sorry about that. I think I tried the full suite of tests for 2.5.0, but I may have only run my new test and the previously failing test on 2.5.2. I'm running the test suite now. If pinv is unstable even with the hermitian=True argument, should we try @ilayn 's solution from Issue #117 (using lstsq instead of pinv)? Alternatively, could we try setting a larger default rtol value?

I can test either solution if you think they are viable.

[tylerjereddy]

It might be good to determine what changed in NumPy 2.5.2 to trigger the change in model performance for us. Maybe check their release notes for anything relevant, bisect on NumPy if needed to find the first commit where we start failing this new test to understand a bit better.

I'm hesitant to change approaches until there's a clearer understanding of the problem.

[tylerjereddy]

This is to allow installation of the new .npy testing files added via i.e., pip install ...

[tylerjereddy]

As far as I can tell, the directory that was previously here no longer exists, so can be replaced.

[tylerjereddy]

This is a simplified/condensed/concise version of the test suggested at https://github.com/lanl/GFDL/pull/120/files#r3797463241.

[tylerjereddy]

@aganguly-lab I've pushed in a commit to switch the calculation of the Moore-Penrose pseudoinverse of a Hermitian matrix from NumPy to the SciPy equivalent, since SciPy binaries currently ship with an older version of OpenBLAS without the stability problem. In my hands, that allows the full testsuite to pass on x86_64 Linux with NumPy 2.4.6, 2.5.1, and 2.5.2.

I'd say this no longer fixes the target issue, and instead just temporarily circumvents it until SciPy releases start using newer OpenBLAS (which is inevitable), so I've adjusted the original comment above to no longer auto-close the matching issue--this is now just a temporary workaround.

[ilayn]

pinv and pinvh are unstable in the sense that the returning singular vectors might swap signs or the their number of vectors depending on how you decide to truncate. That's why using pinv as standalone term in the computations are slightly problematic. You might be having bunch of SVs at the truncation value. Or that particular architecture has a sensitivity in BLAS layer

Typically the common path is using lstsq to have a norm guaranteed projection as the solution without the pitfalls of SVD based choices.

[tylerjereddy]

We already have a codepath that may use least squares via sklearn Ridge, which is indeed generally more stable, both in our own experience and in the generally published literature around RVFL estimators.

I believe @eiviani-lanl hasn't had great luck subbing in lstsq here, but we could try it again maybe. The RVFL literature usually mentions MP, but as long as both are closed form the substitution is probably fine.

[tylerjereddy]

With the diff below the fold applied to this PR branch, to switch all pinv uses to np.linalg.lstsq,

Details

--- a/src/gfdl/model.py
+++ b/src/gfdl/model.py
@@ -16,7 +16,7 @@ from sklearn.linear_model import Ridge
 from sklearn.preprocessing import OneHotEncoder
 from sklearn.utils import column_or_1d
 from sklearn.utils.metaestimators import available_if
-from sklearn.utils.multiclass import check_classification_targets, unique_labels
+from sklearn.utils.multiclass import check_classification_targets, unique_labels, type_of_target
 from sklearn.utils.validation import check_is_fitted, validate_data
 
 from gfdl.activations import resolve_activation
@@ -47,6 +47,7 @@ class GFDL(BaseEstimator):
         # Assumption : X, Y have been pre-processed.
         # X shape: (n_samples, n_features)
         # Y shape: (n_samples, n_classes-1)
+        type_of_target(Y, raise_unknown=True)
         if self.reg_alpha is not None and self.reg_alpha < 0.0:
             raise ValueError("Negative reg_alpha. Expected range : None or [0.0, inf).")
         hidden_layer_sizes = np.asarray(self.hidden_layer_sizes)
@@ -103,7 +104,7 @@ class GFDL(BaseEstimator):
         # If reg_alpha is None, use direct solve using
         # MoorePenrose Pseudo-Inverse, otherwise use ridge regularized form.
         if self.reg_alpha is None:
-            self.coeff_ = np.linalg.pinv(D, rtol=self.rtol) @ Y
+            self.coeff_ = np.linalg.lstsq(D, Y, rcond=self.rtol)[0]
         else:
             ridge = Ridge(alpha=self.reg_alpha, fit_intercept=False)
             ridge.fit(D, Y)
@@ -215,7 +216,7 @@ class GFDL(BaseEstimator):
         if self.reg_alpha is None:
             # NOTE: the equivalent NumPy operation on Hermitian matrices
             # is currently unstable because of OpenBLAS per gh-117
-            self.coeff_ = scipy.linalg.pinvh(self.A, rtol=self.rtol) @ self.B
+            self.coeff_ = np.linalg.lstsq(self.A, self.B, rcond=self.rtol)[0]
         else:
             # scipy.linalg.solve(self.A + reg_mat, self.B)
             # is equivalent to
@@ -616,7 +617,7 @@ class EnsembleGFDL(BaseEstimator):
             # If reg_alpha is None, use direct solve using
             # MoorePenrose Pseudo-Inverse, otherwise use ridge regularized form.
             if self.reg_alpha is None:
-                coeff = np.linalg.pinv(D, rtol=self.rtol) @ Y
+                coeff = np.linalg.lstsq(D, Y, rcond=self.rtol)[0]
             else:
                 ridge = Ridge(alpha=self.reg_alpha, fit_intercept=False)
                 ridge.fit(D, Y)
@@ -704,7 +705,7 @@ class EnsembleGFDL(BaseEstimator):
             if self.reg_alpha is None:
                 # NOTE: the equivalent NumPy operation on Hermitian matrices
                 # is currently unstable because of OpenBLAS per gh-117
-                coef_ = scipy.linalg.pinvh(self.As[i], rtol=self.rtol) @ self.Bs[i]
+                coef_ = np.linalg.lstsq(self.As[i], self.Bs[i], rcond=self.rtol)[0]
             else:
                 # scipy.linalg.solve(self.A + reg_mat, self.B)
                 # is equivalent to

using NumPy 2.5.2 produces the exact same two test failures we've seen previously with test_partial_fit and test_gh_117 with pinv approach, so it really is just back to square one if we do that. Just as for pinv, moving back to the NumPy 2.4.x binary wheel series allows the full testsuite to pass with lstsq as well.

So, as it stands, I'm still standing by the suggestion that, for now, our best bet is to use a solution that doesn't leverage that newer OpenBLAS as done here. If lstsq is genuinely more stable for these very pathological scenarios (and these are basically ridiculous scenarios), it isn't showing with the newer OpenBLAS and we end up with an identical stability issue with np.linalg.lstsq. If we have a scenario where lstsq is provably better (i.e., another regression test to enforce that), that seems totally fine by me, but that isn't the case here--the only viable solution I see is to avoid the new OpenBLAS shipped with some binaries for now.

[ilayn]

Then indeed there is an issue somewhere in the Reference-LAPACK like Martin mentioned.

[tylerjereddy]

I carefully checked a few versions of reference LAPACK on x86_64 Linux and was not able to find evidence of the issue/reproduction there, after a few hours of analysis in #117 (comment). That included reference LAPACK v3.12.1, d134ca0070b0 (latest master of reference lapack), and v3.9.1 (about 5 years ago).