Background
I'm working on a numerical linear algebra package in C#. I'm trying to implement a variety of "principal pivoting" methods to solve optimization problems (specifically linear complementarity problems, but I think similar schemes are used in linear and quadratic programming). These methods require the ability to multiply a column vector against the "principal pivotal transform" of a matrix.
The principal pivot transform (PPT) of a matrix $M$ is defined as:
$$PM^'P^T = \left (\begin{array}{cc} M_{\alpha\alpha}^{-1} &\ -M_{\alpha\alpha}^{-1} M_{\alpha\beta} \\ M_{\beta\alpha}M_{\alpha\alpha}^{-1} &\ M_{\beta\beta} - M_{\beta\alpha} M_{\alpha\alpha}^{-1}M_{\alpha\beta} \end{array} \right )$$
Where $\alpha$ is the set of pivoted row/column indices, and $\beta$ is the set of unpivoted row/column indices. $P$ is a permutation matrix that keeps the pivoted and unpivoted rows/columns together to make the block matrix form above cleaner. $M_{\alpha\beta}$ is a submatrix of $M$ with the rows in $\alpha$ and the columns in $\beta$, and the other submatrices of $M$ are defined similarly. Also, $M_{\alpha\alpha}^{-1} = (M_{\alpha\alpha})^{-1}$ and not $(M^{-1})_{\alpha\alpha}$
...
Each iteration of a principal pivot method swaps a single index from $\alpha$ to $\beta$ or vice-verse. This single pivoting operation can be expressed as a low rank update of the matrix. Say we're pivoting index $r$ and using $M^'$ to represent the updated matrix:
$$\begin{matrix} M_{rr}^' & = & & \frac{1}{M_{rr}} & \\ \\ \\ M_{ir}^' & = & & \frac{M_{ir}}{M_{rr}}, & (i \neq r) \\ \\ \\ M_{ri}^' & = & & \frac{-M_{ri}}{M_{rr}}, & (i \neq r) \\ \\ \\ M_{ij}^' & = & M_{ij} - & \frac{M_{ir} * M_{ri}}{M_{rr}}, & (i \neq r, j \neq r) \end{matrix} $$
This is what I'm currently doing. Unfortunately, this is basically calculating $M_{\alpha\alpha}^{-1}$ using something like Gaussian elimination. I have a few test cases where the results are suffering from numerical problems even on on small problems, even using double precision floating point. To help combat this, if I remove an index from $\alpha$ (that is, if I unpivot an index), I recalculate the matrix from scratch using one pivot at a time with the above scheme. But this makes most of the algorithms quite slow.
Question
I know there are reasonable ways to calculate the product $M_{\alpha\alpha}^{-1} z$ for some column vector $z$ without calculating $M_{\alpha\alpha}^{-1}$ directly (such as LU decomposition or QR decomposition). I could calculate this decomposition fresh every iteration from the current pivot set ($\alpha$). Which means each iteration of a pivoting algorithm would take $O(n^3)$. That's too slow. I could do a low rank update of the LU or QR decomposition, but since most of the algorithms pivot a single index at a time, I lose most of the benefit of these decomposition schemes after a few iterations. I could perform the full decomposition every few iterations, but that again would be too slow.
Is there something more clever that can be done? I think I basically need to find the decomposition of principal submatrices. Is there any literature on the subject? Any other ways to approach the problem?
