4
$\begingroup$

I've written a C# solver for linear least squares problems with inequality constraints. That is, given $A$, $b$, $G$, $h$

$$\min\|Ax-b\|^2\text{ s.t. }Gx\ge h$$

I have a few hand crafted test problems that my solver gives the correct answer for, and now I'd like to throw a gauntlet of randomly generated problems of various ranks at it to make sure there aren't any edge cases I'm missing.

So what I need is a way to determine that a given $b$ vector calculated satisfies the constraints $Gx \ge h$ (which is easy to check for) and that the solution vector can't be improved by perturbing it in a given non-constrained direction. The second part is what I'm at a loss for.

  • 0
    Do you mean "a given $x$ vector" satisfies $Gx \geq h$?2011-10-18
  • 0
    Yeah. The idea is that you're trying to find $x$ given some "hard" constraints and "soft" suggestions.2011-10-18
  • 0
    Why don't you compare the solution your code gives to one from a mature linear algebra library, e.g. [this one](http://en.wikipedia.org/wiki/Galahad_library)?2012-02-16
  • 0
    I'd like to but I couldn't ever find an actual implementation of a solver for this specific problem. I can't tell if the one you link has one. Of course, even if it does, I have to write the bridge code to let me use the fortran libraries from C#, which is a lot of work :/2012-02-17
  • 0
    How did you solve the problem? What method the solver use?2017-02-21
  • 0
    @Royi: Do you mean how did I test for optimality, or how did I write the solver? If for optimality, I put it on hold and went to other problems. If you mean how did I write the solver, you can take a look at my code here: http://svn.darwinbots.com/Darwinbots3/Trunk/Modules/Azimuth/Azimuth/DenseLinearAlgebra/ConstrainedLeastSquares.cs. There's a few abortive attempts in there that are commented out, but the LSI code at the top is what I'm using. The algorithm is from "Solving least squares problems" by Lawson and Hanson. It transforms the problem in to non negative LS and solves that.2017-02-22
  • 0
    Better than code, where is the theoretical background of your solver? Thank You.2017-02-22
  • 0
    @Royi: It's from "Solving least squares problems" by Lawson and Hanson. ISBN 978-0-898713-56-5. Chapter 23: "Linear least squares with linear inequality constraints". It transforms this problem ("Problem LSI") in to a non-negative least squares problem ("problem NNLS"), for which it provides a pivoting algorithm. NNLS here just means a least squares problem where x >= 0. There are other ways to solve the problem, I'm sure. I think there's a transformation in to linear complimentary problem (LCPs), for instance. But I could never get a robust pivoting solver for those working.2017-02-22
  • 0
    I see. I don't have access to this book and I'm looking for a source which describes how to solve it. I could write the dual form which indeed has a form of OLS with Non Negativity Constrain.2017-02-22
  • 0
    @Royi I don't think you'll find a good source without buying a book. Believe me I've looked. Most of the research on pivoting methods was done in the 70s or earlier, and most of it's not online. In contrast, if you wanted to do an interior method, there's lots of resources online around conjugate gradient and gradient descent.2017-02-22
  • 0
    This is the best source I've found that's free: http://infolab.stanford.edu/pub/cstr/reports/cs/tr/69/134/CS-TR-69-134.pdf2017-02-22
  • 0
    I see. I'm not talking specifically about Pivoting but just how to solve Linear Least Squares with Inequality Constraints.2017-02-23
  • 0
    @Royi: Outside my sphere of knowledge. But maybe this will get you pointed in the right direction: http://www2.esm.vt.edu/~zgurdal/COURSES/4084/4084-Docs/LECTURES/GradProj.pdf2017-02-23
  • 0
    @JayLemmon, by the way, Maybe it can be solved easily using Augmented Lagrangian or Penalty Method.2017-02-23
  • 0
    @JayLemmon, See here: https://ipvs.informatik.uni-stuttgart.de/mlr/marc/teaching/13-Optimization/03-constrainedOpt.pdf2017-02-23
  • 0
    Another nice source: https://www.cs.ubc.ca/~schmidtm/MLSS/constrained.pdf2017-02-23
  • 0
    Note to myself - Add full solution.2017-09-04

2 Answers 2

2

As written above, the nice thing about the KKT conditions is that they allow you to check any solution as it was a black box without a reference solver just by validating its $ \lambda $ on the KKT conditions.

Yet while there are many solvers for this I'm adding another approach - The Projected Gradient Descent.

The Problem Statement

$$ \begin{align*} \arg \min_{x} & \quad & \frac{1}{2} \left\| A x - b \right\|_{2}^{2} \\ \text{subject to} & \quad & C x \leq d \end{align*} $$

The Solution

Now, if we had a closed form projection into the set of the Linear Inequality (Convex Polytop / Convex Polyhedron), which is a Linear Inequality Constraints Least Least Square problem by itself, using the Projected Gradient Descent was easy:

  1. Gradient Descent Step.
  2. Project Solution onto the Inequality Set.
  3. Unless converged go to (1).

Yet, another way is to use Alternating Projections.
Since the set:

$$ \mathcal{S} = \left\{ x \mid C x \leq d \right\} $$

Can be decomposed into:

$$ \mathcal{S} = \cap_{i = 1}^{l} {\mathcal{s}}_{i} = \cap_{i = 1}^{l} \left\{ x \mid {c}_{i}^{T} x \leq {d}_{i} \right\} $$

Where $ {c}_{i} $ is the $ i $ -th row of $ C $.

Projection onto Half Space

In the above, each $ \mathcal{s}_{i} $ is an half space which the projection onto is known:

$$ \begin{align*} \arg \min_{x} & \quad & \frac{1}{2} \left\| x - y \right\|_{2}^{2} \\ \text{subject to} & \quad & {c}^{T} x \leq d \end{align*} $$

The solution is given by:

$$ \hat{x} = y - \lambda c, \quad \lambda = \max \left\{ \frac{ {x}^{T} y - d }{ \left\| c \right\|_{2}^{2} }, 0 \right\} $$

Now the solution becomes:

  1. Gradient Descent Step.
  2. Project Solution onto the Inequality Set:

    • Project onto the 1st Half Space.
    • Project onto the 2nd Half Space.
    • ...
    • Project onto the k-th Half Space.
  3. Unless converged go to (1).

The code for the Gradient Descent is given by:

mAA = mA.' * mA;
vAb = mA.' * vB;

vX          = zeros([numCols, 1]);
vObjVal(1)  = hObjFun(vX);

for ii = 2:numIterations

    vG = (mAA * vX) - vAb; %

The Projection (Alternating Projection):

function [ vX ] = ProjectOntoLinearInequality( vY, mC, vD, stopThr )

numConst    = size(mC, 1);
vCNorm      = sum(mC .^ 2, 2);

vX      = vY;
vRes    = (mC * vX) - vD;
maxRes  = max(vRes);

while(maxRes > stopThr)
    for ii = 1:numConst
        paramLambda = max(((mC(ii, :) * vX) - vD(ii)) / vCNorm(ii), 0);
        vX = vX - (paramLambda * mC(ii, :).');
    end

    vRes = (mC * vX) - vD;
    maxRes = max(vRes);

end


end

The result:

enter image description here

The full code (Including validation against CVX) is available on my StackExchange Math Q73712 GitHub Repository.

1

You want the KKT conditions of the problem; since it is convex, a given $x$ is a minimizer if and only there exists Lagrange multipliers $\lambda$ such that $$\begin{align*}A^TAx - A^T b - G^T\lambda &= 0\\ Gx-h &\geq 0\\ \lambda &\geq 0\\ \lambda^T (Gx-h) &= 0\end{align*}.$$

I'm not completely sure of the best way of finding a $\lambda$ certifying the above or of proving one doesn't exist; you can start by using the last equation to split the constraints into an active set ($G_ax-h_a=0$, $\lambda_a \geq 0$) and inactive set $(G_i x - h_i > 0, \lambda_i = 0)$, and deleting the inactive constraints from the above equations. If $A^TA$ is invertible you can then directly solve for $\lambda_a$ and check if all entries are nonnegative. If $A^TA$ is singular I'm not sure how best to proceed; maybe another answer will elaborate.

  • 0
    That's more or less the way I'm solving the system right now, yeah. I have $\lambda$, and from that I can calculate $x$. I guess I could test those 4 conditions from the KKT directly. But I was hoping for a method (maybe something calculus based) that wasn't so closely coupled with how I solved the problem in the first place.2011-10-18
  • 0
    But when the problem and constraints are convex, as in your problem, satisfying the KKT conditions is proof of optimality -- why look for other criteria (and, of course, these equations *are* Calculus-based; they're essentially "take the derivative and set it equal to zero.")2011-10-18
  • 0
    It's partly numerical issues. The KKT involve a lot of multiplications for things that are supposed to cancel out and hit 0. It's also partly that everything is my own code and I'm trying to find ways to test it that touch as many areas of the code base as possible.2011-10-18
  • 0
    What method are you using to solve the problem? As for testing for optimality, the KKT conditions are the way to go.2011-10-29