Retry failed Lanczos eigensoves with jacobi

cpp/include/raft/linalg/detail/eig.cuh

@@ -25,6 +25,12 @@ namespace raft {
 namespace linalg {
 namespace detail {
 
+/** A completed symmetric eigensolve that did not converge. */
+class eigensolver_convergence_error : public raft::exception {
+ public:
+  using raft::exception::exception;
+};
+
 template <typename math_t>
 void eigDC_legacy(raft::resources const& handle,
                   const math_t* in,
@@ -153,15 +159,18 @@ void eigDC(raft::resources const& handle,
   RAFT_CUDA_TRY(cudaGetLastError());
   RAFT_CUSOLVER_TRY(cusolverDnDestroyParams(dn_params));
   int dev_info = d_dev_info.value(stream_new);
-  ASSERT(dev_info == 0,
-         "eig.cuh: eigensolver couldn't converge to a solution. "
-         "This usually occurs when some of the features do not vary enough.");
 
   if (cusolver_version < 110603) {
     // Synchronize the created stream with the original stream before return
     RAFT_CUDA_TRY(cudaEventRecord(sync_event, stream_new));
     RAFT_CUDA_TRY(cudaStreamWaitEvent(stream, sync_event));
   }
+  if (dev_info > 0) {
+    throw eigensolver_convergence_error(
+      "eig.cuh: eigensolver couldn't converge to a solution. "
+      "This usually occurs when some of the features do not vary enough.");
+  }
+  RAFT_EXPECTS(dev_info == 0, "eig.cuh: invalid eigensolver argument (info=%d)", dev_info);
 }
 
 enum EigVecMemUsage { OVERWRITE_INPUT, COPY_INPUT };
@@ -275,7 +284,8 @@ void eigJacobi(raft::resources const& handle,
                math_t* eig_vals,
                cudaStream_t stream,
                math_t tol = 1.e-7,
-               int sweeps = 15)
+               int sweeps = 15,
+               bool check_convergence = false)
 {
   cusolverDnHandle_t cusolverH = resource::get_cusolver_dn_handle(handle);
 
@@ -326,6 +336,14 @@ void eigJacobi(raft::resources const& handle,
 
   RAFT_CUDA_TRY(cudaGetLastError());
   RAFT_CUSOLVER_TRY(cusolverDnDestroySyevjInfo(syevj_params));
+  // Only callers opting in read and check the convergence flag.
+  if (check_convergence) {
+    int info = dev_info.value(stream);
+    if (info > 0) {
+      throw eigensolver_convergence_error("eig.cuh: Jacobi eigensolver couldn't converge to a solution.");
+    }
+    RAFT_EXPECTS(info == 0, "eig.cuh: invalid Jacobi eigensolver argument (info=%d)", info);
+  }
 }
 
 }  // namespace detail

cpp/include/raft/sparse/solver/detail/lanczos.cuh

@@ -63,6 +63,7 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
+#include <limits>
 #include <optional>
 #include <type_traits>
 #include <utility>
@@ -127,6 +128,28 @@ RAFT_KERNEL kernel_clamp_down_vector(T* vec, T threshold, int size)
   if (idx < size) { vec[idx] = (fabs(vec[idx]) < threshold) ? 0 : vec[idx]; }
 }
 
+/** Retry float nonconvergence with Jacobi for the small projected problem.
+ * Input and output must not alias: eigDC overwrites output even on failure.
+ * CUDA/API errors and a failed Jacobi solve propagate to the caller.
+ */
+template <typename T>
+void lanczos_eig_dc_with_fallback(
+  raft::resources const& handle, const T* input, uint32_t n, T* vectors, T* values)
+{
+  cudaStream_t stream = resource::get_cuda_stream(handle).get();
+  try {
+    raft::linalg::eigDC(handle, input, n, n, vectors, values, stream);
+  } catch (raft::linalg::detail::eigensolver_convergence_error const&) {
+    if constexpr (std::is_same_v<T, float>) {
+      // Retry the untouched projected matrix, not the failed eigDC output.
+      raft::linalg::detail::eigJacobi(
+        handle, input, n, n, vectors, values, stream, std::numeric_limits<T>::epsilon(), 100, true);
+    } else {
+      throw;
+    }
+  }
+}
+
 template <typename IndexTypeT, typename ValueTypeT>
 void lanczos_solve_ritz(
   raft::resources const& handle,
@@ -179,11 +202,11 @@ void lanczos_solve_ritz(
                         ncv);
   }
 
-  auto triangular_matrix_view =
-    raft::make_device_matrix_view<const ValueTypeT, uint32_t, raft::col_major>(
-      triangular_matrix.data_handle(), ncv, ncv);
-
-  raft::linalg::eig_dc(handle, triangular_matrix_view, eigenvectors, eigenvalues);
+  lanczos_eig_dc_with_fallback(handle,
+                               triangular_matrix.data_handle(),
+                               static_cast<uint32_t>(ncv),
+                               eigenvectors.data_handle(),
+                               eigenvalues.data_handle());
 
   IndexTypeT nEigVecs = k;
 

cpp/tests/CMakeLists.txt

@@ -252,6 +252,7 @@ if(BUILD_TESTS)
     linalg/eigen_solvers.cu
     lap/lap.cu
     sparse/mst.cu
+    sparse/solver/lanczos_fallback.cu
     sparse/solver/lanczos.cu
     sparse/solver/lanczos_svds.cu
     sparse/solver/randomized_svds.cu

cpp/tests/sparse/solver/lanczos_fallback.cu

@@ -0,0 +1,57 @@
+/*
+ * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+ * SPDX-License-Identifier: Apache-2.0
+ */
+#include <raft/core/resources.hpp>
+#include <raft/sparse/solver/detail/lanczos.cuh>
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cmath>
+
+namespace raft::sparse::solver {
+TEST(LanczosEigensolverFallback, Float32Jacobi)
+{
+  // Reduced cuSOLVER reproducer. Symmetry also makes this column-major input.
+  // clang-format off
+  constexpr std::array<float, 25> input = {
+     4.52996005e-8f,   6.31930868e-13f,  0.0f,             0.0f,             0.0f,
+     6.31930868e-13f, -2.11824158e-9f,  -2.23728129e-8f,  0.0f,             0.0f,
+     0.0f,           -2.23728129e-8f,  -2.36300252e-7f,  4.03643658e-14f,  0.0f,
+     0.0f,            0.0f,            4.03643658e-14f, -2.38418579e-7f,   0.0f,
+     0.0f,            0.0f,            0.0f,             0.0f,             1.0f,
+  };
+  // clang-format on
+  constexpr std::array<float, 5> expected = {
+    -2.38418596e-7f, -2.38418484e-7f, 0.0f, 4.52996005e-8f, 1.0f};
+
+  raft::resources handle;
+  auto stream = raft::resource::get_cuda_stream(handle).get();
+  rmm::device_uvector<float> d_input(25, stream), d_vectors(25, stream), d_values(5, stream);
+  raft::update_device(d_input.data(), input.data(), input.size(), stream);
+
+  try {
+    raft::linalg::eigDC(handle, d_input.data(), 5, 5, d_vectors.data(), d_values.data(), stream);
+    GTEST_SKIP() << "This cuSOLVER version solves the reproducer without a fallback.";
+  } catch (raft::linalg::detail::eigensolver_convergence_error const&) {
+  }
+  ASSERT_NO_THROW(detail::lanczos_eig_dc_with_fallback(
+    handle, d_input.data(), 5, d_vectors.data(), d_values.data()));
+
+  std::array<float, 25> vectors, unchanged;
+  std::array<float, 5> values;
+  raft::update_host(vectors.data(), d_vectors.data(), vectors.size(), stream);
+  raft::update_host(values.data(), d_values.data(), values.size(), stream);
+  raft::update_host(unchanged.data(), d_input.data(), unchanged.size(), stream);
+  raft::resource::sync_stream(handle);
+
+  EXPECT_EQ(input, unchanged);
+  for (int i = 0; i < 5; ++i) {
+    EXPECT_NEAR(values[i], expected[i], 1e-7f);
+  }
+  for (float value : vectors) {
+    EXPECT_TRUE(std::isfinite(value));
+  }
+}
+}  // namespace raft::sparse::solver