GitHub

Frequently Asked Questions

How is LinearSolve.jl compared to just using normal \, i.e. A\b?

Check out this video from JuliaCon 2022 which goes into detail on how and why LinearSolve.jl can be a more general and efficient interface.

Note that if \ is good enough for you, great! We still tend to use \ in the REPL all the time! However, if you're building a package, you may want to consider using LinearSolve.jl for the improved efficiency and ability to choose solvers.

I'm seeing some dynamic dispatches in the default algorithm choice, how do I reduce that?

Make sure you set the OperatorAssumptions to get the full performance, especially the issquare choice as otherwise that will need to be determined at runtime.

I found a faster algorithm that can be used than what LinearSolve.jl chose?

What assumptions are made as part of your method? If your method only works on well-conditioned operators, then make sure you set the WellConditioned assumption in the assumptions. See the OperatorAssumptions page for more details. If using the right assumptions does not improve the performance to the expected state, please open an issue and we will improve the default algorithm.

Python's NumPy/SciPy just calls fast Fortran/C code, why would LinearSolve.jl be any better?

This is addressed in the JuliaCon 2022 video. This happens in a few ways:

  1. The Fortran/C code that NumPy/SciPy uses is actually slow. It's OpenBLAS, a library developed in part by the Julia Lab back in 2012 as a fast open source BLAS implementation. Many open source environments now use this build, including many R distributions. However, the Julia Lab has greatly improved its ability to generate optimized SIMD in platform-specific ways. This, and improved multithreading support (OpenBLAS's multithreading is rather slow), has led to pure Julia-based BLAS implementations which the lab now works on. This includes RecursiveFactorization.jl which generally outperforms OpenBLAS by 2x-10x depending on the platform. It even outperforms MKL for small matrices (<100). LinearSolve.jl uses RecursiveFactorization.jl by default sometimes, but switches to BLAS when it would be faster (in a platform and matrix-specific way).
  2. Standard approaches to handling linear solves re-allocate the pivoting vector each time. This leads to GC pauses that can slow down calculations. LinearSolve.jl has proper caches for fully preallocated no-GC workflows.
  3. LinearSolve.jl makes many other optimizations, like factorization reuse and symbolic factorization reuse, automatic. Many of these optimizations are not even possible from the high-level APIs of things like Python's major libraries and MATLAB.
  4. LinearSolve.jl has a much more extensive set of sparse matrix solvers, which is why you see a major difference (2x-10x) for sparse matrices. Which sparse matrix solver between KLU, UMFPACK, Pardiso, etc. is optimal depends a lot on matrix sizes, sparsity patterns, and threading overheads. LinearSolve.jl's heuristics handle these kinds of issues.

How do I use IterativeSolvers solvers with a weighted tolerance vector?

IterativeSolvers.jl computes the norm after the application of the left preconditioner. Thus, in order to use a vector tolerance weights, one can mathematically hack the system via the following formulation:

import LinearSolve as LS
import LinearAlgebra as LA

n = 2
A = rand(n, n)
b = rand(n)

weights = [1e-1, 1]
precs = Returns((LS.InvPreconditioner(LA.Diagonal(weights)), LA.Diagonal(weights)))

prob = LS.LinearProblem(A, b)
sol = LS.solve(prob, LS.KrylovJL_GMRES(precs))

sol.u
2-element Vector{Float64}:
 -0.4286154920819371
  1.1907995115116115

If you want to use a “real” preconditioner under the norm weights, then one can use ComposePreconditioner to apply the preconditioner after the application of the weights like as follows:

import LinearSolve as LS
import LinearAlgebra as LA

n = 4
A = rand(n, n)
b = rand(n)

weights = rand(n)
realprec = LA.lu(rand(n, n)) # some random preconditioner
Pl = LS.ComposePreconditioner(LS.InvPreconditioner(LA.Diagonal(weights)),
    realprec)
Pr = LA.Diagonal(weights)

prob = LS.LinearProblem(A, b)
sol = LS.solve(prob, LS.KrylovJL_GMRES(precs = Returns((Pl, Pr))))
retcode: Success
u: 4-element Vector{Float64}:
  1.3509988634214545
  0.6649561611284609
 -1.4783106990272532
  0.5115462345772893

Can I use a Krylov solver with a custom array type?

Yes, if similar works for it.

Krylov.jl allocates its workspace as S(undef, n) with S = typeof(b), so a type that cannot be built that way never reaches an iteration. RecursiveArrayTools.ArrayPartition is the common example: it is stored as several separate arrays, so there is no way to know how to split n across them.

For these, LinearSolve builds the workspace with Krylov.KrylovConstructor, which uses similar on the vectors it is given rather than the undef constructor. The solve then runs on the caller's own array type, with no flattening and no copying, and u comes back in the same type it went in as.

ArrayPartition works out of the box. Another array type needs similar and the usual array operations Krylov uses.

A right-hand side that is not an array at all, such as a parameter object, is not supported yet: it has neither the undef constructor nor similar. The intended general answer there is the SciMLStructures.jl interface, since a type implementing it can be canonicalized to a flat buffer and repacked afterwards, which works for containers that are not arrays. Tracked in #1208.

Why does LinearSolve.jl depend on MKL_jll, and how do I stop it from loading?

MKL is the fastest BLAS on most x86 hardware, often by a wide margin, so LinearSolve.jl ships it by default and lets the default algorithm choice pick MKLLUFactorization where it wins. Without it, most installations end up substantially slower.

If you would rather not load it, for example because you only solve small static-array systems where it brings nothing, set the LoadMKL_JLL preference to false:

using Preferences, UUIDs

Preferences.set_preferences!(
    UUID("7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"),  # LinearSolve
    "LoadMKL_JLL" => false; force = true
)

The preference is read when LinearSolve.jl loads, so restart Julia afterwards and let the package recompile. From then on LinearSolve.usemkl is false, MKL_jll is never using'd, and the default algorithm falls back to the next best choice for your matrix, typically LUFactorization or RFLUFactorization. Nothing else about the interface changes and every algorithm you select explicitly keeps working.

Two details worth knowing:

  • The default is already architecture aware. MKL is only considered on x86_64 and i686, and it is off by default on AMD EPYC CPUs, where it does not win. On other architectures such as Apple Silicon it is never loaded regardless of this setting.
  • The preference controls whether LinearSolve.jl loads and uses MKLjll, not whether it is installed. MKLjll stays a declared dependency, so it remains in the dependency graph and Pkg still installs it.

Why does differentiating a solve with Enzyme fail on cache.u?

Take the solution from the value solve! gives back, not from the cache:

sol = solve!(cache)
sum(sol.u)        # differentiable
sum(cache.u)      # errors, "Adjoint case currently not handled"

The reverse rule is written against the returned solution, so reading cache.u after the solve has no derivative attached to it.

Mixing activities is fine for a dense or sparse A: carrying a derivative on A while b is a captured constant, or the other way round, differentiates under a plain Reverse. A structured A, such as a Tridiagonal or a unit triangular, still wants set_runtime_activity for that case.