Incomplete LU Preconditioning and Dropped Fill

An iterative solver converges at a rate set by the spectrum of the matrix it is given. The cheapest way to change that spectrum is to multiply by an approximation of the inverse, and the cheapest useful approximation of the inverse is a factorization that has been deliberately prevented from filling in.

The result is a preconditioner: a matrix M that is close enough to A for M^{-1} A to have a clustered spectrum, and structured enough that applying M^{-1} costs about as much as one matrix-vector product. Incomplete LU is the most widely used construction of that kind, and its central decision is which entries to throw away.

The factorization with a hole in it

A complete LU factorization of a sparse matrix produces L and U whose combined pattern is much larger than that of A. Incomplete LU keeps the factors inside a prescribed pattern S, which is usually the pattern of A itself, and discards every update that would land outside it. The discarded updates are not applied at all, so the factors satisfy L U = A + E, where E has nonzeros only in positions outside S.

The computation is the same as the dense one with the inner loops restricted to the stored pattern. The multipliers are formed from the current values, and the update a_ij -= a_ik * a_kj runs only over the j for which row i has a stored entry. Nothing is allocated during the factorization, and the storage is exactly the storage of A.

/* ILU(0) in place on a CSR matrix with sorted column indices.
   On return values[] holds L below the diagonal and U on and above it. */
int ilu0(int n, const int *row_ptr, const int *col_ind, double *values)
{
    for (int i = 0; i < n; i++) {
        for (int kk = row_ptr[i]; kk < row_ptr[i + 1]; kk++) {
            int k = col_ind[kk];
            if (k >= i) break;                       /* only strictly lower part */
            if (values[kk] == 0.0) continue;
            int dk = -1;                             /* find diagonal of row k */
            for (int p = row_ptr[k]; p < row_ptr[k + 1]; p++)
                if (col_ind[p] == k) { dk = p; break; }
            if (dk < 0 || values[dk] == 0.0) return -1;   /* zero pivot */

            double lik = values[kk] / values[dk];
            values[kk] = lik;
            for (int p = dk + 1; p < row_ptr[k + 1]; p++) {
                int j = col_ind[p];
                for (int q = kk + 1; q < row_ptr[i + 1] && col_ind[q] <= j; q++) {
                    if (col_ind[q] == j) {           /* entry exists in the pattern */
                        values[q] -= lik * values[p];
                        break;
                    }
                }
            }
        }
    }
    return 0;
}

The inner search for the position of column j in row i is what makes this version slower than it needs to be. With sorted indices it is a merge, and the standard implementation walks two pointers instead of restarting the scan. With a precomputed mapping from (i, j) to a position, the whole inner loop becomes a sequence of indexed updates, at the cost of an extra array of size nnz.

Why throwing away fill still works

The dropped entries are exactly the terms that make the true factorization accurate. Discarding them means M = L U is not close to A in any norm that would satisfy a numerical analyst, and yet the preconditioned matrix M^{-1} A = I - M^{-1} E can still have a spectrum clustered near one.

The reason is that E is not an arbitrary perturbation. The dropped updates are products of a multiplier and an entry of U, and those products are small when the matrix is diagonally dominant: the multipliers are bounded by the ratio of off-diagonal to diagonal entries, and the entries of U are of the order of the diagonal. For a matrix that is an M-matrix, meaning symmetric positive definite with nonpositive off-diagonals and a positive diagonal, the incomplete factorization exists for every pattern and the multipliers stay bounded. This is the Meijerink and van der Vorst result, and it is the theoretical foundation of the whole family.

For such matrices, the preconditioned operator is close to the identity in the sense that ||I - M^{-1} A|| is bounded below one in an appropriate norm, and the iteration converges. The bound is not sharp and it does not extend to general matrices, but it explains the empirical observation that incomplete factorization works well on discretized diffusion problems and poorly on strongly convection-dominated ones.

When the matrix is not an M-matrix, the guarantee is gone. The multipliers can grow, the dropped terms can be large, and the preconditioned spectrum can be worse conditioned than the original. The construction still runs; it just may not help.

Levels of fill and threshold dropping

ILU(0) is the most aggressive form: no fill at all beyond the original pattern. The family generalizes by assigning a level to each position, with the original nonzeros at level zero and a fill entry at level i + j + 1 when it is created by the product of entries at levels i and j. ILU(k) keeps all positions up to level k. The pattern grows rapidly with k, and the cost of each increment compounds on a two-dimensional problem, so levels beyond the first few become impractical.

The level scheme is structural and ignores magnitude, which is its weakness. A fill entry created by a product of two large numbers is kept at the same level as one created by two small numbers, and the pattern is decided before any values are known. For matrices with strongly varying coefficients, the levels do not reflect importance.

Threshold-based dropping addresses this. ILUT computes the factorization row by row, generating fill as a complete factorization would, and then discarding entries whose magnitude is below a threshold relative to the norm of the row, keeping at most a fixed number of entries per row. Both parameters matter: the drop tolerance controls which entries are negligible, and the fill limit controls memory. The resulting factor has a predictable storage bound, which the level-based scheme does not, and its quality is usually better for the same memory.

The price is that the pattern is no longer known in advance, so the factorization cannot be done in place and the data structures must grow. Threshold dropping also interacts badly with pivoting: choosing a large entry for stability may move a row, which changes the pattern of the rows below it. Implementations handle this with a limited form of column permutation and accept that the result is a preconditioner, not a factorization with a clean error bound.

The triangular solve is the bottleneck

Applying the preconditioner means solving L y = r and U z = y. Each solve is a forward or backward substitution, and the dependency structure is a chain: the i-th entry of y depends on all previous entries with a nonzero in row i. The parallelism available within a single solve is limited to the entries that are mutually independent, which for a banded matrix is the width of the band, and for an irregular matrix can be very small.

Level scheduling exposes what parallelism exists. The rows are grouped into levels such that no row in a level depends on another row in the same level, and the groups are processed in sequence with the rows within a group processed in parallel. For a matrix from a two-dimensional grid ordered by nested dissection, the number of levels is proportional to the separator depth, and the average group size is small at the top of the tree. The speedup is real but bounded, and it degrades as the problem becomes more irregular.

This is the fundamental disadvantage of incomplete factorization relative to preconditioners that are applied by a matrix-vector product. A polynomial preconditioner, a sparse approximate inverse, or a multigrid cycle can be applied with the same parallelism as the product itself. An incomplete factorization cannot, because the triangular solve is inherently sequential along the elimination order.

The standard compromise partitions the matrix into blocks and applies an independent incomplete factorization within each block, with the coupling between blocks handled by the outer iteration. This is block-Jacobi preconditioning with an ILU subdomain solve. It parallelizes trivially, since the blocks are independent, and it weakens as the number of blocks grows, because the preconditioner no longer captures the long-range coupling. The number of iterations rises with the block count, and there is a crossover beyond which the parallel speedup is consumed by the extra iterations.

Breakdown and the pivot that disappears

An incomplete factorization can fail on a matrix that is perfectly nonsingular. The dropped updates are exactly the terms that would have kept the pivot away from zero, and removing them can leave a zero or negative pivot in a position where the complete factorization had a comfortable one. For an M-matrix this cannot happen, by the existence theorem; for a general matrix it happens routinely.

The usual remedy is a diagonal shift: replace the zero pivot by a small positive value, or add a small multiple of the identity to the whole matrix before factorizing. The shifted preconditioner is not a factorization of anything in particular, but it is defined, and the iteration either converges with it or does not. A shift that is too large destroys the approximation; one that is too small does not fix the breakdown. The choice is usually made adaptively, by increasing the shift until the factorization completes.

A second failure is a preconditioner that is accurate but ill conditioned. If M is nearly singular, applying M^{-1} amplifies the residual, and the iteration can diverge even though the preconditioned matrix has a clustered spectrum on paper. The triangular solves are where this shows up: a tiny diagonal entry in U produces a huge entry in the solution, and the subsequent iterations inherit the error. Estimating the condition number of the factors, by the same kind of norm estimation used for a direct solve, is the diagnostic.

The economics of the setup

The incomplete factorization is computed once and applied once per iteration. Its cost is therefore amortized over the whole solve, and the trade is between setup time and iteration count. A denser factor costs more to compute and more to apply per iteration, but reduces the number of iterations. There is a minimum in total time, and it usually sits at a fill level far below the complete factorization.

The trade shifts when the matrix is factorized many times with the same pattern. In a nonlinear solver, the Jacobian is reassembled at each Newton step but its sparsity pattern is fixed, so the symbolic structure and the allocation can be reused while the values are recomputed. In that setting a more expensive preconditioner can be worth its setup cost, because the setup is a smaller fraction of the total.

The alternative is a direct factorization, which has no iteration count to tune and no convergence to diagnose, at the cost of memory that grows with fill. The choice between them is a memory question more than an arithmetic one, and the fill analysis in Fill-In, Elimination Ordering, and Symbolic Factorization is what predicts it. The direct method’s own structure, which reuses dense kernels on small fronts, is described in Multifrontal Factorization and the Dense Front Matrix.