diff --git a/CHANGELOG.md b/CHANGELOG.md
index bac863ac9..25e8f3888 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
 - Fix CUDA `sparse_to_dense`. ([#386])
+- Improve pseudo inverse spectrum solver. ([#388])
 
 ## [v0.25.2]
 Release date: 2025-02-02
@@ -107,3 +108,4 @@ Release date: 2024-11-13
 [#380]: https://github.com/qutip/QuantumToolbox.jl/issues/380
 [#383]: https://github.com/qutip/QuantumToolbox.jl/issues/383
 [#386]: https://github.com/qutip/QuantumToolbox.jl/issues/386
+[#388]: https://github.com/qutip/QuantumToolbox.jl/issues/388
diff --git a/src/spectrum.jl b/src/spectrum.jl
index 77b690976..14bc88893 100644
--- a/src/spectrum.jl
+++ b/src/spectrum.jl
@@ -125,19 +125,20 @@ function _spectrum(
     _tr = SparseVector(D^2, [1 + n * (D + 1) for n in 0:(D-1)], ones(_CType(L), D)) # same as vec(system_identity_matrix)
     _tr_A = transpose(_tr) * spre(A).data
 
-    cache = nothing
-    I_cache = I(D^2)
+    Id = I(D^2)
+
+    # DO the idx = 1 case
+    ω = ωList[1]
+    cache = init(LinearProblem(L.data - 1im * ω * Id, b), solver.alg, kwargs...)
+    sol = solve!(cache)
+    spec[1] = -2 * real(dot(_tr_A, sol.u))
+    popfirst!(ωList)
     for (idx, ω) in enumerate(ωList)
-        if idx == 1
-            cache = init(LinearProblem(L.data - 1im * ω * I_cache, b), solver.alg, kwargs...)
-            sol = solve!(cache)
-        else
-            cache.A = L.data - 1im * ω * I_cache
-            sol = solve!(cache)
-        end
+        cache.A = L.data - 1im * ω * Id
+        sol = solve!(cache)
 
         # trace over the Hilbert space of system (expectation value)
-        spec[idx] = -2 * real(dot(_tr_A, sol.u))
+        spec[idx+1] = -2 * real(dot(_tr_A, sol.u))
     end
 
     return spec