bpchesney.org

In an earlier post, we saw the benefits of finding the solution to a set of equations with the most diversity, which was found by minimizing the peak-to-sum ratio (PSR).  This post discusses ways to formulate the linear program to arrive at that solution.

Consider the following system of equations.

The A and b matrices for this system of equations are

We focus on finding an initial solution for the linear program using two methods.

Method 1: Elimination and Back-substitution

The first method for solving for the minimum PSR solution starts with Gaussian elimination and Jordanian back-substitution to find the initial solution. From there, the linear program iterates to find a feasible solution with the minimum PSR, that is also an integer solution, greater than 0.  With elimination and back-substitution, [A b] turns into:

This means that 

is a solution to the system of equations, if non-integer solutions are allowed.  However, in this problem, they are not.  The first method takes this solution and then finds the nearest integer vector — not necessarily a solution — and then proceeds to the minimum PSR solution in the feasible set, considering only integer vectors.

For an underdetermined system of equations, using elimination and substitution will not yield a solution, at all.  Instead it will describe the system of equations in terms of certain variables for which values need to be chosen, called “free variables.”  From the choice of free variables, the others are determined and the linear program can iteratively make better and better choices of free variables.  One common initial choice for the the free variables is to set them all to 0.  This, unfortunately, will yield the most sparse solution, which we saw earlier is the opposite of the most diverse solution.  Starting with the sparsest possible solution means the linear program will have to work for as long as possible to get to the most diverse solution.

 Method 2: The Intersection of a Boundary Condition and the Minimum PSR Line

Let’s modify the system of equations slightly to look at them from a different angle.  We’ll turn the “greater-than” to “less-than” and add in another constraint.  This looks more like the investing example from earlier.

Note also that PSR is bounded for nonnegative integers.

The second method for finding the initial solution for the maximum diversity integer program exploits this fact: the PSR of an N-element, nonnegative, integer vector can never be below 1/N.  So if one of our equations has the form:

The minimum PSR integer program can reduce to:

So you can use a sparsity-promoting integer program to find the minimum PSR, maximum diversity solution once you subtract off the evenly distributed vector. Note that this vector does not need to be in the feasible solution space – it can be used as an initial value. Now the minimum PSR integer program can take advantage of the computational efficiency of and the extensive analysis devoted to sparsity-promoting integer programs.  We are looking for a vector that is not very different from the evenly distributed vector, just like with sparsity-promoting programs we look for a vector that is not very much different from 0.

Making Sure All the Units are the Same is a Big Deal

If N >> S, then subtracting off the evenly distributed vector doesn’t do much and we’re pretty much back to trying to solve a system of equations by minimizing the first norm.  This only works if x is sparse, otherwise, its a terrible approximation to the solution of Ax<=b.  However, N >> S is telling us that there isn’t much to distribute, compared to the number of bins we are distributing to.  So we may not get a good, diverse answer anyway.

Another way of looking at the constraint

is that it relates all of the elements of x to one another.  This equation requires that all of the units of x are the same, for the equation to make any sense.  Referring back to the investing example again, the A that gives the system of equations for this example, overall, is actually pretty sparse.  Without the two row vectors of A that are not mostly made up of zeroes, 1^T and p^T, the possible values of x could take on wildly different values, if the values of b varied wildly.  For instance, x₂ could be on the order of 1000 and x₃ could be on the order of 1,000,000.  This would make it very difficult for an integer program that minimizes the first norm, as given in method 2, to arrive at a good answer.

So the choice of starting with a dense row of A and finding its most diverse solution as an initial vector is not arbitrary.  For well-formed problems that seek a most diverse solution, it can actually be a great starting point.  If your system of equations for a problem does not have such a vector then you may not really need the most diverse solution or maybe you have not fully captured the problem.  Usually, a system of equations has most, if not all dense rows.  Then it becomes a question of finding the best one for finding an initial solution to the maximum diversity integer program, which I’m sure is a whole investigation unto itself.

A great example of an application of Gaussian elimination comes from Dr. Gil Strang’s Introduction to Linear Algebra. In it, he discusses solving (AtCA) for graphs and networks. Let’s use the example he formulates with Ohm’s Law and Kirchoff’s Voltage and Current Laws (KVL and KCL, respectively) to illustrate how Gaussian elimination is used to solve problems. It gives a network of resistors, every node a voltage, every pair of nodes connected by an edge, representing current flow. The matrix representation of this problem is

Find x for a given b. Once you know x, which represents the node voltages, you can find y, the currents between the nodes. The matrix A is determined by the connectivity diagram in the example in the book. It’s given below.  The matrix C represents the conductances between each node and is given below as well.

One method for finding the node voltages and currents given the connectivity and conductance matrices is to use Gaussian elimination.  This page shows one way to visualize the solution to this problem, using HTML5′s canvas element.  Dr. Strang does an excellent job of explaining the solution to this problem in his book.  The method can be summarized as follows:

1. Multiply A on the left by C; call this CA
2. Multiply the result of step 1 on the left by A’ (A transpose)
3. Remove the 4th row and 4th column from the result of step2.  Append a column of [1, 0, 0]‘ appended on what’s left.
4. Perform Gaussian elimination on the result of step 3. The node voltages, x, will appear as the 4th column of the 3×4 matrix on which elimination was performed.
5. Multiply x on the left by CA and scale all entries of the result by -1.  The final result is y, the edge currents.

The purpose of step 3 is to set a common potential, ground, as a reference for all the other nodes.  In doing so, the 4th row becomes a redundant combination of the remaining three unknown node voltages and can be eliminated.  The 4th column can be eliminated as well since the common voltage has been set to ground, which is 0 potential.

In step 4, b is introduced into the solution.  This is done by setting the node voltage x1, to some potential, in this case 1, relative to ground (x4).  If x4 is where ground is connected, then x1 is where the external power supply is connected.  The rest of b (x2=0, x3=0) means that x2 and x3 do not have any external voltage source connected to them.

In Python, this solution can be found using scipy and this script.  Using scipy, Gaussian elimination can be performed by breaking down LU factorization into two steps: the forward elimination step and the back substitution step.  To do this, the user must use the “lu_factor” and “lu_solve” functions, which don’t return P, L, U (the permutation matrix, the lower triangular portion and the upper triangular portion).  Instead, it returns L and U combined into one matrix to save space and the pivots.  The pivots aren’t used in this example, but the LU form is, along with b to find x.

In the case where b is a measurement vector, we will want to calculate (AtCA)^-1 and simply do a matrix-vector multiply every time b is updated to find a new x and y.  Performing elimination every time b is updated is costly and unnecessary.  A matrix vector multiply can turn out to be faster.