Dense and selective matrix inverses
An existing LU factorization can be used either to compute every entry of a matrix inverse or only selected entries. Lowercase italic symbols denote scalars. Bold lowercase symbols denote vectors. Bold uppercase symbols denote matrices, including dense blocks.
Dense inverse
For a dense matrix \(\mathbf{A}=\mathbf{L}\mathbf{U}\), its inverse \(\mathbf{X}=\mathbf{A}^{-1}\) satisfies
Substituting the LU factorization gives:
Let \(\mathbf{Y}=\mathbf{U}\mathbf{X}\). First, solve the lower-triangular system by forward substitution:
Then solve the upper-triangular system by backward substitution:
For every column \(\mathbf{e}_j\) of \(\mathbf{I}\), this means first solving \(\mathbf{L}\mathbf{y}_j=\mathbf{e}_j\) by forward substitution and then \(\mathbf{U}\mathbf{x}_j=\mathbf{y}_j\) by backward substitution. The resulting columns \(\mathbf{x}_j\) form the dense inverse \(\mathbf{X}\). In the block-sparse algorithm, the same dense procedure is applied locally to obtain the diagonal contribution \((\mathbf{L}_p\mathbf{U}_p)^{-1}\) of each pivot block.
Selective inverse
The selective inverse computes entries of \(\mathbf{A}^{-1}\) matching the sparsity pattern of an existing LU factor. We use the Takahashi equations for selective inverse in PGM. It was introduced by Takahashi[1] and further described by Erisman and Tinney[2].
The following sections derive the Takahashi equations for scalar-sparse and block-sparse matrices. They then describe their numerical implementation in PGM.
Scalar-sparse matrices
Let \(\mathbf{A}=\mathbf{L}\mathbf{U}\), where \(\mathbf{L}\) is unit lower triangular and \(\mathbf{U}\) is upper triangular, and define
The identities
give the inverse entries at pivot \(p\). Because \((\mathbf{L}^{-1})_{pj}=0\) for \(j>p\),
Similarly, \((\mathbf{U}^{-1})_{kp}=0\) for \(k>p\) and \(l_{pp}=1\), hence
Finally, \((\mathbf{L}^{-1})_{pp}=1\) gives
Intuitively, an upper-triangular entry of \(\mathbf{Z}\) depends on previously computed entries \(z_{mj}(m>p)\) below it. An entry in the lower triangular part depends on previously computed entries \(z_{km}(m>p)\) to its right. We call these dependencies “looking down” and “looking right,” respectively.
These equations are evaluated for \(p=n-1,\ldots,0\). The sums only visit structurally non-zero entries in the filled LU pattern.
Block-sparse matrices
Now consider block sparse matrices, where each non-zero entry is a dense block of size \(b\times b\).
The inverse is again \(\mathbf{Z}=\mathbf{U}^{-1}\mathbf{L}^{-1}\), so the identities
can be expanded block by block.
Because \((\mathbf{L}^{-1})_{pj}=\mathbf{0}\) for \(j>p\),
Similarly, \((\mathbf{U}^{-1})_{kp}=\mathbf{0}\) for \(k>p\), hence
Finally, \((\mathbf{L}^{-1})_{pp}=\mathbf{L}_p^{-1}\) gives
Multiplying from the left by \(\mathbf{U}_p^{-1}\) yields
These are the scalar equations with division replaced by dense inversions.
For block size \(1\times 1\), \(\mathbf{L}_p=[1]\) and \(\mathbf{U}_p=[u_{pp}]\), so the block equations reduce exactly to the scalar equations.
Numerical implementation
Similarities and differences with the LU solver
The numerical implementation consumes the LU factors produced by PGM’s block-sparse LU factorization.
Block-sparse LU factorization and selective inversion have the following similarities:
They use the same filled block-sparse indexing.
They operate in place.
At each pivot step, they both update one pivot block row and block column.
They use the packed dense diagonal factor and triangular solves at each pivot.
The main differences are the direction and result of the sweep:
LU factorization proceeds from the top-left to the bottom-right. At each pivot step, one column of \(\mathbf{L}\) and one row of \(\mathbf{U}\) are computed. This \(\mathbf{L}\) column and \(\mathbf{U}\) row are then used to update the trailing Schur complement.
Selective inversion proceeds from the bottom-right to the top-left. At each pivot step, selected blocks in one column and one row of the inverse are computed. Each pivot step uses previously computed inverse blocks at the bottom-right, which are already final.
Dependency blocks and target blocks
For the use-case of state estimation, the target blocks are often only those in the original y_bus pattern.
In radial networks, the LU pattern and the original y_bus pattern are the same.
In meshed networks, factorization introduces fill-ins.
The dependency block set is then the complete filled LU pattern.
A fill-in cannot be skipped merely because it is not part of the final y_bus target pattern.
Target blocks are extracted only after all pivot updates have completed.
Reverse pivot traversal
Each update at pivot \(p\) refers only to inverse blocks whose row and column indices are greater than \(p\). The sweep therefore starts at the last pivot, i.e, the bottom-rightmost block, and proceeds in reverse order, \(p=n-1,\ldots,0\).
Pivot step order
Each pivot step performs the following operations:
Buffer the packed pivot block, the original \(\mathbf{U}_{pm}\) blocks, and the original \(\mathbf{L}_{mp}\) blocks.
Compute the lower blocks \(\mathbf{Z}_{kp}\).
Compute the upper blocks \(\mathbf{Z}_{pj}\).
Compute \(\mathbf{Z}_{pp}\) last, using the newly computed lower blocks.
This order permits the LU storage to be overwritten in place without losing values needed by the current step.
Why buffers are needed?
Each pivot step overwrites packed LU blocks in data with inverse blocks.
The implementation first buffers every LU value that a later update still needs.
Off-diagonal row and column buffers
u_row copies the contiguous \(\mathbf{U}_{pm}\) blocks.
l_col gathers the scattered \(\mathbf{L}_{mp}\) blocks, and l_indices records their locations in data.
Lower and upper updates write \(\mathbf{Z}_{kp}\) and \(\mathbf{Z}_{pj}\) in place while preserving factors for later sums.
The diagonal update reuses u_row with the newly computed \(\mathbf{Z}_{mp}\) blocks.
Factorized pivot-block buffer
pivot copies the packed \(\mathbf{L}_p\) and \(\mathbf{U}_p\) block and remains factorized throughout the pivot step.
The lower and upper updates use it for triangular solves.
The diagonal update uses it to compute \((\mathbf{L}_p\mathbf{U}_p)^{-1}\).
The pivot location is overwritten with \(\mathbf{Z}_{pp}\) only after these operations finish.
Temporary storage is therefore limited to u_row, l_col, l_indices, and pivot; the selective inverse remains in
data and is produced in place.
Restoring dense block permutations
The sparse LU solver in PGM uses dense full pivoting inside each block. It therefore produces the factorization
where \(\mathbf{P}=\operatorname{blockdiag}(\mathbf{P}_i)\) and \(\mathbf{Q}=\operatorname{blockdiag}(\mathbf{Q}_i)\). The selective inverse sweep operates directly on the stored \(\mathbf{L}\) and \(\mathbf{U}\) factors and computes
on the filled pattern. This is the inverse in the permuted ordering, not yet \(\mathbf{A}^{-1}\). Rearranging the factorization shows that the original ordering is restored by applying \(\mathbf{Q}\) from the left and \(\mathbf{P}\) from the right:
Because \(\mathbf{P}\) and \(\mathbf{Q}\) are block diagonal, this can be done independently for every stored block: