As the sparse matrix for, in particular, the BicubicBSpline grows in size the built in \ operator may fail. From the docstring of the \ operator
When A is sparse, a similar polyalgorithm is used. For indefinite matrices, the LDLt factorization does not use pivoting during the numerical factorization and therefore the procedure can fail even for invertible matrices.
As the matrix for the Bicubic B-spline is indefinite, the lack of pivoting is an issue. The current strategy from #48 is to compute the qr factorization prior to the \ call. This serves to stabilize the computation of the B-splines coefficients, but may not be optimal.
One could add LinearSolve.jl as a dependency which opens the possibility of many sparse solvers. Then we would replace a line like
with
Ax_b = LinearProblem(A, b)
x = solve(Ax_b, selected_solver()).uSome preliminary testing for an approximately 100_000 x 100_000 sparse system for particular choices of selected_solver() are
QRFactorization() perform equivalent to qr(A) \ b, which is not surprising as LinearSolve.jl calls the underlying built in qr factorization from LinearAlgebraUMFPACKFactorization() failsKLUFactorization() works but is about 20x time slower than the QR strategy. The QR takes ≈13.5s and KLU takes ≈242s- Many iterative solver options are available like
KrylovJL_GMRES or KrylovJL_BICGSTAB but testing would need to be done to find a "good" preconditioner and initial guess.
Including the LinearSolve.jl as a dependency gives flexibility with respect to solver selection but also some pre-compilation overhead. Whether or not this is "worth it" remains to be seen and requires more testing. In particular, with respect to Krylov-type methods for very large datasets. However, for modest datasets (up to ≈ 500_000 x 500_000) the qr strategy has proven sufficient.
Finally, below is a spy plot of the sparsity pattern for the approximately 100_000 x 100_000 sparse matrix for the Bicubic B-spline with "not-a-knot" boundary conditions. There is a predominant block-tridiagonal structure but the boundary coupling at the far right may introduce some complexity. Perhaps seeing the structure could help others more experienced with numerical linear algebra guide the choice of a sparse solver.
As the sparse matrix for, in particular, the BicubicBSpline grows in size the built in
\operator may fail. From the docstring of the\operatorAs the matrix for the Bicubic B-spline is indefinite, the lack of pivoting is an issue. The current strategy from #48 is to compute the
qrfactorization prior to the\call. This serves to stabilize the computation of the B-splines coefficients, but may not be optimal.One could add
LinearSolve.jlas a dependency which opens the possibility of many sparse solvers. Then we would replace a line likewith
Some preliminary testing for an approximately
100_000 x 100_000sparse system for particular choices ofselected_solver()areQRFactorization()perform equivalent toqr(A) \ b, which is not surprising asLinearSolve.jlcalls the underlying built inqrfactorization fromLinearAlgebraUMFPACKFactorization()failsKLUFactorization()works but is about 20x time slower than the QR strategy. The QR takes ≈13.5s and KLU takes ≈242sKrylovJL_GMRESorKrylovJL_BICGSTABbut testing would need to be done to find a "good" preconditioner and initial guess.Including the
LinearSolve.jlas a dependency gives flexibility with respect to solver selection but also some pre-compilation overhead. Whether or not this is "worth it" remains to be seen and requires more testing. In particular, with respect to Krylov-type methods for very large datasets. However, for modest datasets (up to ≈500_000 x 500_000) theqrstrategy has proven sufficient.Finally, below is a
spyplot of the sparsity pattern for the approximately100_000 x 100_000sparse matrix for the Bicubic B-spline with "not-a-knot" boundary conditions. There is a predominant block-tridiagonal structure but the boundary coupling at the far right may introduce some complexity. Perhaps seeing the structure could help others more experienced with numerical linear algebra guide the choice of a sparse solver.