diff --git a/include/ceres/solver.h b/include/ceres/solver.h
index 2a1dd18..27198ac 100644
--- a/include/ceres/solver.h
+++ b/include/ceres/solver.h
@@ -691,6 +691,23 @@
     // as its termination type.
     int max_linear_solver_iterations = 500;
 
+    // Maximum number of iterations performed by SCHUR_POWER_SERIES_EXPANSION.
+    // This value controls the maximum number of iterations whether it is used
+    // as a preconditioner or just to initialize the solution for
+    // ITERATIVE_SCHUR.
+    int max_num_spse_iterations = 5;
+
+    // Use SCHUR_POWER_SERIES_EXPANSION to initialize the solution for
+    // ITERATIVE_SCHUR. This option can be set true regardless of what
+    // preconditioner is being used.
+    bool use_power_series_expansion_initialization = false;
+
+    // When use_power_series_expansion_initialization is true, this parameter
+    // along with max_num_spse_iterations controls the number of
+    // SCHUR_POWER_SERIES_EXPANSION iterations performed for initialization. It
+    // is not used to control the preconditioner.
+    double spse_tolerance = 0.1;
+
     // Forcing sequence parameter. The truncated Newton solver uses
     // this number to control the relative accuracy with which the
     // Newton step is computed.
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 52e1e7b..1b07ab3 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -107,9 +107,6 @@
 
   // Use power series expansion to approximate the inversion of Schur complement
   // as a preconditioner.
-  // WARNING! Application of this preconditioner currently is not integrated
-  // into linear solvers, so failure to use it via public API is expected
-  // behaviour.
   SCHUR_POWER_SERIES_EXPANSION,
 
   // Visibility clustering based preconditioners.
diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc
index 92631de..8f7a2a0 100644
--- a/internal/ceres/implicit_schur_complement.cc
+++ b/internal/ceres/implicit_schur_complement.cc
@@ -57,6 +57,7 @@
   b_ = b;
 
   compute_ftf_inverse_ =
+      options_.use_power_series_expansion_initialization ||
       options_.preconditioner_type == JACOBI ||
       options_.preconditioner_type == SCHUR_POWER_SERIES_EXPANSION;
 
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc
index 5f56be2..8be33bc 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -43,6 +43,7 @@
 #include "ceres/implicit_schur_complement.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/linear_solver.h"
+#include "ceres/power_series_expansion_preconditioner.h"
 #include "ceres/preconditioner.h"
 #include "ceres/schur_jacobi_preconditioner.h"
 #include "ceres/triplet_sparse_matrix.h"
@@ -90,9 +91,18 @@
     return summary;
   }
 
-  // Initialize the solution to the Schur complement system to zero.
+  // Initialize the solution to the Schur complement system.
   reduced_linear_system_solution_.resize(schur_complement_->num_rows());
   reduced_linear_system_solution_.setZero();
+  if (options_.use_power_series_expansion_initialization) {
+    PowerSeriesExpansionPreconditioner pse_solver(
+        schur_complement_.get(),
+        options_.max_num_spse_iterations,
+        options_.spse_tolerance);
+    pse_solver.RightMultiplyAndAccumulate(
+        schur_complement_->rhs().data(),
+        reduced_linear_system_solution_.data());
+  }
 
   CreatePreconditioner(A);
   if (preconditioner_ != nullptr) {
@@ -168,6 +178,12 @@
       preconditioner_ = std::make_unique<SparseMatrixPreconditionerWrapper>(
           schur_complement_->block_diagonal_FtF_inverse());
       break;
+    case SCHUR_POWER_SERIES_EXPANSION:
+      // Ignoring the value of spse_tolerance to ensure preconditioner stays
+      // fixed during the iterations of cg.
+      preconditioner_ = std::make_unique<PowerSeriesExpansionPreconditioner>(
+          schur_complement_.get(), options_.max_num_spse_iterations, 0);
+      break;
     case SCHUR_JACOBI:
       preconditioner_ = std::make_unique<SchurJacobiPreconditioner>(
           *A->block_structure(), preconditioner_options);
diff --git a/internal/ceres/iterative_schur_complement_solver_test.cc b/internal/ceres/iterative_schur_complement_solver_test.cc
index 80e388a..a98b295 100644
--- a/internal/ceres/iterative_schur_complement_solver_test.cc
+++ b/internal/ceres/iterative_schur_complement_solver_test.cc
@@ -74,7 +74,9 @@
     num_eliminate_blocks_ = problem->num_eliminate_blocks;
   }
 
-  AssertionResult TestSolver(double* D) {
+  AssertionResult TestSolver(double* D,
+                             PreconditionerType preconditioner_type,
+                             bool use_power_series_expansion_initialization) {
     TripletSparseMatrix triplet_A(
         A_->num_rows(), A_->num_cols(), A_->num_nonzeros());
     A_->ToTripletSparseMatrix(&triplet_A);
@@ -95,7 +97,10 @@
     options.elimination_groups.push_back(num_eliminate_blocks_);
     options.elimination_groups.push_back(0);
     options.max_num_iterations = num_cols_;
-    options.preconditioner_type = SCHUR_JACOBI;
+    options.max_num_spse_iterations = 1;
+    options.use_power_series_expansion_initialization =
+        use_power_series_expansion_initialization;
+    options.preconditioner_type = preconditioner_type;
     IterativeSchurComplementSolver isc(options);
 
     Vector isc_sol(num_cols_);
@@ -119,16 +124,30 @@
   std::unique_ptr<double[]> D_;
 };
 
-TEST_F(IterativeSchurComplementSolverTest, NormalProblem) {
+TEST_F(IterativeSchurComplementSolverTest, NormalProblemSchurJacobi) {
   SetUpProblem(2);
-  EXPECT_TRUE(TestSolver(nullptr));
-  EXPECT_TRUE(TestSolver(D_.get()));
+  EXPECT_TRUE(TestSolver(nullptr, SCHUR_JACOBI, false));
+  EXPECT_TRUE(TestSolver(D_.get(), SCHUR_JACOBI, false));
+}
+
+TEST_F(IterativeSchurComplementSolverTest,
+       NormalProblemSchurJacobiWithPowerSeriesExpansionInitialization) {
+  SetUpProblem(2);
+  EXPECT_TRUE(TestSolver(nullptr, SCHUR_JACOBI, true));
+  EXPECT_TRUE(TestSolver(D_.get(), SCHUR_JACOBI, true));
+}
+
+TEST_F(IterativeSchurComplementSolverTest,
+       NormalProblemPowerSeriesExpansionPreconditioner) {
+  SetUpProblem(5);
+  EXPECT_TRUE(TestSolver(nullptr, SCHUR_POWER_SERIES_EXPANSION, false));
+  EXPECT_TRUE(TestSolver(D_.get(), SCHUR_POWER_SERIES_EXPANSION, false));
 }
 
 TEST_F(IterativeSchurComplementSolverTest, ProblemWithNoFBlocks) {
   SetUpProblem(3);
-  EXPECT_TRUE(TestSolver(nullptr));
-  EXPECT_TRUE(TestSolver(D_.get()));
+  EXPECT_TRUE(TestSolver(nullptr, SCHUR_JACOBI, false));
+  EXPECT_TRUE(TestSolver(D_.get(), SCHUR_JACOBI, false));
 }
 
 }  // namespace internal
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index 4916e41..69e984f 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -165,6 +165,23 @@
     int min_num_iterations = 1;
     int max_num_iterations = 1;
 
+    // Maximum number of iterations performed by SCHUR_POWER_SERIES_EXPANSION.
+    // This value controls the maximum number of iterations whether it is used
+    // as a preconditioner or just to initialize the solution for
+    // ITERATIVE_SCHUR.
+    int max_num_spse_iterations = 5;
+
+    // Use SCHUR_POWER_SERIES_EXPANSION to initialize the solution for
+    // ITERATIVE_SCHUR. This option can be set true regardless of what
+    // preconditioner is being used.
+    bool use_power_series_expansion_initialization = false;
+
+    // When use_power_series_expansion_initialization is true, this parameter
+    // along with max_num_spse_iterations controls the number of
+    // SCHUR_POWER_SERIES_EXPANSION iterations performed for initialization. It
+    // is not used to control the preconditioner.
+    double spse_tolerance = 0.1;
+
     // If possible, how many threads can the solver use.
     int num_threads = 1;
 
diff --git a/internal/ceres/power_series_expansion_preconditioner.cc b/internal/ceres/power_series_expansion_preconditioner.cc
index 1b2dbb9..7a36d92 100644
--- a/internal/ceres/power_series_expansion_preconditioner.cc
+++ b/internal/ceres/power_series_expansion_preconditioner.cc
@@ -34,13 +34,11 @@
 
 PowerSeriesExpansionPreconditioner::PowerSeriesExpansionPreconditioner(
     const ImplicitSchurComplement* isc,
-    const double spse_tolerance,
-    const int min_num_iterations,
-    const int max_num_iterations)
+    const int max_num_spse_iterations,
+    const double spse_tolerance)
     : isc_(isc),
-      spse_tolerance_(spse_tolerance),
-      min_num_iterations_(min_num_iterations),
-      max_num_iterations_(max_num_iterations) {}
+      max_num_spse_iterations_(max_num_spse_iterations),
+      spse_tolerance_(spse_tolerance) {}
 
 PowerSeriesExpansionPreconditioner::~PowerSeriesExpansionPreconditioner() =
     default;
@@ -66,8 +64,7 @@
     isc_->InversePowerSeriesOperatorRightMultiplyAccumulate(
         previous_series_term.data(), series_term.data());
     yref += series_term;
-    if (i >= min_num_iterations_ &&
-        (i >= max_num_iterations_ || series_term.norm() < norm_threshold)) {
+    if (i >= max_num_spse_iterations_ || series_term.norm() < norm_threshold) {
       break;
     }
     std::swap(previous_series_term, series_term);
diff --git a/internal/ceres/power_series_expansion_preconditioner.h b/internal/ceres/power_series_expansion_preconditioner.h
index c5cf32d..a8bb9a6 100644
--- a/internal/ceres/power_series_expansion_preconditioner.h
+++ b/internal/ceres/power_series_expansion_preconditioner.h
@@ -46,9 +46,8 @@
     : public Preconditioner {
  public:
   PowerSeriesExpansionPreconditioner(const ImplicitSchurComplement* isc,
-                                     const double spse_tolerance_,
-                                     const int min_num_iterations_,
-                                     const int max_num_iterations_);
+                                     const int max_num_spse_iterations,
+                                     const double spse_tolerance);
   PowerSeriesExpansionPreconditioner(
       const PowerSeriesExpansionPreconditioner&) = delete;
   void operator=(const PowerSeriesExpansionPreconditioner&) = delete;
@@ -60,9 +59,8 @@
 
  private:
   const ImplicitSchurComplement* isc_;
+  const int max_num_spse_iterations_;
   const double spse_tolerance_;
-  const int min_num_iterations_;
-  const int max_num_iterations_;
 };
 
 }  // namespace ceres::internal
diff --git a/internal/ceres/power_series_expansion_preconditioner_test.cc b/internal/ceres/power_series_expansion_preconditioner_test.cc
index a101831..a291405 100644
--- a/internal/ceres/power_series_expansion_preconditioner_test.cc
+++ b/internal/ceres/power_series_expansion_preconditioner_test.cc
@@ -81,10 +81,9 @@
 TEST_F(PowerSeriesExpansionPreconditionerTest,
        InverseValidPreconditionerToleranceReached) {
   const double spse_tolerance = kEpsilon;
-  const int min_iterations = 1;
-  const int max_iterations = 50;
+  const int max_num_iterations = 50;
   PowerSeriesExpansionPreconditioner preconditioner(
-      isc_.get(), spse_tolerance, min_iterations, max_iterations);
+      isc_.get(), max_num_iterations, spse_tolerance);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
@@ -104,10 +103,10 @@
 
 TEST_F(PowerSeriesExpansionPreconditionerTest,
        InverseValidPreconditionerMaxIterations) {
-  const double spse_tolerance = 1 / kEpsilon;
-  const int num_iterations = 50;
+  const double spse_tolerance = 0;
+  const int max_num_iterations = 50;
   PowerSeriesExpansionPreconditioner preconditioner_fixed_n_iterations(
-      isc_.get(), spse_tolerance, num_iterations, num_iterations);
+      isc_.get(), max_num_iterations, spse_tolerance);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
@@ -129,10 +128,9 @@
 TEST_F(PowerSeriesExpansionPreconditionerTest,
        InverseInvalidBadPreconditionerTolerance) {
   const double spse_tolerance = 1 / kEpsilon;
-  const int min_iterations = 1;
-  const int max_iterations = 50;
+  const int max_num_iterations = 50;
   PowerSeriesExpansionPreconditioner preconditioner_bad_tolerance(
-      isc_.get(), spse_tolerance, min_iterations, max_iterations);
+      isc_.get(), max_num_iterations, spse_tolerance);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
@@ -151,9 +149,9 @@
 TEST_F(PowerSeriesExpansionPreconditionerTest,
        InverseInvalidBadPreconditionerMaxIterations) {
   const double spse_tolerance = kEpsilon;
-  const int num_iterations = 1;
+  const int max_num_iterations = 1;
   PowerSeriesExpansionPreconditioner preconditioner_bad_iterations_limit(
-      isc_.get(), spse_tolerance, num_iterations, num_iterations);
+      isc_.get(), max_num_iterations, spse_tolerance);
 
   Vector x(num_f_cols_), y(num_f_cols_);
   for (int i = 0; i < num_f_cols_; i++) {
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index b07d4f0..71b4537 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -259,12 +259,25 @@
     return false;
   }
 
-  if (options.use_explicit_schur_complement &&
-      options.preconditioner_type != SCHUR_JACOBI) {
-    *error =
-        "use_explicit_schur_complement only supports "
-        "SCHUR_JACOBI as the preconditioner.";
-    return false;
+  if (options.use_explicit_schur_complement) {
+    if (options.preconditioner_type != SCHUR_JACOBI) {
+      *error =
+          "use_explicit_schur_complement only supports "
+          "SCHUR_JACOBI as the preconditioner.";
+      return false;
+    }
+    if (options.use_power_series_expansion_initialization) {
+      *error =
+          "use_explicit_schur_complement does not support "
+          "use_power_series_expansion_initialization.";
+      return false;
+    }
+  }
+
+  if (options.use_power_series_expansion_initialization ||
+      options.preconditioner_type == SCHUR_POWER_SERIES_EXPANSION) {
+    OPTION_GE(max_num_spse_iterations, 1)
+    OPTION_GE(spse_tolerance, 0.0)
   }
 
   if (options.use_mixed_precision_solves) {
diff --git a/internal/ceres/trust_region_preprocessor.cc b/internal/ceres/trust_region_preprocessor.cc
index a5a986f..e9406d1 100644
--- a/internal/ceres/trust_region_preprocessor.cc
+++ b/internal/ceres/trust_region_preprocessor.cc
@@ -199,6 +199,11 @@
       options.max_linear_solver_iterations;
   pp->linear_solver_options.type = options.linear_solver_type;
   pp->linear_solver_options.preconditioner_type = options.preconditioner_type;
+  pp->linear_solver_options.use_power_series_expansion_initialization =
+      options.use_power_series_expansion_initialization;
+  pp->linear_solver_options.spse_tolerance = options.spse_tolerance;
+  pp->linear_solver_options.max_num_spse_iterations =
+      options.max_num_spse_iterations;
   pp->linear_solver_options.visibility_clustering_type =
       options.visibility_clustering_type;
   pp->linear_solver_options.sparse_linear_algebra_library_type =
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc
index c0e3355..1d514de 100644
--- a/internal/ceres/types.cc
+++ b/internal/ceres/types.cc
@@ -81,6 +81,7 @@
     CASESTR(IDENTITY);
     CASESTR(JACOBI);
     CASESTR(SCHUR_JACOBI);
+    CASESTR(SCHUR_POWER_SERIES_EXPANSION);
     CASESTR(CLUSTER_JACOBI);
     CASESTR(CLUSTER_TRIDIAGONAL);
     CASESTR(SUBSET);
@@ -94,6 +95,7 @@
   STRENUM(IDENTITY);
   STRENUM(JACOBI);
   STRENUM(SCHUR_JACOBI);
+  STRENUM(SCHUR_POWER_SERIES_EXPANSION);
   STRENUM(CLUSTER_JACOBI);
   STRENUM(CLUSTER_TRIDIAGONAL);
   STRENUM(SUBSET);