ISAI (Incomplete Sparse Approximate Inverse)#

gko::preconditioner::Isai<IsaiType, ValueType, IndexType> computes a sparse approximation \(M_A \approx A^{-1}\) directly, then applies it as the preconditioner via a single sparse matrix-vector product. The apply step stays free of triangular solves, which is exactly the appeal on hardware where triangular substitution serialises poorly.

The trade-off is that the quality of the preconditioner is bounded by the sparsity pattern it is constrained to. ISAI fixes the pattern of \(M_A\) to match that of \(A\) (or one of its triangular factors), then solves a collection of small dense systems to choose the entries that minimise the local error.

Four ISAI flavours are selectable through the IsaiType template parameter:

`IsaiType`	Approximates	Typical pairing
`isai_type::general`	\(M_A \approx A^{-1}\)	Standalone preconditioner for general \(A\) where a triangular factorisation is unavailable.
`isai_type::lower`	\(M_L \approx L^{-1}\)	Replaces the lower-triangular solve in an ILU / IC pipeline.
`isai_type::upper`	\(M_U \approx U^{-1}\)	Replaces the upper-triangular solve in an ILU pipeline.
`isai_type::spd`	\(M_C^T M_C \approx C^{-1} C^{-T}\) for the Cholesky factor \(C\) of \(A\)	Symmetric variant: the Factorized Sparse Approximate Inverse (FSPAI).

When to use ISAI#

As a replacement for the triangular solves inside an ILU / IC preconditioner when those solves serialise too poorly on the target hardware. Build the factorisation as usual, then plug Isai<lower> and Isai<upper> in as the l_solver / u_solver factories.
As a standalone preconditioner for general \(A\) when no factorisation is available; ISAI’s per-iteration cost is just an SpMV.
For SPD systems (spd flavour) where a triangular-solve-free preconditioner is preferred for parallelism reasons.

ISAI’s preconditioner quality is fundamentally limited by the chosen sparsity pattern; for matrices where the inverse is dense, ISAI cannot compete with a full factorisation regardless of how much computation you throw at the generation step. Profile against IC / ILU with a Krylov inner triangular solve before adopting ISAI.

Construction#

// Standalone general ISAI preconditioner
auto isai = gko::preconditioner::Isai<
                gko::preconditioner::isai_type::general,
                double, gko::int32>::build()
                .on(exec);

auto preconditioner = isai->generate(system_matrix);

// As the triangular solvers inside an ILU preconditioner
auto ilu = gko::preconditioner::Ilu<double, false, gko::int32>::build()
               .with_factorization(
                   gko::factorization::ParIlu<double, gko::int32>::build())
               .with_l_solver(
                   gko::preconditioner::LowerIsai<double, gko::int32>::build())
               .with_u_solver(
                   gko::preconditioner::UpperIsai<double, gko::int32>::build())
               .on(exec);

LowerIsai, UpperIsai, GeneralIsai and SpdIsai are aliases for Isai<isai_type::lower|upper|general|spd, ValueType, IndexType>. Pass LowerIsai / UpperIsai factories to Ilu::with_l_solver / with_u_solver, or a LowerIsai factory to Ic::with_l_solver.

Factory parameters#

Parameter	Type	Purpose
`skip_sorting`	`bool`	Set `true` when the input matrix is already sorted by row and column index to skip the internal sort. Default `false`.
`sparsity_power`	`int`	Use the sparsity pattern of \(A^k\) instead of \(A\) for the approximate inverse. Larger values give denser approximations and better preconditioner quality at the cost of generation time and storage. Default `1`.
`excess_limit`	`size_type`	Rows with more than 32 nonzeros are handled through an “excess” system; when set to \(m > 0\), that system is solved as soon as its size exceeds \(m\), repeatedly, until all excess rows are done. Default `0` (one excess solve).
`excess_solver_factory`	`LinOpFactory`	Solver for the excess system. Defaults to a triangular solver for `lower` / `upper` ISAI and to block-Jacobi-preconditioned GMRES for `general` / `spd` ISAI.
`excess_solver_reduction`	`remove_complex<ValueType>`	Residual-reduction factor of the default GMRES excess solver used for `general` and `spd` ISAI. Only the default GMRES excess solver (`general` / `spd`) reads it. Default `1e-6`.

Implementation notes#

Each row of the approximate inverse is computed from a small local system. Rows with more than 32 nonzeros are handled through the excess system controlled by excess_limit and excess_solver_factory; with sparsity_power > 1 more rows reach that size, which raises generation time.

In the general case, \(\mathrm{ISAI}(A)^T \ne \mathrm{ISAI}(A^T)\); only the spd flavour is transpose-symmetric. If you need the transposed approximate inverse, compute it from \(A^T\) directly.