User jay lemmon - MathOverflow

Linear complementarity problem: principal pivoting algorithm

2011-10-06T19:43:22Z

I'm trying to implement the "Dantzig; van de Panne and Whinston" principal pivoting algorithm for solving symmetric positive semi-definite LCPs from "The Linear Complementarity Problem" book (Algorithm 4.2.11 on page 248 of my edition).

The algorithm looks like this:

Given $(q,M)$ with $M$ being square, symmetric and positive semi-definite and of size $n$, try to find $z$ s.t. $Mz+q = w$ and $z\perp w, w \ge 0 \space z \ge 0$, or indicate if no such $z$ exists.

...

Initialize:

$\alpha = \emptyset$, $\beta = \{ 1, ..., n\}$

Step 1: Breaking ties arbitrarily, find: $r \in arg \space min \{ q_i : i \in \beta \}$

If $q_r \ge 0$, stop. A solution of $(q,M)$ is given by $z_\alpha = q_\alpha, z_\beta = 0$

If $M_{rr} = 0$ and $M_{ir} \ge 0$ for all $i \in \alpha$ stop. There is no solution

Step 2: Use the minimum ratio test to find the index $s \in \alpha \cup \{r\}$ of a blocking variable to $r$. Break ties arbitrarily, but prefer $r$ if there's a choice.

Step 3: Using a simple principal pivot on row/column s to update $M$ and $q$

If $s = r$, transfer $r$ from $\beta$ to $\alpha$ and goto Step 1

If $s \ne r$, transfer $s$ from $\alpha$ to $\beta$ and goto Step 2

...

This algorithm works most of the time, but in cases where I perform $n$ pivots, I'm getting to step 1 and $\beta$ is empty, and I'm not sure how to proceed. Some detective work shows that often one more pivot will arrive at the correct answer, but I have a number of different cases and they all require slight differences once they run out of elements in $\beta$ in order to arrive at the correct answer.

One such case is: $M = \begin{bmatrix} 60 & 26 & 37 \\ 26 & 36 & 29 \\ 37 & 29 & 30 \end{bmatrix}$ and $q = \begin{bmatrix} -1 \\ -1 \\ -1 \end{bmatrix}$

For which I know the answer is: $z = \begin{bmatrix} 0 \\ .004184 \\ .029289 \end{bmatrix}$

The matrix M is symmetric PSD as near as I can tell (all the eigen values are positive), so it should be something that Dantzig's algorithm can process. Stepping through the code the above problem gets pivoted 3 times, and I find that a 4th pivot will arrive at a correct answer.

I'm getting tired of hack and check so I'm wondering if there's something I'm just missing about the algorithm's termination condition, or if this case implies a bug in my implementation.

I can provide definitions for "minimum ratio test", "arg min", and "simple principal pivot" if they aren't familiar terms (though they seem to be pretty universal to all pivoting methods for LCPs).

Pseudoinverse of columns of a matrix

2011-01-18T19:10:51Z

First, some background:

I'm working on an implementation in C# of Lemke's algorithm (for solving linear complementarity problems) based on this Matlab implementation: http://ftp.cs.wisc.edu/math-prog/matlab/lemke.m

The implementation calculates $d = B / b_e$ in a tight inner loop. (Here $/$ is the Matlab "backslash" command, which calculates the Moore-Penrose pseudoinverse least squares "solution"). So as the implementation has it I'd have to decompose a matrix every iteration to compute the answer, which would be way to slow for practical use.

However, looking at it closer, the $B$ matrix is actually a permutation of a subset of the columns of a larger matrix which doesn't change from iteration to iteration, which has the form (in block notation) $[M, -I, c]$, where $M$ is a supplied square matrix, $I$ is the identity matrix of the same size as $M$, and $c$ is a calculable column vector. I'll call this larger matrix $B_s$, and $B_s P = B$, where $P$ is like a permutation matrix except that it actually excludes some columns (so that the result, $B$, is square).

As for $b_e$, it's actually just a column from $B_s$.

Which means that I can effectively pre-calculate all possible values for $d$ before I start the loop. If I let $B_s^{+}$ represent the Moore-Penrose Pseudoinverse of $B_s$ (likewise with $P$), I get:

$d = P^+ B_s^+ b_e$, where $b_e = $ column $e$ from $B_s$.

Now the question

Given a matrix $B_s$ and it's Moore-Penrose pseudoinverse $B_s^+$, is there any special significance to calculating $B_s^+b_e$ where $b_e$ is a column of $B_s$?

If $B_s$ were square and of full rank, the pseudoinverse would just be the normal matrix inverse. In which case:

$b_e = B_s e$ (where $e$ is a column vector with only one entry set to $1$ and all others set to $0$).
$B_s^{+} b_e = B_s^{-1} (B_s e) = (B_s^{-1} B_s) e = (I) e = e$

However, I don't think that's true for the general pseudoinverse of rectangular, possibly rank deficient, matrices.

Bounding the roots of the sum of two polynomials

2010-10-08T00:51:44Z

Suppose I have two polynomials with real coefficients. Suppose I can perform any sort of preprocessing on them I want. I want to be able to pre-emptively say that the sum of the polynomials doesn't have any roots inside a given interval without doing any explicit calculations on the sum itself. False positives (that is, saying there aren't any roots when there are some) would be deal-breaking, but false negatives (reporting there might be roots when there aren't) would be acceptable.

Or to put it more explicitly:

All functions $p_x(t)$ have a form like:

$p_x(t) = a_{n,x} * t^n + a_{n-1, x} * t^{n-1} + ... + a_{1,x} * t + a_{0,x}$

We can define $p_3(t) = p_1(t) + p_2(t)$

I want to determine if $p_3(t)$ might have any roots inside a given interval $[t_{min}, t_{max}]$. But I want to do it only using properties of $p_1(t)$ and $p_2(t)$, their roots, etc. and not anything that would need me to calculate anything for $p_3(t)$, its roots, etc.

Any ideas on how to approach the problem?

EDIT: So some motivation of what I'm doing: I have a large set of polynomials that are related to the path of a point through space over time. I want to find polynomials that intersect sometime in the "near" future, but I don't want to have to do all $\frac{n*(n-1)}{2}$ polynomial-to-polynomial evaluations. So I'm trying to build a "broad phase" that only offers up pairs of polynomials to be solved in a "narrow phase" (ie: actual root finding) if they're "pretty close" to colliding. Whatever the algorithm for the broad phase is, it can't involve iterating over all the polynomial pairs or it defeats the point.

One sort of square-peg-round-hole solution would be to use something like bounding boxes around the polynomials and use a spatial partitioning tree to find where boxes overlap, and then do the root finding on those. But it doesn't handle cases very well where the time interval of interest is quite large, or especially if one of the interval ends is infinity or negative infinity.

So I wanted to explore it from another direction and see if I can come up with something that works better.

Comment by Jay Lemmon

2011-10-06T20:31:38Z

I'm not asking to find the bug, I'm asking how you're supposed to handle a specific edge case (the case where you perform $n$ "major cycles" and $\beta = \nullset$) of this algorithm. The algorithm as written in the book doesn't handle it. The proof of termination doesn't seem to either, at least not directly. I don't have an instructor; I'm not in university. Non-students do read textbooks, on occasion.

Comment by Jay Lemmon

2011-10-06T19:47:32Z

From the FAQ: "MathOverflow's primary goal is for users to ask and answer research level math questions, the sorts of questions you come across when you're writing or reading articles or graduate level books." The book sure seems graduate level to me.

Comment by Jay Lemmon

2011-01-19T19:02:05Z

Ah actually I answered my own question. I can't distribute the pseudoinverse from $B^+ = (B_s P)^+$ since $P P^T \not= I$. Back to the drawing board :/

Comment by Jay Lemmon

2011-01-19T17:44:31Z

Why would this be inappropriate if I can essentially cache everything before I start pivoting? Or is there some reason I can't?

Comment by Jay Lemmon

2010-10-08T17:40:08Z

This looks interesting, but I don't follow well enough. Using something like Cauchy bounds you always end up with an upper bound on the absolute value of possible roots $|t|$ that is greater than 1. So it can't ever be $<1/a$. And I definitely don't see where you're building $Q(t)$ from.

Comment by Jay Lemmon

2010-10-08T17:18:55Z

Your blockquote solution would be the naive way, yeah. Anything to improve on that would be good. I'll play with the non-blockquote solution you gave and see if it works for what I have in mind.

Comment by Jay Lemmon

2010-10-08T17:02:54Z

Yep, you guessed the motivation. Cookie for you ;)