3
$\begingroup$

It was suggested I should try posting this question here from Mathematics

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 the Principal Pivot algorithm I have 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$), and probably get a pretty robust method. But then each iteration of the 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 only every few iterations, but then the question becomes how often should I do so? And it still would be quite 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?

$\endgroup$
  • $\begingroup$ Is $M$ dense, positive-definite and/or symmetric in the general case? $\endgroup$ – Aron Ahmadia Apr 23 '12 at 7:53
  • $\begingroup$ In my case it's symmetric and PSD. It's actually formed from the normal equations of some rank deficient dense matrices. I know there are some iterative schemes that might work well here, but I wanted to explore pivoting methods. $\endgroup$ – Jay Lemmon Apr 23 '12 at 16:26
  • 1
    $\begingroup$ If it's SPD, then you can take the Cholesky decomposition, and there is a fair amount of literature on how to stably update the Cholesky decomposition using low-rank factors as you are suggesting. This is very commonly used in implementations of simplectical solvers, so I would not be surprised if it turns out to be the most efficient for you. $\endgroup$ – Aron Ahmadia Apr 24 '12 at 7:51
  • $\begingroup$ Ah, that's actually an interesting idea. I'll look in to it. $\endgroup$ – Jay Lemmon Apr 24 '12 at 20:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.